Add support for NIfTI2 header

src/NIfTI.jl

@@ -17,26 +17,26 @@ An `N`-dimensional NIfTI volume, with raw data of type
 updated to be consistent with the raw volume. 
 
 # Members
-- `header`: a `NIfTI1Header`
+- `header`: a `NIfTIHeader`
 - `extensions`: a Vector of `NIfTIExtension`s 
 - `raw`: Raw data of type `R`
 """
 struct NIVolume{T<:Number,N,R} <: AbstractArray{T,N}
-    header::NIfTI1Header
+    header::NIfTIHeader
     extensions::Vector{NIfTIExtension}
     raw::R
 end
-NIVolume(header::NIfTI1Header, extensions::Vector{NIfTIExtension}, raw::R) where {R}=
+NIVolume(header::NIfTIHeader, extensions::Vector{NIfTIExtension}, raw::R) where {R} =
     niupdate(new(header, extensions, raw))
 
-NIVolume(header::NIfTI1Header, extensions::Vector{NIfTIExtension}, raw::AbstractArray{T,N}) where {T<:Number,N} =
-    NIVolume{typeof(one(T)*1f0+1f0),N,typeof(raw)}(header, extensions, raw)
-NIVolume(header::NIfTI1Header, raw::AbstractArray{T,N}) where {T<:Number,N} =
-    NIVolume{typeof(one(T)*1f0+1f0),N,typeof(raw)}(header, NIfTIExtension[], raw)
+NIVolume(header::NIfTIHeader, extensions::Vector{NIfTIExtension}, raw::AbstractArray{T,N}) where {T<:Number,N} =
+    NIVolume{typeof(one(T) * 1.0f0 + 1.0f0),N,typeof(raw)}(header, extensions, raw)
+NIVolume(header::NIfTIHeader, raw::AbstractArray{T,N}) where {T<:Number,N} =
+    NIVolume{typeof(one(T) * 1.0f0 + 1.0f0),N,typeof(raw)}(header, NIfTIExtension[], raw)
 
-NIVolume(header::NIfTI1Header, extensions::Vector{NIfTIExtension}, raw::AbstractArray{Bool,N}) where {N} =
+NIVolume(header::NIfTIHeader, extensions::Vector{NIfTIExtension}, raw::AbstractArray{Bool,N}) where {N} =
     NIVolume{Bool,N,typeof(raw)}(header, extensions, raw)
-NIVolume(header::NIfTI1Header, raw::AbstractArray{Bool,N}) where {N} =
+NIVolume(header::NIfTIHeader, raw::AbstractArray{Bool,N}) where {N} =
     NIVolume{Bool,N,typeof(raw)}(header, NIfTIExtension[], raw)
 
 header(x::NIVolume) = getfield(x, :header)
@@ -45,15 +45,15 @@ include("coordinates.jl")
 
 # Always in ms
 """
-    time_step(header::NIfTI1Header)
+    time_step(header::NIfTIHeader)
 
