From 0b1d1bb2ce0d5e6e00b981d1e007094e8b7986be Mon Sep 17 00:00:00 2001 From: Yujing Huang Date: Fri, 6 Dec 2024 13:04:06 -0500 Subject: [PATCH 1/6] add shear components to VOL_GEOM 1. introduce new tag TAG_GCAMORPH_GEOM_PLUSSHEAR to output TAG_GCAMORPH_GEOM+shear components 2. add 'shearless' option to VOL_GEOM IO functions to handle both TAG_GCAMORPH_GEOM and TAG_GCAMORPH_GEOM_PLUSSHEAR * notes: this commit is just the preparation for adding shear components to VOL_GEOM attached in warp (mgz, nii.gz) * todo: 1. for backward compatibility, both TAG_GCAMORPH_GEOM_PLUSSHEAR and TAG_GCAMORPH_GEOM are read and written; TAG_GCAMORPH_GEOM_PLUSSHEAR takes preference if both tags exist 2. take shears into consideration when computing vox2world matrix using VOL_GEOM --- surfa/io/fsio.py | 1 + surfa/io/utils.py | 31 ++++++++++++++++++++++--------- 2 files changed, 23 insertions(+), 9 deletions(-) diff --git a/surfa/io/fsio.py b/surfa/io/fsio.py index 6552bf3..dc74767 100644 --- a/surfa/io/fsio.py +++ b/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 diff --git a/surfa/io/utils.py b/surfa/io/utils.py index 7d57e36..beda899 100644 --- a/surfa/io/utils.py +++ b/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): From a3b3f45b8b7d13cebd19c16b1261a8aae4164fed Mon Sep 17 00:00:00 2001 From: Yujing Huang Date: Fri, 6 Dec 2024 13:27:16 -0500 Subject: [PATCH 2/6] add shear components to VOL_GEOM 1. add 'shearless' option to FSNifti1Extension.getlen_gcamorph_geom() 2. add 'verbose' attribute to class FSNifti1Extension to turn debug message on/off --- surfa/io/fsnifti1extension.py | 49 +++++++++++++++++++++++------------ 1 file changed, 33 insertions(+), 16 deletions(-) diff --git a/surfa/io/fsnifti1extension.py b/surfa/io/fsnifti1extension.py index 9f0e0df..d731889 100644 --- a/surfa/io/fsnifti1extension.py +++ b/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 @@ -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, @@ -193,7 +200,8 @@ def write(self, fileobj, content, countbytesonly=False): 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 +214,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 +229,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 +241,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,7 +251,8 @@ 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) @@ -250,7 +262,8 @@ def write(self, fileobj, content, countbytesonly=False): 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 (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, np.ravel(content.ras_xform['rotation'], order='F'), '>f4') @@ -261,7 +274,8 @@ def write(self, fileobj, content, countbytesonly=False): 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 +298,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 +379,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 +398,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) From 19503d772d97e09713165286852d431d0efd7781 Mon Sep 17 00:00:00 2001 From: Yujing Huang Date: Fri, 6 Dec 2024 14:21:00 -0500 Subject: [PATCH 3/6] add shear components to VOL_GEOM for backward compatibility, both TAG_GCAMORPH_GEOM_PLUSSHEAR and TAG_GCAMORPH_GEOM are read and written; TAG_GCAMORPH_GEOM_PLUSSHEAR takes preference if both tags exist --- surfa/io/framed.py | 25 ++++++++++++++++++++ surfa/io/fsnifti1extension.py | 44 ++++++++++++++++++++++++++++------- 2 files changed, 61 insertions(+), 8 deletions(-) diff --git a/surfa/io/framed.py b/surfa/io/framed.py index 4713032..799ff7a 100644 --- a/surfa/io/framed.py +++ b/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,19 @@ 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 + fsio.write_tag(file, fsio.tags.gcamorph_geom_plusshear) + 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') diff --git a/surfa/io/fsnifti1extension.py b/surfa/io/fsnifti1extension.py index d731889..9cea23c 100644 --- a/surfa/io/fsnifti1extension.py +++ b/surfa/io/fsnifti1extension.py @@ -106,6 +106,13 @@ def read(self, fileobj, esize, offset=0): (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): @@ -196,6 +203,26 @@ 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 @@ -426,14 +453,15 @@ class Tags: It is a subset of tag IDs defined in freesurfer/include/tags.h """ - 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 + ras_xform = 8 # TAG_RAS_XFORM + 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: From 1bf9b0ecadb5ac57458ef10fc3ebddc977df194b Mon Sep 17 00:00:00 2001 From: Yujing Huang Date: Fri, 6 Dec 2024 15:53:33 -0500 Subject: [PATCH 4/6] add shear components to VOL_GEOM (commit ) make TAG_GCAMORPH_GEOM_PLUSSHEAR length-less for mgz warp to be consistent with TAG_GCAMORPH_GEOM --- surfa/io/fsio.py | 1 + 1 file changed, 1 insertion(+) diff --git a/surfa/io/fsio.py b/surfa/io/fsio.py index dc74767..99bc687 100644 --- a/surfa/io/fsio.py +++ b/surfa/io/fsio.py @@ -47,6 +47,7 @@ class tags: tags.old_real_ras, tags.old_colortable, tags.gcamorph_geom, + tags.gcamorph_geom_plusshear, tags.gcamorph_type, tags.gcamorph_labels ] From 5fea7298bd851894c02bebf71b70fd42eaa7359c Mon Sep 17 00:00:00 2001 From: Yujing Huang Date: Thu, 12 Dec 2024 12:10:36 -0500 Subject: [PATCH 5/6] add shear components to VOL_GEOM revert previous commit for the data to be backward compatible, TAG_GCAMORPH_GEOM_PLUSSHEAR needs to have a length --- surfa/io/framed.py | 4 +++- surfa/io/fsio.py | 1 - 2 files changed, 3 insertions(+), 2 deletions(-) diff --git a/surfa/io/framed.py b/surfa/io/framed.py index 799ff7a..1a20a30 100644 --- a/surfa/io/framed.py +++ b/surfa/io/framed.py @@ -488,7 +488,9 @@ def save(self, arr, filename, intent=FramedArrayIntents.mri): fname=arr.metadata.get('target-fname', '')) # fsio.tags.gcamorph_geom_plusshear - fsio.write_tag(file, 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), diff --git a/surfa/io/fsio.py b/surfa/io/fsio.py index 99bc687..dc74767 100644 --- a/surfa/io/fsio.py +++ b/surfa/io/fsio.py @@ -47,7 +47,6 @@ class tags: tags.old_real_ras, tags.old_colortable, tags.gcamorph_geom, - tags.gcamorph_geom_plusshear, tags.gcamorph_type, tags.gcamorph_labels ] From eba77bb2413120f980d2b534f55bd752aaf059f6 Mon Sep 17 00:00:00 2001 From: Yujing Huang Date: Fri, 13 Dec 2024 14:07:27 -0500 Subject: [PATCH 6/6] remove TAG_RAS_XFORM output The tag is not necessary for Surfa. * TAG_RAS_XFORM is output by Freesurfer for FS nifti1 header extension only to store shearless rotation and center parameters. It is to prevent precision loss from sform/qform decompose and compose when nii is read and written in Freesurfer. --- surfa/io/fsnifti1extension.py | 24 ++++-------------------- 1 file changed, 4 insertions(+), 20 deletions(-) diff --git a/surfa/io/fsnifti1extension.py b/surfa/io/fsnifti1extension.py index 9cea23c..64df459 100644 --- a/surfa/io/fsnifti1extension.py +++ b/surfa/io/fsnifti1extension.py @@ -92,13 +92,6 @@ 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): @@ -284,18 +277,6 @@ def write(self, fileobj, content, countbytesonly=False): 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 - 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, 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 @@ -451,12 +432,15 @@ 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 gcamorph_geom_plusshear = 15 # information output under gcamorph_geom + shear components