Merge pull request #47 from yhuang43/master

add shear components to VOL_GEOM

surfa/io/framed.py

@@ -352,6 +352,16 @@ def load(self, filename, atype):
                     arr.metadata['target-valid'] = valid
                     arr.metadata['target-fname'] = fname
 
+                # gcamorph src & trg geoms (mgz warp)
+                elif tag == fsio.tags.gcamorph_geom_plusshear:
+                    arr.source, valid, fname = read_geom(file, shearless=False)
+                    arr.metadata['source-valid'] = valid
+                    arr.metadata['source-fname'] = fname
+
+                    arr.target, valid, fname = read_geom(file, shearless=False)
+                    arr.metadata['target-valid'] = valid
+                    arr.metadata['target-fname'] = fname
+
                 # gcamorph meta (mgz warp: int int float)
                 elif tag == fsio.tags.gcamorph_meta:
                     arr.format = read_bytes(file, dtype='>i4')
@@ -463,8 +473,10 @@ def save(self, arr, filename, intent=FramedArrayIntents.mri):
             write_bytes(file, arr.metadata.get('field-strength', 0.0), '>f4')
 
             # gcamorph geom and gcamorph meta for mgz warp
+            # output both fsio.tags.gcamorph_geom and fsio.tags.gcamorph_geom_plusshear
             if intent == FramedArrayIntents.warpmap:
                 # gcamorph src & trg geoms (mgz warp)
+                # fsio.tags.gcamorph_geom
                 fsio.write_tag(file, fsio.tags.gcamorph_geom)
                 write_geom(file,
                            geom=arr.source,
@@ -475,6 +487,21 @@ def save(self, arr, filename, intent=FramedArrayIntents.mri):
                            valid=arr.metadata.get('target-valid', True),
                            fname=arr.metadata.get('target-fname', ''))
 
+                # fsio.tags.gcamorph_geom_plusshear
+                # gcamorph_geom_plusshear has a length, datalength needs to be consistent with write_geom()
+                datalength = 1200
+                fsio.write_tag(file, fsio.tags.gcamorph_geom_plusshear, datalength)
+                write_geom(file,
+                           geom=arr.source,
+                           valid=arr.metadata.get('source-valid', True),
+                           fname=arr.metadata.get('source-fname', ''),
+                           shearless=False)
+                write_geom(file,
+                           geom=arr.target,
+                           valid=arr.metadata.get('target-valid', True),
+                           fname=arr.metadata.get('target-fname', ''),
+                           shearless=False)
+
                 # gcamorph meta (mgz warp: int int float)
                 fsio.write_tag(file, fsio.tags.gcamorph_meta, 12)
                 write_bytes(file, arr.format, dtype='>i4')

surfa/io/fsio.py

@@ -19,6 +19,7 @@ class tags:
     gcamorph_labels = 12
     gcamorph_meta = 13    # introduced for mgz warpfield
     gcamorph_affine = 14  # introduced for mgz warpfield (m3z outputs the same information under xform)
+    gcamorph_geom_plusshear = 15  # information output under gcamorph_geom + shear components
     old_surf_geom = 20
     surf_geom = 21
     surf_dataspace = 22   # surface tag - surface [x y z] space

surfa/io/fsnifti1extension.py

@@ -15,15 +15,17 @@ class FSNifti1Extension:
     This class handles Freesurfer Nifti1 header extension IO.
 
     Class variables:
