Skip to content

Commit

Permalink
Merge branch 'develop' into main
Browse files Browse the repository at this point in the history
  • Loading branch information
dfsp-spirit committed Jan 27, 2021
2 parents 3692f4a + 1716358 commit 1b7c3ed
Show file tree
Hide file tree
Showing 16 changed files with 409 additions and 66 deletions.
17 changes: 12 additions & 5 deletions README.md
Original file line number Diff line number Diff line change
Expand Up @@ -4,7 +4,15 @@ Handling of structural neuroimaging file formats for [Julia](https://julialang.o

## About

Some basic packages for reading neuroimaging data files are available from [JuliaNeuroscience](https://github.com/JuliaNeuroscience), e.g., NIFTI volume and GIFTI mesh support. In this NeuroFormats package, we will provide an API, similar to that of the [freesurferformats R package](https://github.com/dfsp-spirit/freesurferformats), for reading structural neuroimaging data files in Julia. The focus will be on surface-based data, as produced by [FreeSurfer](https://freesurfer.net).
In the NeuroFormats package, we will provide an API, similar to that of the [freesurferformats R package](https://github.com/dfsp-spirit/freesurferformats), for reading structural neuroimaging data files in Julia. The focus will be on surface-based data, as produced by [FreeSurfer](https://freesurfer.net).

Note that some basic packages for reading neuroimaging data files are available from [JuliaNeuroscience](https://github.com/JuliaNeuroscience), e.g., NIFTI volume and GIFTI mesh support.

## Features (WIP)

* Read and write FreeSurfer per-vertex data in curv format (like `subject/surf/lh.thickness`)
* Read brain meshes in FreeSurfer binary mesh format (like `subject/surf/lh.white`)


## Installation

Expand All @@ -16,17 +24,16 @@ cd NeuroFormats.jl
julia --project=.
```

## Usage
## Usage Example

This is a very early package version, please keep in mind that the API is not very stable yet. If you still want to give it a try already, here is what you can do:
Please keep in mind that the API is not very stable yet. If you still want to give it a try already, here is an example for what you can do (after following the installation steps above):

```julia
using NeuroFormats
curv_file = "path/to/subjects_dir/subjectX/surf/lh.thickness" # Adapt path to your data.
thickness = readcurv(curv_file) # An Array{Float32, 1} with your cortical thickness data.
thickness = read_curv(curv_file) # An Array{Float32, 1} with your cortical thickness data.
```


## Continuous integration results:

[![Build Status](https://travis-ci.org/dfsp-spirit/NeuroFormats.jl.svg?branch=main)](https://travis-ci.org/dfsp-spirit/NeuroFormats.jl) Travis CI under Linux
9 changes: 9 additions & 0 deletions examples/export_mesh.jl
Original file line number Diff line number Diff line change
@@ -0,0 +1,9 @@
# Illustrates how to load a FreeSurfer brain mesh and export it to Wavefront Object format.
# You can open the result in Blender or other standard 3D modeling software.

using NeuroFormats

BRAIN_MESH_FILE = joinpath(ENV["HOME"], "develop/NeuroFormats.jl/test/data/subjects_dir/subject1/surf/lh.white")
surface = read_fs_surface(BRAIN_MESH_FILE)
export_to_obj(joinpath(ENV["HOME"], "brain.obj"), surface.mesh)

12 changes: 12 additions & 0 deletions src/FreeSurfer.jl
Original file line number Diff line number Diff line change
@@ -0,0 +1,12 @@
module FreeSurfer

using Printf

export read_fs_int24, interpret_fs_int24, read_curv, write_curv, read_fs_surface, num_vertices, num_faces, export_to_obj

include("./utils.jl")
include("./fs_common.jl")
include("./fs_curv.jl")
include("./fs_surface.jl")

end
5 changes: 4 additions & 1 deletion src/NeuroFormats.jl
Original file line number Diff line number Diff line change
Expand Up @@ -3,7 +3,10 @@ module NeuroFormats

using Reexport

include("./fs_curv.jl")
include("./utils.jl")
export read_variable_length_string

include("./FreeSurfer.jl")
using .FreeSurfer
@reexport using .FreeSurfer

Expand Down
29 changes: 29 additions & 0 deletions src/fs_common.jl
Original file line number Diff line number Diff line change
@@ -0,0 +1,29 @@
# Common functions for reading FreeSurfer data files.

import Base.hton, Base.ntoh


function read_fs_int24(io::IO; endian::AbstractString = "little")
if ! (endian in ["little", "big"])
error("Parameter 'endian' must be one of 'little' or 'big'.")
end

sub_values::Array{Int64,1} = zeros(3)
endian_func = (endian == "little" ? ltoh : ntoh)

b1::Int64 = Int64(endian_func(read(io, UInt8)))
b2::Int64 = Int64(endian_func(read(io, UInt8)))
b3::Int64 = Int64(endian_func(read(io, UInt8)))
fs_int24 = b1 << 16 + b2 << 8 + b3
return fs_int24
end


""" Interpret 3 single-byte unsigned integers as a single integer, as used in several FreeSurfer file formats. """
function interpret_fs_int24(b1::Integer, b2::Integer, b3::Integer)
#@printf("b1=%d, b2=%d, b3=%d (types: %s, %s, %s)\n", b1, b2, b3, typeof(b1), typeof(b2), typeof(b1));
fs_int24::Int64 = Int64(b1) << 16 + Int64(b2) << 8 + Int64(b3)
return fs_int24
end


65 changes: 50 additions & 15 deletions src/fs_curv.jl
Original file line number Diff line number Diff line change
@@ -1,14 +1,9 @@

module FreeSurfer

#using LinearAlgebra
#using Reexport
using Printf

import Base.getindex, Base.size, Base.length, Base.reinterpret, Base.hton
export readcurv, curv_header
import Base.getindex, Base.size, Base.length, Base.reinterpret, Base.hton, Base.ntoh, Base.show


""" Models the header section of a file in Curv format. """
mutable struct CurvHeader
curv_magic_b1::UInt8
curv_magic_b2::UInt8
Expand All @@ -18,17 +13,21 @@ mutable struct CurvHeader
values_per_vertex::Int32
end


""" Models the structure of a file in Curv format. """
mutable struct Curv
header::CurvHeader
data::Array{Float32, 1}
end

Base.show(io::IO, x::Curv) = @printf("FreeSurfer per-vertex brain surface data for %d vertices.\n", Base.length(x.data))

const CURV_MAGIC_HDR = 16777215::Int64
const CURV_HDR_SIZE = sizeof(CurvHeader)


""" Read header from a Curv file """
function readcurv_header(io::IO)
function read_curv_header(io::IO)
header = CurvHeader(UInt8(hton(read(io,UInt8))), UInt8(hton(read(io,UInt8))), UInt8(hton(read(io,UInt8))), hton(read(io,Int32)), hton(read(io,Int32)), hton(read(io,Int32)))
if !(header.curv_magic_b1 == 0xff && header.curv_magic_b2 == 0xff && header.curv_magic_b3 == 0xff)
error("This is not a binary FreeSurfer Curv file: header magic code mismatch.")
Expand All @@ -38,21 +37,22 @@ end


"""
readcurv(file)
read_curv(file::AbstractString; with_header::Bool=false)
Read per-vertex data for brain meshes from the Curv file `file`. The file must be in FreeSurfer binary `Curv` format, like `lh.thickness`.
See also: [`write_curv`](@ref)
# Examples
```julia-repl
julia> curv = readcurv("~/study1/subject1/surf/lh.thickness")
julia> curv = read_curv("~/study1/subject1/surf/lh.thickness")
```
"""
function readcurv(file::AbstractString; with_header::Bool=false)
function read_curv(file::AbstractString; with_header::Bool=false)
file_io = open(file, "r")
header = readcurv_header(file_io)
header = read_curv_header(file_io)

@printf("Loaded curv header with data for %d vertices, now at fh position %d.\n", header.num_vertices, Base.position(file_io))
per_vertex_data = reinterpret(Float32, read(file_io, sizeof(Float32) * header.num_vertices))
per_vertex_data::Array{Float32,1} = reinterpret(Float32, read(file_io, sizeof(Float32) * header.num_vertices))
per_vertex_data .= ntoh.(per_vertex_data)

close(file_io)
Expand All @@ -66,4 +66,39 @@ function readcurv(file::AbstractString; with_header::Bool=false)
end


end # module
"""
write_curv(file::AbstractString, curv_data::Vector{<:Number})
Write a numeric vector to a binary file in FreeSurfer Curv format. The data will be coverted to Float32.
This function is typically used to write surface-based neuroimaging data, like per-vertex cortical thickness
measurements from a reconstructed brain mesh.
See also: [`read_curv`](@ref)
# Examples
```julia-repl
julia> write_curv("~/study1/subject1/surf/lh.thickness", convert(Array{Float32}, zeros(100)))
```
"""
function write_curv(file::AbstractString, curv_data::Vector{<:Number})
curv_data = convert(Vector{Float32}, curv_data)
header = CurvHeader(0xff, 0xff, 0xff, length(curv_data), 0, 1)
file_io = open(file, "w")

# Write header
write(file_io, ntoh(header.curv_magic_b1))
write(file_io, ntoh(header.curv_magic_b2))
write(file_io, ntoh(header.curv_magic_b3))
write(file_io, ntoh(header.num_vertices))
write(file_io, ntoh(header.num_faces))
write(file_io, ntoh(header.values_per_vertex))

# Write data
for idx in eachindex(curv_data)
write(file_io, ntoh(curv_data[idx]))
end

close(file_io)
end

104 changes: 104 additions & 0 deletions src/fs_surface.jl
Original file line number Diff line number Diff line change
@@ -0,0 +1,104 @@
# Functions for reading FreeSurfer brain surface meshes.

import Base.show

const TRIS_MAGIC_FILE_TYPE_NUMBER = 16777214

""" Models the header section of a file in FreeSurfer Surface format. The files are big endian. """
mutable struct FsSurfaceHeader
magic_b1::UInt8
magic_b2::UInt8
magic_b3::UInt8
info_line::AbstractString
num_vertices::Int32
num_faces::Int32
end

num_vertices(sh::FsSurfaceHeader) = sh.num_vertices
num_faces(sh::FsSurfaceHeader) = sh.num_faces


""" Models a trimesh. Vertices are defined by their xyz coordinates, and faces are given as indices into the vertex array. """
struct BrainMesh
vertices::Array{Float32, 2} # vertex xyz coords
faces::Array{Int32, 2} # indices of the 3 vertices forming the face / polygon / triangle
end

num_vertices(bm::BrainMesh) = Base.length(bm.vertices) / 3
num_faces(bm::BrainMesh) = Base.length(bm.faces) / 3
Base.show(io::IO, x::BrainMesh) = @printf("Brain mesh with %d vertices and %d faces.\n", num_vertices(x), num_faces(x))

""" Models FreeSurfer Surface file. """
struct FsSurface
header::FsSurfaceHeader
mesh::BrainMesh
end

num_vertices(fsf::FsSurface) = num_vertices(fsf.mesh)
num_faces(fsf::FsSurface) = num_faces(fsf.mesh)
Base.show(io::IO, x::FsSurface) = @printf("FreeSurfer brain mesh with %d vertices and %d faces.\n", num_vertices(x), num_faces(x))


""" Read header from a FreeSurfer brain surface file """
# For fixed length strings, we could do: my_line = bytestring(readbytes(fh, 4)) I guess.
function read_fs_surface_header(io::IO)
header = FsSurfaceHeader(UInt8(hton(read(io, UInt8))), UInt8(hton(read(io, UInt8))), UInt8(hton(read(io, UInt8))),
read_variable_length_string(io),
Int32(hton(read(io, Int32))), Int32(hton(read(io, Int32))))
magic = interpret_fs_int24(header.magic_b1, header.magic_b2, header.magic_b3)
if magic != TRIS_MAGIC_FILE_TYPE_NUMBER
error("This is not a supported binary FreeSurfer Surface file: header magic code mismatch.")
end
header
end


"""
read_fs_surface(file::AbstractString)
Read a brain surface model represented as a mesh from a file in FreeSurfer binary surface format. Such a file typically represents a single hemisphere.
# Examples
```julia-repl
julia> mesh = read_fs_surface("~/study1/subject1/surf/lh.white")
```
"""
function read_fs_surface(file::AbstractString)
file_io = open(file, "r")
header = read_fs_surface_header(file_io)

vertices_raw::Array{Float32,1} = reinterpret(Float32, read(file_io, sizeof(Float32) * header.num_vertices * 3))
vertices_raw .= ntoh.(vertices_raw)
vertices::Array{Float32,2} = Base.reshape(vertices_raw, (3, Base.length(vertices_raw)÷3))'

faces_raw::Array{Int32,1} = reinterpret(Int32, read(file_io, sizeof(Int32) * header.num_faces * 3))
faces_raw .= ntoh.(faces_raw)
faces::Array{Int32,2} = Base.reshape(faces_raw, (3, Base.length(faces_raw)÷3))'

close(file_io)

surface = FsSurface(header, BrainMesh(vertices, faces))
surface
end


""" Export a brain mesh to a Wavefront Object File. """
function export_to_obj(file:: AbstractString, bm::BrainMesh)
open(file, "w") do f
for row_idx in 1:size(bm.vertices, 1)
row = bm.vertices[row_idx, :]
vertex_rep = @sprintf("v %f %f %f\n", row[1], row[2], row[3])
write(f, vertex_rep)
end
for row_idx in 1:size(bm.faces, 1)
row = bm.faces[row_idx, :]
face_rep = @sprintf("f %d %d %d\n", row[1]+1, row[2]+1, row[3]+1)
write(f, face_rep)
end
end
end


""" Export the mesh of a FreeSurfer surface to a Wavefront Object File. """
export_to_obj(file:: AbstractString, x::FsSurface) = export_to_obj(file, x.mesh)

22 changes: 22 additions & 0 deletions src/utils.jl
Original file line number Diff line number Diff line change
@@ -0,0 +1,22 @@
# General utility functions used throughout the project.


#https://stackoverflow.com/questions/61264545/read-null-terminated-string-from-byte-vector-in-julia


""" Read a variable length, 0-terminated C-style string from a binary file. The trailing zero will be read if consume_zero is set, but not included in the result. """
function read_variable_length_string(io::IO; consume_zero::Bool = false)
res_string = ""
while (! eof(io))
char = read(io, UInt8)
if char == 0
if ! consume_zero
Base.seek(io, Base.position(io) - 1)
end
break
else
res_string = res_string * Char(char)
end
end
res_string
end
Binary file added test/data/subjects_dir/subject1/surf/lh.tinysurface
Binary file not shown.
Binary file added test/data/subjects_dir/subject1/surf/lh.white
Binary file not shown.
Loading

0 comments on commit 1b7c3ed

Please sign in to comment.