-Get the TR **in ms** from a `NIfTI1Header`.
+Get the TR **in ms** from a `NIfTIHeader`.
 """
-time_step(header::NIfTI1Header) =
-    header.pixdim[5] * TIME_UNIT_MULTIPLIERS[header.xyzt_units >> 3]
+time_step(header::NIfTIHeader) =
+    header.pixdim[5] * TIME_UNIT_MULTIPLIERS[header.xyzt_units>>3]
 
 # Gets the size of a type in bits
-nibitpix(t::Type) = Int16(sizeof(t)*8)
+nibitpix(t::Type) = Int16(sizeof(t) * 8)
 nibitpix(::Type{Bool}) = Int16(1)
 
 # Constructor
@@ -68,36 +68,36 @@ function NIVolume(
     glmin::Integer=Int16(0),
 
     # The frequency encoding, phase encoding and slice dimensions.
-    dim_info::NTuple{3, Integer}=(0, 0, 0),
+    dim_info::NTuple{3,Integer}=(0, 0, 0),
     # Describes data contained in the file; for valid values, see
     # http://nifti.nimh.nih.gov/nifti-1/documentation/nifti1fields/nifti1fields_pages/group__NIfTI1__INTENT__CODES.html
-    intent_p1::Real=0f0, intent_p2::Real=0f0, intent_p3::Real=0f0,
+    intent_p1::Real=0.0f0, intent_p2::Real=0.0f0, intent_p3::Real=0.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=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=UInt8(18),
+    voxel_size::NTuple{3,Real}=(1.0f0, 1.0f0, 1.0f0), time_step::Real=0.0f0, xyzt_units=UInt8(18),
     # Slope and intercept by which volume shoudl be scaled
-    scl_slope::Real=1f0, scl_inter::Real=0f0,
+    scl_slope::Real=1.0f0, scl_inter::Real=0.0f0,
     # These describe how data should be scaled when displayed on the screen. They are probably
     # rarely used
-    cal_max::Real=0f0, cal_min::Real=0f0,
+    cal_max::Real=0.0f0, cal_min::Real=0.0f0,
     # The amount of time taken to acquire a slice
-    slice_duration::Real=0f0,
+    slice_duration::Real=0.0f0,
     # Indicates a non-zero start point for time axis
-    toffset::Real=0f0,
+    toffset::Real=0.0f0,
 
     # "Any text you like"
     descrip::AbstractString="",
     # Name of auxiliary file
     aux_file::AbstractString="",
 
     # Parameters for Method 2. See the NIfTI spec
-    qfac::Float32=0f0, quatern_b::Real=0f0, quatern_c::Real=0f0, quatern_d::Real=0f0,
-    qoffset_x::Real=0f0, qoffset_y::Real=0f0, qoffset_z::Real=0f0,
+    qfac::Float32=0.0f0, quatern_b::Real=0.0f0, quatern_c::Real=0.0f0, quatern_d::Real=0.0f0,
+    qoffset_x::Real=0.0f0, qoffset_y::Real=0.0f0, qoffset_z::Real=0.0f0,
     # Orientation matrix for Method 3
-    orientation::Union{Matrix{Float32},Nothing}=nothing) where {T <: Number}
+    orientation::Union{Matrix{Float32},Nothing}=nothing) where {T<:Number}
     local t
     if isempty(raw)
         raw = Int16[]
@@ -111,7 +111,7 @@ function NIVolume(
     end
 
     method2 = qfac != 0 || quatern_b != 0 || quatern_c != 0 || quatern_d != 0 || qoffset_x != 0 ||
-        qoffset_y != 0 || qoffset_z != 0
+              qoffset_y != 0 || qoffset_z != 0
     method3 = orientation !== nothing
 
     if method2 && method3
@@ -132,29 +132,28 @@ function NIVolume(
     end
 
     NIVolume(NIfTI1Header(SIZEOF_HDR1, string_tuple(data_type, 10), string_tuple(db_name, 18), extents, session_error,
-        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,
-        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, :]...,),
-        (orientation[2, :]...,), (orientation[3, :]...,), string_tuple(intent_name, 16), NP1_MAGIC), extensions, raw)
+            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,
+            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, :]...,),
+            (orientation[2, :]...,), (orientation[3, :]...,), string_tuple(intent_name, 16), NP1_MAGIC), extensions, raw)
 end
 
 # Validates the header of a volume and updates it to match the volume's contents
 function niupdate(vol::NIVolume{T,N,R}) where {T,N,R}
-    vol.header.sizeof_hdr = SIZEOF_HDR1
     vol.header.dim = to_dim_i16(size(vol.raw))
     vol.header.datatype = eltype_to_int16(eltype(R))
     vol.header.bitpix = nibitpix(T)
-    vol.header.vox_offset = isempty(vol.extensions) ? Int32(352) :
-        Int32(mapreduce(esize, +, vol.extensions) + SIZEOF_HDR1)
+    vol.header.vox_offset = isempty(vol.extensions) ? vol.header.sizeof_hdr + 4 :
+                            Int32(mapreduce(esize, +, vol.extensions) + vol.header.sizeof_hdr)
     vol
 end
 
@@ -188,7 +187,7 @@ end
 Write a NIVolume to a file specified by `path`.
 """
 function niwrite(path::AbstractString, vol::NIVolume)
-    if split(path,".")[end] == "gz"
+    if split(path, ".")[end] == "gz"
         io = open(path, "w")
         stream = GzipCompressorStream(io)
         write(stream, vol)
@@ -203,11 +202,19 @@ end
 
 # Read header from a NIfTI file
 function read_header(io::IO)
-    header, swapped = read(io, NIfTI1Header)
-    if header.sizeof_hdr != SIZEOF_HDR1
-        error("This is not a NIfTI-1 file")
+    header_size = peek(io, Int32)
+    if header_size == SIZEOF_HDR1 || bswap(header_size) == SIZEOF_HDR1
+        header, swapped = read(io, NIfTI1Header)
+        #error("This is not a NIfTI-1 file")
+        return header, swapped
     end
-    header, swapped
+
+    if header_size == SIZEOF_HDR2 || bswap(header_size) == SIZEOF_HDR2
+        header, swapped = read(io, NIfTI2Header)
+        #error("This is not a NIfTI-1 file")
+        return header, swapped
+    end
+    error("This is not a NIfTI file")
 end
 
 # Look for a gzip header in an IOStream
