Merge pull request #19 from yhuang43/master

support mgz warp

surfa/core/framed.py

@@ -5,6 +5,18 @@
 from surfa.core.labels import LabelLookup
 
 
+# mgz now has its intent encoded in the version number
+# version = (intent & 0xff ) << 8) | MGH_VERSION
+# MGH_VERSION = 1
+class FramedArrayIntents:
+    unknown     = -1
+    mri         = 0
+    label       = 1
+    shape       = 2
+    warpmap     = 3
+    warpmap_inv = 4
+
+
 class FramedArray:
 
     def __init__(self, basedim, data, labels=None, metadata=None):
@@ -264,7 +276,8 @@ def _shape_changed(self):
         """
         pass
 
-    def save(self, filename, fmt=None):
+    # optional parameter to specify FramedArray intent, default is MRI data
+    def save(self, filename, fmt=None, intent=FramedArrayIntents.mri):
         """
         Write array to file.
 
@@ -276,7 +289,7 @@ def save(self, filename, fmt=None):
             Optional file format to force.
         """
         from surfa.io.framed import save_framed_array
-        save_framed_array(self, filename, fmt=fmt)
+        save_framed_array(self, filename, fmt=fmt, intent=intent)
 
     def min(self, nonzero=False, frames=False):
         """

surfa/io/framed.py

@@ -9,13 +9,16 @@
 from surfa import ImageGeometry
 from surfa.core.array import pad_vector_length
 from surfa.core.framed import FramedArray
+from surfa.core.framed import FramedArrayIntents
 from surfa.image.framed import FramedImage
 from surfa.io import fsio
 from surfa.io import protocol
 from surfa.io.utils import read_int
 from surfa.io.utils import write_int
 from surfa.io.utils import read_bytes
 from surfa.io.utils import write_bytes
+from surfa.io.utils import read_volgeom
+from surfa.io.utils import write_volgeom
 from surfa.io.utils import check_file_readability
 
 
@@ -117,7 +120,8 @@ def load_framed_array(filename, atype, fmt=None):
     return iop().load(filename, atype)
 
 
-def save_framed_array(arr, filename, fmt=None):
+# optional parameter to specify FramedArray intent, default is MRI data
+def save_framed_array(arr, filename, fmt=None, intent=FramedArrayIntents.mri):
     """
     Save a `FramedArray` object to file.
 
@@ -141,7 +145,11 @@ def save_framed_array(arr, filename, fmt=None):
             raise ValueError(f'unknown file format {fmt}')
         filename = iop.enforce_extension(filename)
 
-    iop().save(arr, filename)
+    # pass intent if iop() is an instance of MGHArrayIO
+    if (isinstance(iop(), MGHArrayIO)):
+       iop().save(arr, filename, intent=intent)
+    else:
+        iop().save(arr, filename)
 
 
 def framed_array_from_4d(atype, data):
@@ -231,8 +239,8 @@ def load(self, filename, atype):
         fopen = gzip.open if filename.lower().endswith('gz') else open
         with fopen(filename, 'rb') as file:
 
-            # skip version tag
-            file.read(4)
+            # read version number, retrieve intent
+            intent = read_bytes(file, '>i4', 1) >> 8 & 0xff
 
             # read shape and type info
             shape = read_bytes(file, '>u4', 4)
@@ -283,6 +291,7 @@ def load(self, filename, atype):
             if isinstance(arr, FramedImage):
                 arr.geom.update(**geom_params)
                 arr.metadata.update(scan_params)
+                arr.metadata['intent'] = intent
 
             # read metadata tags
             while True:
@@ -312,13 +321,28 @@ def load(self, filename, atype):
                 elif tag == fsio.tags.fieldstrength:
                     arr.metadata['field-strength'] = read_bytes(file, dtype='>f4')
 
+                # gcamorph src & trg geoms (mgz warp)
+                elif tag == fsio.tags.gcamorph_geom:
+                    # read src vol geom
+                    arr.metadata['gcamorph_volgeom_src'] = read_volgeom(file)
+
+                    # read trg vol geom
+                    arr.metadata['gcamorph_volgeom_trg'] = read_volgeom(file)
+
+                # gcamorph meta (mgz warp: int int float)
+                elif tag == fsio.tags.gcamorph_meta:
+                    arr.metadata['warpfield_dtfmt']  = read_bytes(file, dtype='>i4')
+                    arr.metadata['gcamorph_spacing'] = read_bytes(file, dtype='>i4')
+                    arr.metadata['gcamorph_exp_k']   = read_bytes(file, dtype='>f4')
+
                 # skip everything else
                 else:
                     file.read(length)
 
         return arr
 
-    def save(self, arr, filename):
+    # optional parameter to specify FramedArray intent, default is MRI data
+    def save(self, arr, filename, intent=FramedArrayIntents.mri):
         """
         Write array to a MGH/MGZ file.
 
@@ -364,7 +388,8 @@ def save(self, arr, filename):
             shape[-1] = arr.nframes
 
             # begin writing header
-            write_bytes(file, 1, '>u4')  # version
+            version = ((intent & 0xff) << 8) | 1  # encode intent in version
+            write_bytes(file, version, '>u4')  # version
             write_bytes(file, shape, '>u4')  # shape
             write_bytes(file, dtype_id, '>u4')  # MGH data type
             write_bytes(file, 1, '>u4')  # DOF
@@ -412,6 +437,19 @@ def save(self, arr, filename):
             fsio.write_tag(file, fsio.tags.fieldstrength, 4)
             write_bytes(file, arr.metadata.get('field-strength', 0.0), '>f4')
 
