From b5c579cc970953f3c6ab5f958b664a117118ef49 Mon Sep 17 00:00:00 2001 From: Yujing Huang Date: Mon, 25 Mar 2024 16:23:32 -0400 Subject: [PATCH 1/3] handle Freesurfer Nifti1 header extension 1. implement class FSNifti1Extension (io.fsnifti1extension.py) to handle Freesurfer Nifti1 header extension IO 2. fix tr handling in NiftiArrayIO for converting between mgz and nii (freesurfer saves tr as msec internally) 3. set MGHArrayIO intent using metadata['intent'] value --- surfa/io/framed.py | 72 +++++ surfa/io/fsnifti1extension.py | 593 ++++++++++++++++++++++++++++++++++ surfa/io/utils.py | 2 +- 3 files changed, 666 insertions(+), 1 deletion(-) create mode 100644 surfa/io/fsnifti1extension.py diff --git a/surfa/io/framed.py b/surfa/io/framed.py index f801e56..7c99880 100644 --- a/surfa/io/framed.py +++ b/surfa/io/framed.py @@ -409,6 +409,7 @@ def save(self, arr, filename, intent=FramedArrayIntents.mri): shape[-1] = arr.nframes # begin writing header + intent = arr.metadata.get('intent', intent) version = ((intent & 0xff) << 8) | 1 # encode intent in version write_bytes(file, version, '>u4') # version write_bytes(file, shape, '>u4') # shape @@ -546,6 +547,54 @@ def load(self, filename, atype): arr.geom.voxsize = arr.geom.voxsize * 1000 elif spatial_units_code == self.units_to_code['um']: arr.geom.voxsize = arr.geom.voxsize * 0.001 + + # freesurfer saves tr as msec internally + time_units_factor = 0.0 + if (arr.metadata['frame_units'] == 'sec'): + time_units_factor = 1000.0 + elif (arr.metadata['frame_units'] == 'msec'): + time_units_factor = 1.0 + elif (arr.metadata['frame_units'] == 'usec'): + time_units_factor = 0.001 + + arr.metadata['tr'] = arr.metadata['frame_dim'] * time_units_factor + + # handle nifti1 header extension + niiextsions = nii.header.extensions + if (not niiextsions): + print("[DEBUG] NiftiArrayIO.load(): no header extensions found!") + return arr + + # try to find freesurfer nifti1 header extension + index = -1 + ecodes = niiextsions.get_codes() + for ind, ecode in enumerate(ecodes): + if (ecode == 14): + index = ind + break + + if (index < 0): + return arr + + # found freesurfer nifti1 header extension + import io + from surfa.io.fsnifti1extension import FSNifti1Extension + + # retrieve freesurfer header extension data in bytes + fsext = nii.header.extensions[index] + esize = fsext.get_sizeondisk() + fsext_bytedata = fsext.get_content() + + # create file-like object from bytes data + bytes_reader = io.BytesIO(fsext_bytedata) + + # read freesurfer header extension data into FSNifti1Extension.Content object + fsextreader = FSNifti1Extension() + fsextcontent = fsextreader.read(bytes_reader, esize) + + # update arr metadata + fsextcontent.update_framedimage(arr) + return arr def save(self, arr, filename): @@ -581,6 +630,10 @@ def save(self, arr, filename): # initialize spatial and temporal spacing nii.header['pixdim'][:] = 1 nii.header['pixdim'][4] = arr.metadata.get('frame_dim', 1) + + tr = arr.metadata.get('tr') + if (tr is not None): + nii.header['pixdim'][4] = tr / 1000.0 # for now we pretty much have to enforce spatial units of mm # and if frame units isn't specified, fallback to seconds @@ -604,6 +657,25 @@ def save(self, arr, filename): nii.set_qform(arr.geom.vox2world.matrix, arr.metadata.get('qform_code', 1)) nii.header['pixdim'][1:4] = arr.geom.voxsize.astype(np.float32) + # add freesurfer nifti1 header extension + import io + from surfa.io.fsnifti1extension import FSNifti1Extension + + # create FSNifti1Extension.Content object from arr + fsextcontent = FSNifti1Extension.Content(arr) + + # create file-like object + bytes_writer = io.BytesIO() + + # write freesurfer header extension data to buffer + fsext = FSNifti1Extension() + esize = fsext.write(bytes_writer, fsextcontent) + + # convert FSNifti1Extension.Content to bytes + fsext_bytesdata = bytes_writer.getvalue() + niiext = self.nib.nifti1.Nifti1Extension(14, fsext_bytesdata) + nii.header.extensions.append(niiext) + # write self.nib.save(nii, filename) diff --git a/surfa/io/fsnifti1extension.py b/surfa/io/fsnifti1extension.py new file mode 100644 index 0000000..8c1fd4f --- /dev/null +++ b/surfa/io/fsnifti1extension.py @@ -0,0 +1,593 @@ +import os +import gzip +import warnings +import numpy as np + +from surfa.io import fsio +from surfa.io import utils as iou +from surfa.transform.geometry import ImageGeometry +from surfa.core.labels import LabelLookup +from surfa.core.framed import FramedArrayIntents + + +class FSNifti1Extension: + """ + This class handles Freesurfer Nifti1 header extension IO. + + Class variables: + _content: FSNifti1Extension.Content + """ + + def __init__(self): + """ + FSNifti1Extension Constructor + """ + + # initialization + self._content = FSNifti1Extension.Content() + + + def read(self, fileobj, esize, offset=0): + """ + Read Freesurfer Nifti1 header extension data from the file-like object. + + Parameters + ---------- + fileobj : file-like object + opened file-like object + esize : int + nifti1 header extension size including sizes of esize and ecode + offset : int + offset to Freesurfer Nifti1 Header Extension + + Returns + ------- + content : FSNifti1Extension.Content + """ + + if (offset > 0): + fileobj.seek(offset) # seek to Freesurfer Header Extension + + # freesurfer nifti header extension data is in big endian + # the first 4 bytes are as following: + # endian ">" (1 byte), intent (unsigned short, 2 bytes), version (1 byte) + fsexthdr = fileobj.read(4) + self.content.endian = fsexthdr[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}') + + # process Freesurfer Nifti1 extension tag data + tagdatalen = esize - 12 # exclude esize, ecode, fsexthdr + len_tagheader = 12 # tagid (4 bytes), data-length (8 bytes) + while True: + # 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 (tag == 0): + break + + # embedded lookup table (TAG_OLD_COLORTABLE = 1) + elif (tag == FSNifti1Extension.Tags.old_colortable): + self.content.labels = fsio.read_binary_lookup_table(fileobj) + + # command history (TAG_CMDLINE = 3) + elif (tag == FSNifti1Extension.Tags.history): + history = fileobj.read(length).decode('utf-8').rstrip('\x00') + if (self.content.history): + self.content.history.append(history) + else: + self.content.history = [history] + + # dof (TAG_DOF = 7) + elif (tag == FSNifti1Extension.Tags.dof): + 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) + (self.content.warpmeta['target-geom'], self.content.warpmeta['target-valid'], self.content.warpmeta['target-fname']) = iou.read_geom(fileobj) + + # gcamorph meta (warp: int int float) (TAG_GCAMORPH_META = 13) + elif (tag == FSNifti1Extension.Tags.gcamorph_meta): + if (not self.content.warpmeta): + self.content.warpmeta = {} + + self.content.warpmeta['format'] = iou.read_bytes(fileobj, dtype='>i4') + self.content.warpmeta['spacing'] = iou.read_bytes(fileobj, dtype='>i4') + self.content.warpmeta['exp_k'] = iou.read_bytes(fileobj, dtype='>f4') + + # scan_parameters (TAG_SCAN_PARAMETERS = 45) + elif (tag == FSNifti1Extension.Tags.scan_parameters): + self.content.scan_parameters = {} + self.content.scan_parameters['te'] = iou.read_bytes(fileobj, dtype='>f4') + self.content.scan_parameters['ti'] = iou.read_bytes(fileobj, dtype='>f4') + self.content.scan_parameters['flip_angle'] = iou.read_bytes(fileobj, dtype='>f8') # ??? double 8 bytes + self.content.scan_parameters['field_strength'] = iou.read_bytes(fileobj, dtype='>f4') + + len_pedir = length - 20 + self.content.scan_parameters['pedir'] = fileobj.read(len_pedir).decode('utf-8').rstrip('\x00') + + # skip everything else + else: + fileobj.read(length) + + # check if we reach the end + tagdatalen -= (len_tagheader + length) + if (tagdatalen < len_tagheader): + print(f'[DEBUG] FSNifti1Extension.read(): remaining taglen = {tagdatalen:6d}') + break; + + return self.content + + + def write(self, fileobj, content, countbytesonly=False): + """ + Write Freesurfer Nifti1 header extension data saved to the file-like object. + + Parameters + ---------- + fileobj : file-like object + opened file-like object + content : FSNifti1Extension.Content + Data to output or count bytes + countbytesonly : boolean + if True, count total bytes to output w/o writing to fileobj + + Returns + ------- + num_bytes : long + total bytes output or counted + """ + + # freesurfer nifti header extension data is in big endian + # the first 4 bytes are as following: + # endian ">" (1 byte), intent (unsigned short, 2 bytes), version (1 byte) + if (not countbytesonly): + fileobj.write(content.endian.encode('utf-8')) + iou.write_int(fileobj, content.intent, size=2, byteorder='big') + iou.write_int(fileobj, content.version, size=1, byteorder='big') + + num_bytes = 4 + print(f'[DEBUG] FSNifti1Extension.write(): dlen = {num_bytes:6d}') + + # tag data + addtaglength = 12 + if (content.intent == FramedArrayIntents.warpmap): + # gcamorph src & trg geoms (warp) (TAG_GCAMORPH_GEOM = 10) + tag = FSNifti1Extension.Tags.gcamorph_geom + length = 1176 # 2 * (4 * sizeof(int) + 15 * sizeof(float) + 512), see freesurfer/utils/fstagsio.cpp + 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_geom(fileobj, + geom=content.warpmeta['source-geom'], + valid=content.warpmeta.get('source-valid', True), + fname=content.warpmeta.get('source-fname', '')) + iou.write_geom(fileobj, + geom=content.warpmeta['target-geom'], + valid=content.warpmeta.get('target-valid', True), + fname=content.warpmeta.get('target-fname', '')) + + # 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 (not countbytesonly): + FSNifti1Extension.write_tag(fileobj, tag, length) + iou.write_bytes(fileobj, content.warpmeta['format'], dtype='>i4') + iou.write_bytes(fileobj, content.warpmeta.get('spacing', 1), dtype='>i4') + iou.write_bytes(fileobj, content.warpmeta.get('exp_k', 0.0), dtype='>f4') + + return num_bytes + + + # lookup table (TAG_OLD_COLORTABLE = 1) + if (content.labels): + 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 (not countbytesonly): + FSNifti1Extension.write_tag(fileobj, tag, length) + fsio.write_binary_lookup_table(fileobj, content.labels) + + # history (TAG_CMDLINE = 3) + if (content.history): + for hist in content.history: + 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 (not countbytesonly): + FSNifti1Extension.write_tag(fileobj, tag, length) + fileobj.write(hist.encode('utf-8')) + + # dof (TAG_DOF = 7) + 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 (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']) + 1 + 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, content.scan_parameters['te'], '>f4') + iou.write_bytes(fileobj, content.scan_parameters['ti'], '>f4') + iou.write_bytes(fileobj, content.scan_parameters['flip_angle'], '>f8') # ??? double 8 bytes + iou.write_bytes(fileobj, content.scan_parameters['field_strength'], '>f4') + fileobj.write(content.scan_parameters['pedir'].encode('utf-8')) + + return num_bytes + + + @staticmethod + def read_tag(fileobj): + """ + Reads the next FreeSurfer tag ID and associated data length. + + Parameters + ---------- + fileobj : file-like object + opened file-like object + + Returns + ------- + tag : FSNifti1Extension.Tags + Tag ID + length : long + data length of tagged data + """ + + tag = iou.read_int(fileobj) + length = iou.read_int(fileobj, size=8) + + return (tag, length) + + + @staticmethod + def write_tag(fileobj, tag, length): + """ + Write tag ID and associated data-length + + Parameters + ---------- + fileobj : file-like object + opened file-like object + tag : FSNifti1Extension.Tags + Tag ID + length : long + data length + """ + + iou.write_int(fileobj, tag) + iou.write_int(fileobj, length, size=8) + + + @staticmethod + def getlen_labels(labels): + """ + Calculate total bytes that will be written for the labels + + Parameters + ---------- + labels : LabelLookup + Labels used for calculation + + Returns + ------- + num_bytes : int + total bytes that will be written for the input labels + """ + + num_bytes = 12 # version, nentries, len(fname) + # io.fsio.write_binary_lookup_table() doesn't write fname, len(fname)=0 + num_bytes += 4 # num_entries + for index, element in labels.items(): + num_bytes += 8 # structure id, len(structure-name)+1 + num_bytes += len(element.name) + 1 # structure name + num_bytes += 16 # ri, gi, bi, t-ai + + return num_bytes + + + @property + def content(self): + return self._content + @content.setter + def content(self, content): + self._content = content + + + class Tags: + """ + Freesurfer Nifti1 header extension data is in big endian and has the following format: + endian ">" (1 byte), intent (2 bytes), version (1 byte) + tagid-1 data-length-1 tag-data-1 + tagid-2 data-length-2 tag-data-2 + ... + + This class defines the tags recognized in surfa. + 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 + + + class Content: + """ + Data class represents Freesurfer Nifti1 header extension data. + + Class variables: + _endian: big endian = '>', or little endian = '<' + _intent: FramedArrayIntents + _version: int + + _dof: int + _scan_parameters: dict { + 'te' : float, + 'ti' : float, + 'flip_angle' : double, + 'field_strength' : float, + 'pedir' : string + } + _ras_xform: dict { + 'rotation' : array(3, 3, float) + 'center' : array(1, 3, float) + } + _warpmeta: dict { + 'format' : Warp.Format + 'spacing' : int + 'exp_k' : float + 'source-valid' : int + 'source-fname' : string + 'source-geom' : ImageGeometry + 'target-valid' : int + 'target-fname' : string + 'target-geom' : ImageGeometry + } + _history: list[string] + _labels: LabelLookup + """ + + def __init__(self, framedimage=None): + """ + FSNifti1Extension.Content Constructor + """ + + self._endian = '>' + self._intent = FramedArrayIntents.mri + self._version = 1 + + self._dof = 1 + self._scan_parameters = None + self._ras_xform = None + self._warpmeta = None + self._history = None + self._labels = None + + if (framedimage): + self._from_framedimage(framedimage) + + + def update_framedimage(self, framedimage): + """ + Update input FramedImage metadata + + Parameters + ---------- + framedimage : FramedImage + input FramedImage to update + """ + + # update input framedimage metadata + framedimage.metadata['intent'] = self.intent + + if (self.intent == FramedArrayIntents.warpmap): + # gcamorph src & trg geoms (mgz warp) + framedimage.source = self.warpmeta['source-geom'] + framedimage.metadata['source-valid'] = self.warpmeta['source-valid'] + framedimage.metadata['source-fname'] = self.warpmeta['source-fname'] + + framedimage.target = self.warpmeta['target-geom'] + framedimage.metadata['target-valid'] = self.warpmeta['target-valid'] + framedimage.metadata['target-fname'] = self.warpmeta['target-fname'] + + # gcamorph meta (mgz warp: int int float) + framedimage.format = self.warpmeta['format'] + framedimage.metadata['spacing'] = self.warpmeta['spacing'] + framedimage.metadata['exp_k'] = self.warpmeta['exp_k'] + + return + + # update image-specific information + if (self.ras_xform): + framedimage.geom.update( + rotation=self.ras_xform['rotation'], + center=self.ras_xform['center'] + ) + + if (self.scan_parameters): + framedimage.metadata['phase-encode-direction'] = self.scan_parameters['pedir'] + framedimage.metadata['field-strength'] = self.scan_parameters['field_strength'] + + scan_params = {} + scan_params['fa'] = self.scan_parameters['flip_angle'] + scan_params['te'] = self.scan_parameters['te'] + scan_params['ti'] = self.scan_parameters['ti'] + framedimage.metadata.update(scan_params) + + if (self.history): + framedimage.metadata['history'] = self.history + + if (self.labels): + framedimage.labels = self.labels + + + def _from_framedimage(self, framedimage): + """ + Create FSNifti1Extension.Content object from input FramedImage metadata + + Parameters + ---------- + framedimage : FramedImage + input FramedImage + """ + + self.endian = '>' + self.version = 1 + self.intent = framedimage.metadata.get('intent', FramedArrayIntents.mri) + if (isinstance(self.intent, np.int_)): + self.intent = self.intent.item() # convert numpy int to python int + + self.dof = 1 + + if (self.intent == FramedArrayIntents.warpmap): + if (not self.warpmeta): + self.warpmeta = {} + + # gcamorph src & trg geoms (mgz warp) + self.warpmeta['source-geom'] = framedimage.source + self.warpmeta['source-valid'] = framedimage.metadata['source-valid'] + self.warpmeta['source-fname'] = framedimage.metadata['source-fname'] + + self.warpmeta['target-geom'] = framedimage.target + self.warpmeta['target-valid'] = framedimage.metadata['target-valid'] + self.warpmeta['target-fname'] = framedimage.metadata['target-fname'] + + # gcamorph meta (mgz warp: int int float) + self.warpmeta['format'] = framedimage.format + self.warpmeta['spacing'] = framedimage.metadata['spacing'] + self.warpmeta['exp_k'] = framedimage.metadata['exp_k'] + + return + + # update ras_xform + self.ras_xform = dict( + rotation = framedimage.geom.rotation, + center = framedimage.geom.center + ) + + # update scan_parameters + self.scan_parameters = dict( + pedir = framedimage.metadata.get('phase-encode-direction', 'UNKNOWN'), + field_strength = framedimage.metadata.get('field-strength'), + flip_angle = framedimage.metadata.get('fa', 0), + te = framedimage.metadata.get('te', 0), + ti = framedimage.metadata.get('ti', 0) + ) + + if (framedimage.metadata.get('history')): + self.history = framedimage.metadata['history'] + + if (framedimage.labels): + self.labels = framedimage.labels + + + @property + def endian(self): + return self._endian + @endian.setter + def endian(self, endian): + self._endian = endian + + @property + def intent(self): + return self._intent + @intent.setter + def intent(self, intent): + self._intent = intent + + @property + def version(self): + return self._version + @version.setter + def version(self, version): + self._version = version + + @property + def dof(self): + return self._dof + @dof.setter + def dof(self, dof): + self._dof = dof + + @property + def scan_parameters(self): + return self._scan_parameters + @scan_parameters.setter + def scan_parameters(self, scan_parameters): + self._scan_parameters = scan_parameters + + @property + def ras_xform(self): + return self._ras_xform + @ras_xform.setter + def ras_xform(self, ras_xform): + self._ras_xform = ras_xform + + @property + def warpmeta(self): + return self._warpmeta + @warpmeta.setter + def warpmeta(self, warpmeta): + self._warpmeta = warpmeta + + @property + def history(self): + return self._history + @history.setter + def history(self, history): + self._history = history + + @property + def labels(self): + return self._labels + @labels.setter + def labels(self, labels): + self._labels = labels + + + diff --git a/surfa/io/utils.py b/surfa/io/utils.py index 7352030..2c51507 100644 --- a/surfa/io/utils.py +++ b/surfa/io/utils.py @@ -143,7 +143,7 @@ def write_geom(file, geom, valid=True, fname=''): Parameters ---------- - file : BufferedReader + file : BufferedWriter Opened file buffer. geom : ImageGeometry Image geometry. From 9553feb094c964a28445654952e9d9f6d8725673 Mon Sep 17 00:00:00 2001 From: Yujing Huang Date: Tue, 26 Mar 2024 14:00:11 -0400 Subject: [PATCH 2/3] handle Freesurfer Nifti1 header extension different fname format for TAG_GCAMORPH_GEOM when output as nifti1 header extension: instead of padded 512 bytes, it has a variable length: fname-length fname if fname-length = 0, no fname output follows (io.utils.read_geom()/write_geom()) notes: 1. The ctab binary from surfa is slightly different from freesurfer output. Surfa doesn't write ctab filename, it sets the filename len to 0. See surfa.io.fsio.write_binary_lookup_table() and freesurfer/utils/colortab.cpp::znzCTABwriteIntoBinaryV2() 2. There is a bug in nibabel.nifti1.Nifti1Extensions.from_fileobj() that will truncate extension data if there are null characters or zeros in the data. The bug is in nibabel/nifti1.py line 629: 'evalue = evalue.rstrip(b'\x00')'. --- surfa/io/fsio.py | 4 ++++ surfa/io/fsnifti1extension.py | 42 ++++++++++++++++++++++++++++++----- surfa/io/utils.py | 39 +++++++++++++++++++++++++++----- 3 files changed, 74 insertions(+), 11 deletions(-) diff --git a/surfa/io/fsio.py b/surfa/io/fsio.py index fe598c2..ae1a2e1 100644 --- a/surfa/io/fsio.py +++ b/surfa/io/fsio.py @@ -123,6 +123,8 @@ def read_binary_lookup_table(file): version = iou.read_bytes(file, '>i4') max_id = iou.read_bytes(file, '>i4') file_name_size = iou.read_bytes(file, '>i4') + # if the file comes from surfa.io.fsio.write_binary_lookup_table(), file_name_size = 0. + # if the file comes from freesurfer/utils/colortab.cpp::znzCTABwriteIntoBinaryV2(), file_name_size > 0. file.read(file_name_size).decode('utf-8') total = iou.read_bytes(file, '>i4') @@ -154,6 +156,8 @@ def write_binary_lookup_table(file, labels): iou.write_bytes(file, -2, '>i4') iou.write_bytes(file, max(labels) + 1, '>i4') iou.write_bytes(file, 0, '>i4') + # filename length is set to 0, here filename associated with Labels is not output. + # this is different from freesurfer/utils/colortab.cpp::znzCTABwriteIntoBinaryV2(). iou.write_bytes(file, len(labels), '>i4') for index, element in labels.items(): iou.write_bytes(file, index, '>i4') diff --git a/surfa/io/fsnifti1extension.py b/surfa/io/fsnifti1extension.py index 8c1fd4f..1bfb9d5 100644 --- a/surfa/io/fsnifti1extension.py +++ b/surfa/io/fsnifti1extension.py @@ -100,9 +100,8 @@ def read(self, fileobj, esize, offset=0): 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) - (self.content.warpmeta['target-geom'], self.content.warpmeta['target-valid'], self.content.warpmeta['target-fname']) = iou.read_geom(fileobj) - + (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 meta (warp: int int float) (TAG_GCAMORPH_META = 13) elif (tag == FSNifti1Extension.Tags.gcamorph_meta): if (not self.content.warpmeta): @@ -170,8 +169,11 @@ def write(self, fileobj, content, countbytesonly=False): addtaglength = 12 if (content.intent == FramedArrayIntents.warpmap): # gcamorph src & trg geoms (warp) (TAG_GCAMORPH_GEOM = 10) + source_fname = content.warpmeta.get('source-fname', '') + target_fname = content.warpmeta.get('target-fname', '') + tag = FSNifti1Extension.Tags.gcamorph_geom - length = 1176 # 2 * (4 * sizeof(int) + 15 * sizeof(float) + 512), see freesurfer/utils/fstagsio.cpp + 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 (not countbytesonly): @@ -179,11 +181,13 @@ def write(self, fileobj, content, countbytesonly=False): iou.write_geom(fileobj, geom=content.warpmeta['source-geom'], valid=content.warpmeta.get('source-valid', True), - fname=content.warpmeta.get('source-fname', '')) + fname=source_fname, + niftiheaderext=True) iou.write_geom(fileobj, geom=content.warpmeta['target-geom'], valid=content.warpmeta.get('target-valid', True), - fname=content.warpmeta.get('target-fname', '')) + fname=target_fname, + niftiheaderext=True) # gcamorph meta (warp: int int float) (TAG_GCAMORPH_META = 13) tag = FSNifti1Extension.Tags.gcamorph_meta @@ -327,6 +331,32 @@ def getlen_labels(labels): return num_bytes + @staticmethod + def getlen_gcamorph_geom(fname_source, fname_target): + """ + Calculate total bytes that will be written for the labels. + + Parameters + ---------- + fname_source : string + File name associated with the source geometry + fname_target : string + File name associated with the target geometry + + Returns + ------- + num_bytes : int + total bytes that will be written for the source and target geometry + """ + + # See freesurfer/utils/fstagsio.cpp::getlen_gcamorph_geom(). + num_bytes = 2 * 80 + num_bytes += len(fname_source) + num_bytes += len(fname_target) + + return num_bytes + + @property def content(self): return self._content diff --git a/surfa/io/utils.py b/surfa/io/utils.py index 2c51507..7d57e36 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): +def read_geom(file, niftiheaderext=False): """ Read an image geometry from a binary file buffer. See VOL_GEOM.read() in mri.h. @@ -125,6 +125,8 @@ def read_geom(file): True if the geometry is valid. str File name associated with the geometry. + niftiheaderext : bool + If True, write for nifti header extension. """ valid = bool(read_bytes(file, '>i4', 1)) geom = ImageGeometry( @@ -133,11 +135,26 @@ def read_geom(file): rotation=read_bytes(file, '>f4', 9).reshape((3, 3), order='F'), center=read_bytes(file, '>f4', 3), ) - fname = file.read(512).decode('utf-8').rstrip('\x00') + + len_fname_max = 512 + if (not niftiheaderext): + fname = file.read(len_fname_max).decode('utf-8').rstrip('\x00') + else: + # variable length fname, if length = 0, no fname output follows + # read the first len_max bytes, skip the rest + len_fname = read_bytes(file, '>i4', 1) + to_read = len_fname + if (len_fname > len_fname_max): + to_read = len_fname_max + + fname = file.read(to_read).decode('utf-8').rstrip('\x00') + remaining = len_fname - to_read + if (remaining > 0): + file.read(remaining) return geom, valid, fname -def write_geom(file, geom, valid=True, fname=''): +def write_geom(file, geom, valid=True, fname='', niftiheaderext=False): """ Write an image geometry to a binary file buffer. See VOL_GEOM.write() in mri.h. @@ -151,6 +168,8 @@ def write_geom(file, geom, valid=True, fname=''): True if the geometry is valid. fname : str File name associated with the geometry. + niftiheaderext : bool + If True, write for nifti header extension. """ write_bytes(file, valid, '>i4') @@ -160,5 +179,15 @@ def write_geom(file, geom, valid=True, fname=''): write_bytes(file, np.ravel(rotation, order='F'), '>f4') write_bytes(file, center, '>f4') - # right-pad with '/x00' to 512 bytes - file.write(fname[:512].ljust(512, '\x00').encode('utf-8')) + len_fname_max = 512 + if (not niftiheaderext): + # right-pad with '/x00' to 512 bytes + file.write(fname[:len_fname_max].ljust(len_fname_max, '\x00').encode('utf-8')) + else: + # variable length fname, if length = 0, no fname output follows + len_fname = len(fname) + if (len_fname > len_fname_max): + len_fname = len_fname_max + write_bytes(file, len_fname, '>i4') + if (len_fname > 0): + file.write(fname.encode('utf-8')) From cc9cdf2ba20290b7ecb9a7a7dc3259c6f9703dbb Mon Sep 17 00:00:00 2001 From: Yujing Huang Date: Mon, 1 Apr 2024 12:52:54 -0400 Subject: [PATCH 3/3] handle Freesurfer Nifti1 header extension write TAG_END_NIIHDREXTENSION at the end of extension data to avoid the data to be truncated: TAG_END_NIIHDREXTENSION (-1) data-length (1) '*' note: If the extension data has trailing null characters or zeros at the end, nibabel.nifti1.Nifti1Extension.get_content() will truncate the data. See https://github.com/nipy/nibabel/blob/master/nibabel/nifti1.py#L629C1-L630C1, line 629: 'evalue = evalue.rstrip(b'\x00')' --- surfa/io/fsnifti1extension.py | 33 +++++++++++++++++++++++++++++++++ 1 file changed, 33 insertions(+) diff --git a/surfa/io/fsnifti1extension.py b/surfa/io/fsnifti1extension.py index 1bfb9d5..e217feb 100644 --- a/surfa/io/fsnifti1extension.py +++ b/surfa/io/fsnifti1extension.py @@ -199,6 +199,18 @@ def write(self, fileobj, content, countbytesonly=False): iou.write_bytes(fileobj, content.warpmeta['format'], dtype='>i4') iou.write_bytes(fileobj, content.warpmeta.get('spacing', 1), dtype='>i4') iou.write_bytes(fileobj, content.warpmeta.get('exp_k', 0.0), dtype='>f4') + + # write TAG_END_NIIHDREXTENSION at the end of extension data to avoid the data to be truncated: + # TAG_END_NIIHDREXTENSION (-1) data-length (1) '*' + # this needs to be the last tag. + 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 (not countbytesonly): + FSNifti1Extension.write_tag(fileobj, tag, length) + extrachar = '*' + fileobj.write(extrachar.encode('utf-8')) return num_bytes @@ -258,6 +270,26 @@ def write(self, fileobj, content, countbytesonly=False): iou.write_bytes(fileobj, content.scan_parameters['field_strength'], '>f4') fileobj.write(content.scan_parameters['pedir'].encode('utf-8')) + # end_data (TAG_END_NIIHDREXTENSION = -1) + """ + write TAG_END_NIIHDREXTENSION at the end of extension data to avoid the data to be truncated: + TAG_END_NIIHDREXTENSION (-1) data-length (1) '*' + this needs to be the last tag. + + If the extension data has trailing null characters or zeros at the end, + nibabel.nifti1.Nifti1Extension.get_content() will truncate the data. + See https://github.com/nipy/nibabel/blob/master/nibabel/nifti1.py#L629C1-L630C1, + line 629: 'evalue = evalue.rstrip(b'\x00')' + """ + 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 (not countbytesonly): + FSNifti1Extension.write_tag(fileobj, tag, length) + extrachar = '*' + fileobj.write(extrachar.encode('utf-8')) + return num_bytes @@ -384,6 +416,7 @@ class Tags: gcamorph_geom = 10 # TAG_GCAMORPH_GEOM gcamorph_meta = 13 # TAG_GCAMORPH_META scan_parameters = 45 # TAG_SCAN_PARAMETERS + end_data = -1 # TAG_END_NIIHDREXTENSION class Content: