Skip to content

Commit

Permalink
phasdim, freqdim, slicedim support
Browse files Browse the repository at this point in the history
  • Loading branch information
Tokazama committed Mar 24, 2022
1 parent 76c84c8 commit 06d4ce3
Show file tree
Hide file tree
Showing 4 changed files with 171 additions and 111 deletions.
2 changes: 1 addition & 1 deletion Project.toml
Original file line number Diff line number Diff line change
@@ -1,6 +1,6 @@
name = "NIfTI"
uuid = "a3a9e032-41b5-5fc4-967a-a6b7a19844d3"
version = "0.5.6"
version = "0.5.7"

[deps]
Base64 = "2a0f44e3-6c83-55bd-87e4-b1978d98bd5f"
Expand Down
102 changes: 3 additions & 99 deletions src/NIfTI.jl
Original file line number Diff line number Diff line change
@@ -1,6 +1,3 @@
# NIfTI.jl
# Methods for reading NIfTI MRI files in Julia

module NIfTI

using CodecZlib, Mmap, MappedArrays, TranscodingStreams
Expand All @@ -11,101 +8,7 @@ export NIVolume, niread, niwrite, voxel_size, time_step, vox, getaffine, setaffi
include("parsers.jl")
include("extensions.jl")
include("volume.jl")

function define_packed(ty::DataType)
packed_offsets = cumsum([sizeof(x) for x in ty.types])
sz = pop!(packed_offsets)
pushfirst!(packed_offsets, 0)

@eval begin
function Base.read(io::IO, ::Type{$ty})
bytes = read!(io, Array{UInt8}(undef, $sz...))
hdr = $(Expr(:new, ty, [:(unsafe_load(convert(Ptr{$(ty.types[i])}, pointer(bytes)+$(packed_offsets[i])))) for i = 1:length(packed_offsets)]...,))
if hdr.sizeof_hdr == ntoh(Int32(348))
return byteswap(hdr), true
end
hdr, false
end
function Base.write(io::IO, x::$ty)
bytes = UInt8[]
for name in fieldnames($ty)
append!(bytes, reinterpret(UInt8, [getfield(x,name)]))
end
write(io, bytes)
$sz
end
end
nothing
end

mutable struct NIfTI1Header
sizeof_hdr::Int32

data_type::NTuple{10,UInt8}
db_name::NTuple{18,UInt8}
extents::Int32
session_error::Int16
regular::Int8

dim_info::Int8
dim::NTuple{8,Int16}
intent_p1::Float32
intent_p2::Float32
intent_p3::Float32
intent_code::Int16
datatype::Int16
bitpix::Int16
slice_start::Int16
pixdim::NTuple{8,Float32}
vox_offset::Float32
scl_slope::Float32
scl_inter::Float32
slice_end::Int16
slice_code::Int8
xyzt_units::Int8
cal_max::Float32
cal_min::Float32
slice_duration::Float32
toffset::Float32

glmax::Int32
glmin::Int32

descrip::NTuple{80,UInt8}
aux_file::NTuple{24,UInt8}

qform_code::Int16
sform_code::Int16
quatern_b::Float32
quatern_c::Float32
quatern_d::Float32
qoffset_x::Float32
qoffset_y::Float32
qoffset_z::Float32

srow_x::NTuple{4,Float32}
srow_y::NTuple{4,Float32}
srow_z::NTuple{4,Float32}

intent_name::NTuple{16,UInt8}

magic::NTuple{4,UInt8}
end
define_packed(NIfTI1Header)

# byteswapping

function byteswap(hdr::NIfTI1Header)
for fn in fieldnames(typeof(hdr))
val = getfield(hdr, fn)
if isa(val, Number) && sizeof(val) > 1
setfield!(hdr, fn, ntoh(val))
elseif isa(val, NTuple) && sizeof(eltype(val)) > 1
setfield!(hdr, fn, map(ntoh, val))
end
end
hdr
end
include("headers.jl")

