Add explicit doc strings and use better internal dimension queries

Author: Tokazama

docs/src/index.md

@@ -3,9 +3,17 @@
 CurrentModule = NIfTI
 ```
 
-# NeuroGraphs
+# NIfTI
 
 Documentation for [NIfTI](https://github.com/JuliaNeuroscience/NIfTI.jl).
 
 ```@docs
+NIfTI.freqdim
+NIfTI.phasedim
+NIfTI.slicedim
+NIfTI.slice_start
+NIfTI.slice_end
+NIfTI.slice_duration
+NIfTI.sdims
+NIfTI.voxel_size
 ```

src/NIfTI.jl

@@ -43,15 +43,6 @@ header(x::NIVolume) = getfield(x, :header)
 
 include("coordinates.jl")
 
-
-"""
-    voxel_size(header::NIfTI1Header)
-
-Get the voxel size **in mm** from a `NIfTI1Header`.
-"""
-voxel_size(header::NIfTI1Header) =
-    [header.pixdim[i] * SPATIAL_UNIT_MULTIPLIERS[header.xyzt_units & Int8(3)] for i = 2:min(header.dim[1], 3)+1]
-
 # Always in ms
 """
     time_step(header::NIfTI1Header)
@@ -61,30 +52,6 @@ Get the TR **in ms** from a `NIfTI1Header`.
 time_step(header::NIfTI1Header) =
     header.pixdim[5] * TIME_UNIT_MULTIPLIERS[header.xyzt_units >> 3]
 
-function to_dim_info(dim_info::Tuple{Integer,Integer,Integer})
-    if dim_info[1] > 3 || dim_info[1] < 0
-        error("Invalid frequency dimension $(dim_info[1])")
-    elseif dim_info[2] > 3 || dim_info[2] < 0
-        error("Invalid phase dimension $(dim_info[2])")
-    elseif dim_info[3] > 3 || dim_info[3] < 0
-        error("Invalid slice dimension $(dim_info[3])")
-    end
-
-    return Int8(dim_info[1] | (dim_info[2] << 2) | (dim_info[3] << 4))
-end
-
-# Returns or sets dim_info as a tuple whose values are the frequency, phase, and slice dimensions
-function dim_info(header::NIfTI1Header)
-    return (
-        header.dim_info & int8(3),
-        (header.dim_info >> 2) & int8(3),
-        (header.dim_info >> 4) & int8(3)
-    )
-end
-function dim_info(header::NIfTI1Header, dim_info::Tuple{T, T, T}) where {T<:Integer}
-    header.dim_info = to_dim_info(dim_info)
-end
-
 # Gets the size of a type in bits
 nibitpix(t::Type) = Int16(sizeof(t)*8)
 nibitpix(::Type{Bool}) = Int16(1)
@@ -107,10 +74,10 @@ function NIVolume(
     intent_p1::Real=0f0, intent_p2::Real=0f0, intent_p3::Real=0f0,
     intent_code::Integer=Int16(0), intent_name::AbstractString="",
     # Information about which slices were acquired, in case the volume has been padded
-    slice_start::Integer=Int16(0), slice_end::Integer=Int16(0), slice_code::Int8=Int8(0),
+    slice_start::Integer=Int16(0), slice_end::Integer=Int16(0), slice_code=UInt8(0),
     # The size of each voxel and the time step. These are formulated in mm unless xyzt_units is
     # explicitly specified
-    voxel_size::NTuple{3, Real}=(1f0, 1f0, 1f0), time_step::Real=0f0, xyzt_units::Int8=Int8(18),
+    voxel_size::NTuple{3, Real}=(1f0, 1f0, 1f0), time_step::Real=0f0, xyzt_units=UInt8(18),
     # Slope and intercept by which volume shoudl be scaled
     scl_slope::Real=1f0, scl_inter::Real=0f0,
     # These describe how data should be scaled when displayed on the screen. They are probably
@@ -168,8 +135,12 @@ function NIVolume(
         regular, to_dim_info(dim_info), to_dim_i16(size(raw)), intent_p1, intent_p2,
         intent_p3, intent_code, eltype_to_int16(t), nibitpix(t),
         slice_start, (qfac, voxel_size..., time_step, 0, 0, 0), 352,
-        scl_slope, scl_inter, slice_end, slice_code,
-        xyzt_units, cal_max, cal_min, slice_duration,
+        scl_slope,
+        scl_inter,
+        slice_end,
+        UInt8(slice_code),
+        UInt8(xyzt_units),
+        cal_max, cal_min, slice_duration,
         toffset, glmax, glmin, string_tuple(descrip, 80), string_tuple(aux_file, 24), (method2 || method3),
         method3, quatern_b, quatern_c, quatern_d,
         qoffset_x, qoffset_y, qoffset_z, (orientation[1, :]...,),

src/headers.jl

@@ -28,13 +28,13 @@ end
 mutable struct NIfTI1Header
     sizeof_hdr::Int32
 
-    data_type::NTuple{10,UInt8}
-    db_name::NTuple{18,UInt8}
+    data_type::NTuple{10, UInt8}
+    db_name::NTuple{18, UInt8}
     extents::Int32
     session_error::Int16
     regular::Int8
 
-    dim_info::Int8
+    dim_info::UInt8
     dim::NTuple{8,Int16}
     intent_p1::Float32
     intent_p2::Float32
@@ -48,8 +48,8 @@ mutable struct NIfTI1Header
     scl_slope::Float32
     scl_inter::Float32
     slice_end::Int16
-    slice_code::Int8
-    xyzt_units::Int8
+    slice_code::UInt8
+    xyzt_units::UInt8
     cal_max::Float32
     cal_min::Float32
     slice_duration::Float32
@@ -58,8 +58,8 @@ mutable struct NIfTI1Header
     glmax::Int32
     glmin::Int32
 
-    descrip::NTuple{80,UInt8}
-    aux_file::NTuple{24,UInt8}
+    descrip::NTuple{80, UInt8}
+    aux_file::NTuple{24, UInt8}
 
     qform_code::Int16
     sform_code::Int16
@@ -70,13 +70,13 @@ mutable struct NIfTI1Header
     qoffset_y::Float32
     qoffset_z::Float32
 
-    srow_x::NTuple{4,Float32}
-    srow_y::NTuple{4,Float32}
-    srow_z::NTuple{4,Float32}
+    srow_x::NTuple{4, Float32}
+    srow_y::NTuple{4, Float32}
+    srow_z::NTuple{4, Float32}
 
-    intent_name::NTuple{16,UInt8}
+    intent_name::NTuple{16, UInt8}
 
-    magic::NTuple{4,UInt8}
+    magic::NTuple{4, UInt8}
 end
 define_packed(NIfTI1Header)
 
@@ -95,53 +95,125 @@ function byteswap(hdr::NIfTI1Header)
 end
 
 """
-    freqdim(img)::Int
+    NIfTI.slice_start(x)::Int
 
-Returns the frequency dimension associated with with `img`. `img` is an image or a collection of
-image related metadata. 
+Which slice corresponds to the first slice acquired during MRI acquisition (i.e. not padded slices).
 """
-freqdim(x::NIfTI1Header) = Int(getfield(x, :dim_info) & 0x03)
-freqdim(x) = freqdim(header(x))
+slice_start(x::NIfTI1Header) = Int(getfield(x, :slice_start)) + 1
+slice_start(x) = slice_start(header(x))
 
 """
-    phasedim(img)::Int
+    NIfTI.slice_end(x)::Int
 
-Returns the phase dimension associated with `img`. `img` is an image or a collection of
-image related metadata.
+Which slice corresponds to the last slice acquired during MRI acquisition (i.e. not padded slices).
 """
-phasedim(x::NIfTI1Header) = Int((getfield(x, :dim_info) >> 2) & 0x03)
-phasedim(x) = phasedim(header(x))
+slice_end(x::NIfTI1Header) = Int(getfield(x, :slice_end)) + 1
+slice_end(x) = slice_end(header(x))
 
 """
-    slicedim(img)::Int
+    NIfTI.slice_duration(x)
 
-Returns the slice dimension associated with `img`. `img` is an image or a collection of
-image related metadata.
+Time to acquire one slice.
 """
-slicedim(x::NIfTI1Header) = Int((getfield(x, :dim_info) >> 4) & 0x03)
-slicedim(x) = slicedim(header(x))
+slice_duration(x::NIfTI1Header) = getfield(x, :slice_duration)
+slice_duration(x) = slice_duration(header(x))
 
 """
-    slice_start(x)::Int
+    NIfTI.sdims(img)
 
-Which slice corresponds to the first slice acquired during MRI acquisition (i.e. not padded slices).
+Return the number of spatial dimensions in the image.
 """
-slice_start(x::NIfTI1Header) = Int(getfield(x, :slice_start)) + 1
-slice_start(x) = slice_start(header(x))
+sdims(x::NIfTI1Header) = min(Int(getfield(getfield(x, :dim), 1)), 3)
+sdims(x) = sdims(header(x))
+
+# Conversion factors to mm/ms
+# http://nifti.nimh.nih.gov/nifti-1/documentation/nifti1fields/nifti1fields_pages/xyzt_units.html
+# 1 => NIfTI_UNITS_METER
+# 2 => NIfTI_UNITS_MM
+# 3 => NIfTI_UNITS_MICRON
+function to_spatial_multiplier(xyzt_units::UInt8)
+    i = xyzt_units & 0x03
+    if i === 0x01
+        return 1000.0f0
+    elseif i === 0x02
+        return 1.0f0
+    else
+        return 0.001f0
+    end
+end
 
 """
-    slice_end(x)::Int
+    NIfTI.voxel_size(header::NIfTI1Header)
 
-Which slice corresponds to the last slice acquired during MRI acquisition (i.e. not padded slices).
+Get the voxel size **in mm** from a `NIfTI1Header`.
 """
-slice_end(x::NIfTI1Header) = Int(getfield(x, :slice_end)) + 1
-slice_end(x) = slice_end(header(x))
+function voxel_size(x::NIfTI1Header)
+    sd = sdims(x)
+    if sd === 0
+        return ()
+    else
+        sconvert = to_spatial_multiplier(getfield(x, :xyzt_units))
+        if sd === 3
+            pd = getfield(x, :pixdim)
+            return (
+                getfield(pd, 2) * sconvert,
+                getfield(pd, 3) * sconvert,
+                getfield(pd, 4) * sconvert
+            )
+        elseif sd === 2
+            pd = getfield(x, :pixdim)
+            return (getfield(pd, 2) * sconvert, getfield(pd, 3) * sconvert,)
+        else # sd === 1
+            pd = getfield(x, :pixdim)
+            return (getfield(pd, 2) * sconvert,)
+        end
+    end
+end
+
+function to_dim_info(dim_info::Tuple{Integer,Integer,Integer})
+    if dim_info[1] > 3 || dim_info[1] < 0
+        error("Invalid frequency dimension $(dim_info[1])")
+    elseif dim_info[2] > 3 || dim_info[2] < 0
+        error("Invalid phase dimension $(dim_info[2])")
+    elseif dim_info[3] > 3 || dim_info[3] < 0
+        error("Invalid slice dimension $(dim_info[3])")
+    end
+
+    return UInt8(dim_info[1] | (dim_info[2] << 2) | (dim_info[3] << 4))
+end
+
+# Returns or sets dim_info as a tuple whose values are the frequency, phase, and slice dimensions
+function dim_info(hdr::NIfTI1Header)
+    return (hdr.dim_info & 0x03, (hdr.dim_info >> 0x02) & 0x03, (hdr.dim_info >> 0x04) & 0x03)
+end
+function dim_info(header::NIfTI1Header, dim_info::Tuple{T, T, T}) where {T<:Integer}
+    header.dim_info = to_dim_info(dim_info)
+end
 
 """
-    slice_duration(x)
+    NIfTI.freqdim(img)::Int
 
-Time to acquire one slice.
+Returns the frequency dimension associated with with `img`. `img` is an image or a collection of
+image related metadata. 
 """
-slice_duration(x::NIfTI1Header) = getfield(x, :slice_duration)
-slice_duration(x) = slice_duration(header(x))
+freqdim(x::NIfTI1Header) = Int(getfield(x, :dim_info) & 0x03)
+freqdim(x) = freqdim(header(x))
+
+"""
+    NIfTI.phasedim(img)::Int
+
+Returns the phase dimension associated with `img`. `img` is an image or a collection of
+image related metadata.
+"""
+phasedim(x::NIfTI1Header) = Int((getfield(x, :dim_info) >> 0x02) & 0x03)
+phasedim(x) = phasedim(header(x))
+
+"""
+    NIfTI.slicedim(img)::Int
+
+Returns the slice dimension associated with `img`. `img` is an image or a collection of
+image related metadata.
+"""
+slicedim(x::NIfTI1Header) = Int((getfield(x, :dim_info) >> 0x04) & 0x03)
+slicedim(x) = slicedim(header(x))
 

src/parsers.jl

@@ -1,4 +1,8 @@
-
+const SPATIAL_UNIT_MULTIPLIERS = [
+    1000,   # 1 => NIfTI_UNITS_METER
+    1,      # 2 => NIfTI_UNITS_MM
+    0.001   # 3 => NIfTI_UNITS_MICRON
+]
 
 const SIZEOF_HDR1 = Int32(348)
 const SIZEOF_HDR2 = Int32(540)
@@ -200,13 +204,6 @@ function string_tuple(x::String, n::Int)
 end
 string_tuple(x::AbstractString) = string_tuple(bytestring(x))
 
-# Conversion factors to mm/ms
-# http://nifti.nimh.nih.gov/nifti-1/documentation/nifti1fields/nifti1fields_pages/xyzt_units.html
-const SPATIAL_UNIT_MULTIPLIERS = [
-    1000,   # 1 => NIfTI_UNITS_METER
-    1,      # 2 => NIfTI_UNITS_MM
-    0.001   # 3 => NIfTI_UNITS_MICRON
-]
 const TIME_UNIT_MULTIPLIERS = [
     1000,   # NIfTI_UNITS_SEC
     1,      # NIfTI_UNITS_MSEC

test/runtests.jl

@@ -37,7 +37,11 @@ function image_tests(fname, mmap)
 
     # Header
     @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 all(isapprox.(voxel_size(img.header), (2.0, 2.0, 2.2)))
+    vs1, vs2, vs3 = voxel_size(img.header)
+    @test isapprox(vs1, Float32(2.0))
+    @test isapprox(vs2, Float32(2.0))
+    @test isapprox(vs3, Float32(2.2))
     @test size(img) == (128, 96, 24, 2)
 
     # Content