+      _verbose:      if it is True, output debug information
       _content:      FSNifti1Extension.Content
     """
 
-    def __init__(self):
+    def __init__(self, verbose=False):
         """
         FSNifti1Extension Constructor
         """
 
         # initialization
+        self._verbose = verbose
         self._content = FSNifti1Extension.Content()
 
 
@@ -56,8 +58,9 @@ def read(self, fileobj, esize, offset=0):
         self.content.intent = int.from_bytes(fsexthdr[1:3], byteorder='big')
         self.content.version = fsexthdr[3]
 
-        print(f'[DEBUG] FSNifti1Extension.read(): esize = {esize:6d}')
-        print(f'[DEBUG] FSNifti1Extension.read(): endian = \'{self.content.endian:c}\', intent = {self.content.intent:d}, version = {self.content.version:d}')
+        if (self._verbose):
+            print(f'[DEBUG] FSNifti1Extension.read(): esize = {esize:6d}')
+            print(f'[DEBUG] FSNifti1Extension.read(): endian = \'{self.content.endian:c}\', intent = {self.content.intent:d}, version = {self.content.version:d}')
 
         # process Freesurfer Nifti1 extension tag data
         tagdatalen = esize - 12  # exclude esize, ecode, fsexthdr
@@ -66,7 +69,8 @@ def read(self, fileobj, esize, offset=0):
             # read tagid (4 bytes), data-length (8 bytes)
             (tag, length) = FSNifti1Extension.read_tag(fileobj)
 
-            print(f'[DEBUG] FSNifti1Extension.read(): remaining taglen = {tagdatalen:6d} (tag = {tag:2d}, length = {length:5d})')
+            if (self._verbose): 
+                print(f'[DEBUG] FSNifti1Extension.read(): remaining taglen = {tagdatalen:6d} (tag = {tag:2d}, length = {length:5d})')
 
             if (tag == 0):
                 break
@@ -88,20 +92,20 @@ def read(self, fileobj, esize, offset=0):
                 dof = iou.read_int(fileobj, length)
                 self.content.dof = dof
 
-            # ras_xform (TAG_RAS_XFORM = 8)
-            elif (tag == FSNifti1Extension.Tags.ras_xform):
-                self.content.ras_xform = dict(
-                    rotation = iou.read_bytes(fileobj, '>f4', 9).reshape((3, 3), order='F'),
-                    center   = iou.read_bytes(fileobj, '>f4', 3)
-                )
-
             # gcamorph src & trg geoms (warp) (TAG_GCAMORPH_GEOM = 10)
             elif (tag == FSNifti1Extension.Tags.gcamorph_geom):
                 if (not self.content.warpmeta):
                     self.content.warpmeta = {}
 
                 (self.content.warpmeta['source-geom'], self.content.warpmeta['source-valid'], self.content.warpmeta['source-fname']) = iou.read_geom(fileobj, niftiheaderext=True)
                 (self.content.warpmeta['target-geom'], self.content.warpmeta['target-valid'], self.content.warpmeta['target-fname']) = iou.read_geom(fileobj, niftiheaderext=True)
+            # gcamorph src & trg geoms (warp) (TAG_GCAMORPH_GEOM_PLUSSHEAR = 15)
+            elif (tag == FSNifti1Extension.Tags.gcamorph_geom_plusshear):
+                if (not self.content.warpmeta):
+                    self.content.warpmeta = {}
+
+                (self.content.warpmeta['source-geom'], self.content.warpmeta['source-valid'], self.content.warpmeta['source-fname']) = iou.read_geom(fileobj, niftiheaderext=True, shearless=False)
+                (self.content.warpmeta['target-geom'], self.content.warpmeta['target-valid'], self.content.warpmeta['target-fname']) = iou.read_geom(fileobj, niftiheaderext=True, shearless=False)
             # gcamorph meta (warp: int int float) (TAG_GCAMORPH_META = 13)
             elif (tag == FSNifti1Extension.Tags.gcamorph_meta):
                 if (not self.content.warpmeta):
@@ -129,7 +133,8 @@ def read(self, fileobj, esize, offset=0):
             # check if we reach the end
             tagdatalen -= (len_tagheader + length)
             if (tagdatalen < len_tagheader):
-                print(f'[DEBUG] FSNifti1Extension.read(): remaining taglen = {tagdatalen:6d}')
+                if (self._verbose):
+                    print(f'[DEBUG] FSNifti1Extension.read(): remaining taglen = {tagdatalen:6d}')
                 break;
 
         return self.content
@@ -163,7 +168,8 @@ def write(self, fileobj, content, countbytesonly=False):
             iou.write_int(fileobj, content.version, size=1, byteorder='big')
 
         num_bytes = 4
-        print(f'[DEBUG] FSNifti1Extension.write():              dlen = {num_bytes:6d}')
+        if (self._verbose):
+            print(f'[DEBUG] FSNifti1Extension.write():              dlen = {num_bytes:6d}')
 
         # tag data
         addtaglength = 12
@@ -175,7 +181,8 @@ def write(self, fileobj, content, countbytesonly=False):
             tag = FSNifti1Extension.Tags.gcamorph_geom
             length = FSNifti1Extension.getlen_gcamorph_geom(source_fname, target_fname)
             num_bytes += length + addtaglength
-            print(f'[DEBUG] FSNifti1Extension.write(): +{length:5d}, +{addtaglength:d}, dlen = {num_bytes:6d}, TAG = {tag:2d}')
+            if (self._verbose):
+                print(f'[DEBUG] FSNifti1Extension.write(): +{length:5d}, +{addtaglength:d}, dlen = {num_bytes:6d}, TAG = {tag:2d}')
             if (not countbytesonly):
                 FSNifti1Extension.write_tag(fileobj, tag, length)
                 iou.write_geom(fileobj,
@@ -189,11 +196,32 @@ def write(self, fileobj, content, countbytesonly=False):
                            fname=target_fname,
                            niftiheaderext=True)
 
+            tag = FSNifti1Extension.Tags.gcamorph_geom_plusshear
+            length = FSNifti1Extension.getlen_gcamorph_geom(source_fname, target_fname, shearless=False)
+            num_bytes += length + addtaglength
+            if (self._verbose):
+                print(f'[DEBUG] FSNifti1Extension.write(): +{length:5d}, +{addtaglength:d}, dlen = {num_bytes:6d}, TAG = {tag:2d}')
+            if (not countbytesonly):
+                FSNifti1Extension.write_tag(fileobj, tag, length)
+                iou.write_geom(fileobj,
+                           geom=content.warpmeta['source-geom'],
+                           valid=content.warpmeta.get('source-valid', True),
+                           fname=source_fname,
+                           niftiheaderext=True,
+                           shearless=False)
+                iou.write_geom(fileobj,
+                           geom=content.warpmeta['target-geom'],
+                           valid=content.warpmeta.get('target-valid', True),
+                           fname=target_fname,
+                           niftiheaderext=True,
+                           shearless=False)
+
             # gcamorph meta (warp: int int float) (TAG_GCAMORPH_META = 13)
             tag = FSNifti1Extension.Tags.gcamorph_meta
             length = 12
             num_bytes += length + addtaglength
-            print(f'[DEBUG] FSNifti1Extension.write(): +{length:5d}, +{addtaglength:d}, dlen = {num_bytes:6d}, TAG = {tag:2d}')
+            if (self._verbose):
+                print(f'[DEBUG] FSNifti1Extension.write(): +{length:5d}, +{addtaglength:d}, dlen = {num_bytes:6d}, TAG = {tag:2d}')
             if (not countbytesonly):
                 FSNifti1Extension.write_tag(fileobj, tag, length)
                 iou.write_bytes(fileobj, content.warpmeta['format'], dtype='>i4')
@@ -206,7 +234,8 @@ def write(self, fileobj, content, countbytesonly=False):
             tag = FSNifti1Extension.Tags.end_data
             length = 1  # extra char '*'
             num_bytes += length + addtaglength
-            print(f'[DEBUG] FSNifti1Extension.write(): +{length:5d}, +{addtaglength:d}, dlen = {num_bytes:6d}, TAG = {tag:2d}')
+            if (self._verbose):
+                print(f'[DEBUG] FSNifti1Extension.write(): +{length:5d}, +{addtaglength:d}, dlen = {num_bytes:6d}, TAG = {tag:2d}')
             if (not countbytesonly):
                 FSNifti1Extension.write_tag(fileobj, tag, length)
                 extrachar = '*'
@@ -220,7 +249,8 @@ def write(self, fileobj, content, countbytesonly=False):
             tag = FSNifti1Extension.Tags.old_colortable
             length = FSNifti1Extension.getlen_labels(content.labels)
             num_bytes += length + addtaglength
-            print(f'[DEBUG] FSNifti1Extension.write(): +{length:5d}, +{addtaglength:d}, dlen = {num_bytes:6d}, TAG = {tag:2d}')
+            if (self._verbose):
+                print(f'[DEBUG] FSNifti1Extension.write(): +{length:5d}, +{addtaglength:d}, dlen = {num_bytes:6d}, TAG = {tag:2d}')
             if (not countbytesonly):
                 FSNifti1Extension.write_tag(fileobj, tag, length)
                 fsio.write_binary_lookup_table(fileobj, content.labels)
@@ -231,7 +261,8 @@ def write(self, fileobj, content, countbytesonly=False):
                 tag = FSNifti1Extension.Tags.history
                 length = len(hist)
                 num_bytes += length + addtaglength
-                print(f'[DEBUG] FSNifti1Extension.write(): +{length:5d}, +{addtaglength:d}, dlen = {num_bytes:6d}, TAG = {tag:2d}')
+                if (self._verbose):
+                    print(f'[DEBUG] FSNifti1Extension.write(): +{length:5d}, +{addtaglength:d}, dlen = {num_bytes:6d}, TAG = {tag:2d}')
                 if (not countbytesonly):
                     FSNifti1Extension.write_tag(fileobj, tag, length)
                     fileobj.write(hist.encode('utf-8'))
@@ -240,28 +271,19 @@ def write(self, fileobj, content, countbytesonly=False):
         tag = FSNifti1Extension.Tags.dof
         length = 4
         num_bytes += length + addtaglength
-        print(f'[DEBUG] FSNifti1Extension.write(): +{length:5d}, +{addtaglength:d}, dlen = {num_bytes:6d}, TAG = {tag:2d}')
+        if (self._verbose):
+            print(f'[DEBUG] FSNifti1Extension.write(): +{length:5d}, +{addtaglength:d}, dlen = {num_bytes:6d}, TAG = {tag:2d}')
         if (not countbytesonly):
             FSNifti1Extension.write_tag(fileobj, tag, length)
             iou.write_int(fileobj, content.dof, size=4)
 
-        # ras_xform (TAG_RAS_XFORM = 8)
-        if (content.ras_xform):
-            tag = FSNifti1Extension.Tags.ras_xform
-            length = 48
-            num_bytes += length + addtaglength
-            print(f'[DEBUG] FSNifti1Extension.write(): +{length:5d}, +{addtaglength:d}, dlen = {num_bytes:6d}, TAG = {tag:2d}')
-            if (not countbytesonly):
-                FSNifti1Extension.write_tag(fileobj, tag, length)
-                iou.write_bytes(fileobj, np.ravel(content.ras_xform['rotation'], order='F'), '>f4')
-                iou.write_bytes(fileobj, content.ras_xform['center'], '>f4')
-
         # scan_parameters (TAG_SCAN_PARAMETERS = 45)
         if (content.scan_parameters):
             tag = FSNifti1Extension.Tags.scan_parameters
             length = 20 + len(content.scan_parameters['pedir'])
             num_bytes += length + addtaglength
-            print(f'[DEBUG] FSNifti1Extension.write(): +{length:5d}, +{addtaglength:d}, dlen = {num_bytes:6d}, TAG = {tag:2d}')
+            if (self._verbose):
+                print(f'[DEBUG] FSNifti1Extension.write(): +{length:5d}, +{addtaglength:d}, dlen = {num_bytes:6d}, TAG = {tag:2d}')
             if (not countbytesonly):
                 FSNifti1Extension.write_tag(fileobj, tag, length)
                 iou.write_bytes(fileobj, content.scan_parameters['te'], '>f4')
@@ -284,7 +306,8 @@ def write(self, fileobj, content, countbytesonly=False):
         tag = FSNifti1Extension.Tags.end_data
         length = 1  # extra char '*'
         num_bytes += length + addtaglength
-        print(f'[DEBUG] FSNifti1Extension.write(): +{length:5d}, +{addtaglength:d}, dlen = {num_bytes:6d}, TAG = {tag:2d}')
+        if (self._verbose):
+            print(f'[DEBUG] FSNifti1Extension.write(): +{length:5d}, +{addtaglength:d}, dlen = {num_bytes:6d}, TAG = {tag:2d}')
         if (not countbytesonly):
             FSNifti1Extension.write_tag(fileobj, tag, length)
             extrachar = '*'
@@ -364,7 +387,7 @@ def getlen_labels(labels):
 
 
     @staticmethod
-    def getlen_gcamorph_geom(fname_source, fname_target):
+    def getlen_gcamorph_geom(fname_source, fname_target, shearless=True):
         """
         Calculate total bytes that will be written for the labels.
 