"""
NIVolume{T<:Number,N,R} <: AbstractArray{T,N}
Expand Down Expand Up @@ -136,6 +39,8 @@ NIVolume(header::NIfTI1Header, raw::AbstractArray{Bool,N}) where {N} =
NIVolume{Bool,N,typeof(raw)}(header, NIfTIExtension[], raw)


header(x::NIVolume) = getfield(x, :header)

include("coordinates.jl")


Expand Down Expand Up @@ -393,4 +298,3 @@ length(f::NIVolume) = length(f.raw)
lastindex(f::NIVolume) = lastindex(f.raw)

end

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


function define_packed(ty::DataType)
packed_offsets = cumsum([sizeof(x) for x in ty.types])
sz = pop!(packed_offsets)
pushfirst!(packed_offsets, 0)

@eval begin
function Base.read(io::IO, ::Type{$ty})
bytes = read!(io, Array{UInt8}(undef, $sz...))
hdr = $(Expr(:new, ty, [:(unsafe_load(convert(Ptr{$(ty.types[i])}, pointer(bytes)+$(packed_offsets[i])))) for i = 1:length(packed_offsets)]...,))
if hdr.sizeof_hdr == ntoh(Int32(348))
return byteswap(hdr), true
end
hdr, false
end
function Base.write(io::IO, x::$ty)
bytes = UInt8[]
for name in fieldnames($ty)
append!(bytes, reinterpret(UInt8, [getfield(x,name)]))
end
write(io, bytes)
$sz
end
end
nothing
end

mutable struct NIfTI1Header
sizeof_hdr::Int32

data_type::NTuple{10,UInt8}
db_name::NTuple{18,UInt8}
extents::Int32
session_error::Int16
regular::Int8

dim_info::Int8
dim::NTuple{8,Int16}
intent_p1::Float32
intent_p2::Float32
intent_p3::Float32
intent_code::Int16
datatype::Int16
bitpix::Int16
slice_start::Int16
pixdim::NTuple{8,Float32}
vox_offset::Float32
scl_slope::Float32
scl_inter::Float32
slice_end::Int16
slice_code::Int8
xyzt_units::Int8
cal_max::Float32
cal_min::Float32
slice_duration::Float32
toffset::Float32

glmax::Int32
glmin::Int32

descrip::NTuple{80,UInt8}
aux_file::NTuple{24,UInt8}

qform_code::Int16
sform_code::Int16
quatern_b::Float32
quatern_c::Float32
quatern_d::Float32
qoffset_x::Float32
qoffset_y::Float32
qoffset_z::Float32

srow_x::NTuple{4,Float32}
srow_y::NTuple{4,Float32}
srow_z::NTuple{4,Float32}

intent_name::NTuple{16,UInt8}

magic::NTuple{4,UInt8}
end
define_packed(NIfTI1Header)

# byteswapping

function byteswap(hdr::NIfTI1Header)
for fn in fieldnames(typeof(hdr))
val = getfield(hdr, fn)
if isa(val, Number) && sizeof(val) > 1
setfield!(hdr, fn, ntoh(val))
elseif isa(val, NTuple) && sizeof(eltype(val)) > 1
setfield!(hdr, fn, map(ntoh, val))
end
end
hdr
end

"""
freqdim(x)::Int
Returns the frequency dimension.
"""
freqdim(x::NIfTI1Header) = Int(getfield(x, :dim_info) & 0x03)
freqdim(x) = freqdim(header(x))

"""
phasedim(x)::Int
Returns the phase dimension.
"""
phasedim(x::NIfTI1Header) = Int((getfield(x, :dim_info) >> 2) & 0x03)
phasedim(x) = phasedim(header(x))

"""
slicedim(x)::Int
Returns the slice dimension.
"""
slicedim(x::NIfTI1Header) = Int((getfield(x, :dim_info) >> 4) & 0x03)
slicedim(x) = slicedim(header(x))

"""
slice_start(x)::Int
Which slice corresponds to the first slice acquired during MRI acquisition (i.e. not padded slices).
"""
slice_start(x::NIfTI1Header) = Int(getfield(x, :slice_start)) + 1
slice_start(x) = slice_start(header(x))

"""
slice_end(x)::Int
Which slice corresponds to the last slice acquired during MRI acquisition (i.e. not padded slices).
"""
slice_end(x::NIfTI1Header) = Int(getfield(x, :slice_end)) + 1
slice_end(x) = slice_end(header(x))

"""
slice_duration(x)
Time to acquire one slice
"""
slice_duration(x::NIfTI1Header) = getfield(x, :slice_duration)
slice_duration(x) = slice_duration(header(x))

33 changes: 22 additions & 11 deletions test/runtests.jl
Original file line number Diff line number Diff line change
Expand Up @@ -33,34 +33,34 @@ extractto(GZIPPED_HDR, HDR)
extractto(joinpath(dirname(@__FILE__), "data/example4d.img.gz"), IMG)

function image_tests(fname, mmap)
file = niread(fname, mmap=mmap)
img = niread(fname, mmap=mmap)

# Header
@test time_step(file.header) == 2000000 # Actually an error in the file AFAIK
@test voxel_size(file.header) Float32[2.0, 2.0, 2.2]
@test size(file) == (128, 96, 24, 2)
@test time_step(img.header) == 2000000 # Actually an error in the file AFAIK
@test voxel_size(img.header) Float32[2.0, 2.0, 2.2]
@test size(img) == (128, 96, 24, 2)

# Content
@test file.raw[65, 49, 13, :][:] == [265, 266]
@test vox(file, 64, 48, 12, :)[:] == [265, 266]
@test vox(file, 69, 56, 13, :)[:] == [502, 521]
@test img.raw[65, 49, 13, :][:] == [265, 266]
@test vox(img, 64, 48, 12, :)[:] == [265, 266]
@test vox(img, 69, 56, 13, :)[:] == [502, 521]

@assert maximum(file) == maximum(file.raw)
@assert maximum(img) == maximum(img.raw)

@test getaffine(file) [
@test getaffine(img) [
-2.0 6.714715653593746e-19 9.081024511081715e-18 117.8551025390625
6.714715653593746e-19 1.9737114906311035 -0.35552823543548584 -35.72294235229492
8.25548088896093e-18 0.3232076168060303 2.171081781387329 -7.248798370361328
0.0 0.0 0.0 1.0
]

@test NIfTI.get_qform(file) [
@test NIfTI.get_qform(img) [
-2.0 7.75482f-26 -6.93824f-27 117.855
7.75482f-26 1.97371 -0.355528 -35.7229
6.30749f-27 0.323208 2.17108 -7.2488
0.0 0.0 0.0 1.0
]
@test NIfTI.orientation(file) == (:right, :posterior, :inferior)
@test NIfTI.orientation(img) == (:right, :posterior, :inferior)
end

image_tests(NII, false)
Expand All @@ -70,6 +70,16 @@ image_tests(HDR, true)
image_tests(GZIPPED_NII, false)
image_tests(GZIPPED_HDR, false)

@testset "Header Field Accessors" begin
img = niread(GZIPPED_NII)
@test NIfTI.freqdim(img) == 1
@test NIfTI.phasedim(img) == 2
@test NIfTI.slicedim(img) == 3
@test NIfTI.slice_start(img) == 1
@test NIfTI.slice_end(img) == 24
@test NIfTI.slice_duration(img) == 0
end

@test_throws ErrorException niread(GZIPPED_NII; mmap=true)
@test_throws ErrorException niread(GZIPPED_HDR; mmap=true)

Expand Down Expand Up @@ -124,3 +134,4 @@ GC.gc() # closes mmapped files
@test NIfTI._dir2ori(1.0, 0.0, 0.0,
0.0, -1.0, 0.0,
0.0, 0.0, -1.0) == (:left, :anterior, :superior)

0 comments on commit 06d4ce3

Please sign in to comment.