+            # gcamorph geom and gcamorph meta for mgz warp
+            if (intent == FramedArrayIntents.warpmap):
+                # gcamorph src & trg geoms (mgz warp)
+                fsio.write_tag(file, fsio.tags.gcamorph_geom)
+                write_volgeom(file, arr.metadata['gcamorph_volgeom_src'])
+                write_volgeom(file, arr.metadata['gcamorph_volgeom_trg'])
+
+                # gcamorph meta (mgz warp: int int float)
+                fsio.write_tag(file, fsio.tags.gcamorph_meta, 12)
+                write_bytes(file, arr.metadata.get('warpfield_dtfmt', 0),    dtype='>i4')
+                write_bytes(file, arr.metadata.get('gcamorph_spacing', 0.0), dtype='>i4')
+                write_bytes(file, arr.metadata.get('gcamorph_exp_k', 0.0),   dtype='>f4')
+
             # write history tags
             for hist in arr.metadata.get('history', []):
                 fsio.write_tag(file, fsio.tags.history, len(hist))

surfa/io/fsio.py

@@ -17,8 +17,13 @@ class tags:
     gcamorph_geom = 10
     gcamorph_type = 11
     gcamorph_labels = 12
+    gcamorph_meta = 13    # introduced for mgz warpfield
+    gcamorph_affine = 14  # introduced for mgz warpfield (m3z outputs the same information under xform)
     old_surf_geom = 20
     surf_geom = 21
+    surf_dataspace = 22   # surface tag - surface [x y z] space
+    surf_matrixdata = 23  # surface tag - transform matrix going from surf_dataspace to surf_transformedspace
+    surf_transformedspace = 24  # surface tag - surface [x y z] space after applying the transform matrix
     old_xform = 30
     xform = 31
     group_avg_area = 32
@@ -27,14 +32,22 @@ class tags:
     pedir = 41
     mri_frame = 42
     fieldstrength = 43
+    orig_ras2vox  = 44   # ???
 
 
 # these are FreeSurfer tags that have a
 # buffer with hardcoded length
+# notes:
+# 1. there seems to be some special handling of 'old_xform' for backwards compatibility
+# 2. 'xform' is used in both m3z and mgz, but with different formats.
+#    it is lengthless for m3z, but it has a data-length for mgz
 lengthless_tags = [
     tags.old_surf_geom,
     tags.old_real_ras,
     tags.old_colortable,
+    tags.gcamorph_geom,
+    tags.gcamorph_type,
+    tags.gcamorph_labels
 ]
 
 

surfa/io/utils.py

@@ -98,3 +98,68 @@ def write_bytes(file, value, dtype):
         Datatype to save as.
     """
     file.write(np.asarray(value).astype(dtype, copy=False).tobytes())
+
+
+# read VOL_GEOM
+# also see VOL_GEOM.read() in mri.h 
+def read_volgeom(file):
+    volgeom = dict(
+        valid  = read_bytes(file, '>i4', 1),
+
+        width  = read_bytes(file, '>i4', 1),
+        height = read_bytes(file, '>i4', 1),
+        depth  = read_bytes(file, '>i4', 1),
+
+        xsize  = read_bytes(file, '>f4', 1),
+        ysize  = read_bytes(file, '>f4', 1),
+        zsize  = read_bytes(file, '>f4', 1),
+
+        x_r    = read_bytes(file, '>f4', 1),
+        x_a    = read_bytes(file, '>f4', 1),
+        x_s    = read_bytes(file, '>f4', 1),
+        y_r    = read_bytes(file, '>f4', 1),
+        y_a    = read_bytes(file, '>f4', 1),
+        y_s    = read_bytes(file, '>f4', 1),
+        z_r    = read_bytes(file, '>f4', 1),
+        z_a    = read_bytes(file, '>f4', 1),
+        z_s    = read_bytes(file, '>f4', 1),
+
+        c_r    = read_bytes(file, '>f4', 1),
+        c_a    = read_bytes(file, '>f4', 1),
+        c_s    = read_bytes(file, '>f4', 1),
+
+        fname  = file.read(512).decode('utf-8').rstrip('\x00')
+        )
+    return volgeom
+
+
+# output VOL_GEOM
+# also see VOL_GEOM.write() in mri.h
+def write_volgeom(file, volgeom):
+    write_bytes(file, volgeom['valid'], '>i4')
+
+    write_bytes(file, volgeom['width'], '>i4')
+    write_bytes(file, volgeom['height'], '>i4')
+    write_bytes(file, volgeom['depth'], '>i4')
+
+    write_bytes(file, volgeom['xsize'], '>f4')
+    write_bytes(file, volgeom['ysize'], '>f4')
+    write_bytes(file, volgeom['zsize'], '>f4')
+
+    write_bytes(file, volgeom['x_r'], '>f4')
+    write_bytes(file, volgeom['x_a'], '>f4')
+    write_bytes(file, volgeom['x_s'], '>f4')
+    write_bytes(file, volgeom['y_r'], '>f4')
+    write_bytes(file, volgeom['y_a'], '>f4')
+    write_bytes(file, volgeom['y_s'], '>f4')
+    write_bytes(file, volgeom['z_r'], '>f4')
+    write_bytes(file, volgeom['z_a'], '>f4')
+    write_bytes(file, volgeom['z_s'], '>f4')
+
+    write_bytes(file, volgeom['c_r'], '>f4')
+    write_bytes(file, volgeom['c_a'], '>f4')
+    write_bytes(file, volgeom['c_s'], '>f4')
+
+    # output 512 bytes padded with '/x00'
+    file.write(volgeom['fname'].ljust(512, '\x00').encode('utf-8'))
+