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)

2 comments on commit 06d4ce3

@Tokazama
Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@JuliaRegistrator
Copy link

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Registration pull request created: JuliaRegistries/General/57225

After the above pull request is merged, it is recommended that a tag is created on this repository for the registered package version.

This will be done automatically if the Julia TagBot GitHub Action is installed, or can be done manually through the github interface, or via:

git tag -a v0.5.7 -m "<description of version>" 06d4ce372ffc1a3ce31953ecad8273afa6a6c868
git push origin v0.5.7

Please sign in to comment.