@@ -221,7 +228,7 @@ function isgz(io::IO)
             @debug "reading the file resulted in an EOF error and \nthe end of the file was read.\nNo more data was available to read from the filestream.\nIt is likely that the file was corrupted or is empty (0 bytes)."
             rethrow(err)
         end
-    end 
+    end
 end
 
 """
@@ -232,9 +239,9 @@ Read a NIfTI file to a NIVolume. Set `mmap=true` to memory map the volume.
 function niread(file::AbstractString; mmap::Bool=false, mode::AbstractString="r")
     io = niopen(file, mode)
     hdr, swapped = read_header(io)
-    ex = read_extensions(io, hdr.vox_offset - 352)
+    ex = read_extensions(io, hdr.vox_offset - hdr.sizeof_hdr - 4)
 
-    if hdr.magic === NP1_MAGIC
+    if hdr.magic === NP1_MAGIC || hdr.magic == NP2_MAGIC
         vol = read_volume(io, to_eltype(hdr.datatype), to_dimensions(hdr.dim), mmap)
     else
         vol = read_volume(niopen(hdr_to_img(file), mode), to_eltype(hdr.datatype), to_dimensions(hdr.dim), mmap)
@@ -250,8 +257,8 @@ end
 # Allow file to be indexed like an array, but with indices yielding scaled data
 @inline getindex(f::NIVolume{T}, idx::Vararg{Int}) where {T} =
     f.header.scl_slope == zero(typeof(f.header.scl_slope)) ?
-        getindex(f.raw, idx...,) :
-        getindex(f.raw, idx...,) * f.header.scl_slope + f.header.scl_inter
+    getindex(f.raw, idx...,) :
+    getindex(f.raw, idx...,) * f.header.scl_slope + f.header.scl_inter
 
 add1(x::Union{AbstractArray{T},T}) where {T<:Integer} = x + 1
 add1(::Colon) = Colon()

src/coordinates.jl

@@ -19,7 +19,7 @@ end
 end
 
 get_sform(x::NIVolume) = get_sform(x.header)
-function get_sform(hdr::NIfTI1Header)
+function get_sform(hdr::NIfTIHeader)
     if hdr.sform_code > 0
         return @inbounds Float32[
             hdr.srow_x[1]  hdr.srow_x[2]  hdr.srow_x[3]  hdr.srow_x[4]
@@ -34,7 +34,7 @@ function get_sform(hdr::NIfTI1Header)
 end
 
 get_qform(x::NIVolume) = get_qform(x.header)
-function get_qform(hdr::NIfTI1Header)
+function get_qform(hdr::NIfTIHeader)
     if hdr.qform_code <= 0
         return @inbounds Float32[
             hdr.pixdim[2]   0              0              0
@@ -84,11 +84,11 @@ getaffine(x::NIVolume) = getaffine(x.header)
 
 # Convert a NIfTI header to a 4x4 affine transformation matrix
 """
-    getaffine(hdr::NIfTI1Header)
+    getaffine(hdr::NIfTIHeader)
 
 Gets a 4x4 affine transformation matrix from a header's sform
 """
-function getaffine(hdr::NIfTI1Header)
+function getaffine(hdr::NIfTIHeader)
     if hdr.sform_code > 0
         return get_sform(hdr)
     else
@@ -97,11 +97,11 @@ function getaffine(hdr::NIfTI1Header)
 end
 
 """
-    function setaffine(h::NIfTI1Header, affine::Array{T,2}) where {T}
+    function setaffine(h::NIfTIHeader, affine::Array{T,2}) where {T}
 
-Set the affine of a `NIfTI1Header` to 4x4 affine matrix `affine`
+Set the affine of a `NIfTIHeader` to 4x4 affine matrix `affine`
 """
-function setaffine(h::NIfTI1Header, affine::Array{T,2}) where {T}
+function setaffine(h::NIfTIHeader, affine::Array{T,2}) where {T}
     size(affine, 1) == size(affine, 2) == 4 ||
         error("affine matrix must be 4x4")
     affine[4, 1] == affine[4, 2] == affine[4, 3] == 0 && affine[4, 4] == 1 ||
@@ -127,7 +127,7 @@ end
 Returns a tuple providing the orientation of a NIfTI image.
 """
 orientation(x) = orientation(x.header)
-function orientation(hdr::NIfTI1Header)
+function orientation(hdr::NIfTIHeader)
     if hdr.sform_code > 0
         return @inbounds _dir2ori(
             hdr.srow_x[1], hdr.srow_x[2], hdr.srow_x[3],