@@ -383,6 +406,8 @@ def getlen_gcamorph_geom(fname_source, fname_target):
 
         # See freesurfer/utils/fstagsio.cpp::getlen_gcamorph_geom().
         num_bytes = 2 * 80
+        if (not shearless):
+            num_bytes += 2 * 12  # 2 * 3 * 4 bytes (3 shear float numbers each)
         num_bytes += len(fname_source)
         num_bytes += len(fname_target)
 
@@ -407,16 +432,20 @@ class Tags:
 
         This class defines the tags recognized in surfa.
         It is a subset of tag IDs defined in freesurfer/include/tags.h
+
+        TAG_RAS_XFORM is also output by Freesurfer for FS nifti1 header extension only to store shearless rotation and center parameters.
+        This is to prevent precision loss from sform/qform decompose and compose when nii is read and written in Freesurfer.
+        The tag is not necessary for Surfa.
         """
 
-        old_colortable  = 1    # TAG_OLD_COLORTABLE
-        history         = 3    # TAG_CMDLINE
-        dof             = 7    # TAG_DOF
-        ras_xform       = 8    # TAG_RAS_XFORM
-        gcamorph_geom   = 10   # TAG_GCAMORPH_GEOM
-        gcamorph_meta   = 13   # TAG_GCAMORPH_META
-        scan_parameters = 45   # TAG_SCAN_PARAMETERS
-        end_data        = -1   # TAG_END_NIIHDREXTENSION
+        old_colortable  = 1           # TAG_OLD_COLORTABLE
+        history         = 3           # TAG_CMDLINE
+        dof             = 7           # TAG_DOF
+        gcamorph_geom   = 10          # TAG_GCAMORPH_GEOM
+        gcamorph_meta   = 13          # TAG_GCAMORPH_META
+        gcamorph_geom_plusshear = 15  # information output under gcamorph_geom + shear components
+        scan_parameters = 45          # TAG_SCAN_PARAMETERS
+        end_data        = -1          # TAG_END_NIIHDREXTENSION
 
 
     class Content:

surfa/io/utils.py

@@ -108,7 +108,7 @@ def write_bytes(file, value, dtype):
     file.write(np.asarray(value).astype(dtype, copy=False).tobytes())
 
 
-def read_geom(file, niftiheaderext=False):
+def read_geom(file, niftiheaderext=False, shearless=True):
     """
     Read an image geometry from a binary file buffer. See VOL_GEOM.read() in mri.h.
 
@@ -129,12 +129,21 @@ def read_geom(file, niftiheaderext=False):
         If True, write for nifti header extension.
     """
     valid = bool(read_bytes(file, '>i4', 1))
-    geom = ImageGeometry(
-        shape=read_bytes(file, '>i4', 3).astype(int),
-        voxsize=read_bytes(file, '>f4', 3),
-        rotation=read_bytes(file, '>f4', 9).reshape((3, 3), order='F'),
-        center=read_bytes(file, '>f4', 3),
-    )
+    if (shearless):
+        geom = ImageGeometry(
+            shape=read_bytes(file, '>i4', 3).astype(int),
+            voxsize=read_bytes(file, '>f4', 3),
+            rotation=read_bytes(file, '>f4', 9).reshape((3, 3), order='F'),
+            center=read_bytes(file, '>f4', 3),
+        )
+    else:
+        geom = ImageGeometry(
+            shape=read_bytes(file, '>i4', 3).astype(int),
+            voxsize=read_bytes(file, '>f4', 3),
+            rotation=read_bytes(file, '>f4', 9).reshape((3, 3), order='F'),
+            center=read_bytes(file, '>f4', 3),
+            shear=read_bytes(file, '>f4', 3)
+        )        
 
     len_fname_max = 512
     if (not niftiheaderext):
@@ -154,7 +163,7 @@ def read_geom(file, niftiheaderext=False):
     return geom, valid, fname
 
 
-def write_geom(file, geom, valid=True, fname='', niftiheaderext=False):
+def write_geom(file, geom, valid=True, fname='', niftiheaderext=False, shearless=True):
     """
     Write an image geometry to a binary file buffer. See VOL_GEOM.write() in mri.h.
 
@@ -173,11 +182,15 @@ def write_geom(file, geom, valid=True, fname='', niftiheaderext=False):
     """
     write_bytes(file, valid, '>i4')
 
-    voxsize, rotation, center = geom.shearless_components()
+    voxsize, rotation, center, shear = geom.voxsize, geom.rotation, geom.center, geom.shear
+    if (shearless):
+        voxsize, rotation, center = geom.shearless_components()
     write_bytes(file, geom.shape, '>i4')
     write_bytes(file, voxsize, '>f4')
     write_bytes(file, np.ravel(rotation, order='F'), '>f4')
     write_bytes(file, center, '>f4')
+    if (not shearless):
+        write_bytes(file, shear, '>f4')
 
     len_fname_max = 512
     if (not niftiheaderext):