From 2c6c91cc15265916e4d5d68990b43ce6e7f815d7 Mon Sep 17 00:00:00 2001 From: Yujing Huang Date: Mon, 18 Dec 2023 09:10:32 -0500 Subject: [PATCH] support mgz warp Freesurfer now supports mgz warp file, which follows mgz format with these tags and data: TAG_GCAMORPH_GEOM : gcamorph image (source) geom & gcamorph atlas (target) geom TAG_GCAMORPH_META : data format (int), spacing (int), exp_k (double) TAG_GCAMORPH_LABELS : gcamorph label data TAG_GCAMORPH_AFFINE : gcamorph m_affine matrix The version number in mgz is now intent encoded: version = (intent & 0xff ) << 8) | MGH_VERSION MGH_VERSION = 1 Changes included in this commit: 1. add read/write TAG_GCAMORPH_GEOM and TAG_GCAMORPH_META in Surfa MGHArrayIO class 2. retrieve intent when MGHArrayIO::load(), encode intent in version number when MGHArrayIO:save() 3. mgz warp related information is saved as the following fields in the FramedArray _metadata dict: intent, gcamorph_volgeom_src, gcamorph_volgeom_trg, warpfield_dtfmt, gcamorph_spacing, gcamorph_exp_k --- surfa/core/framed.py | 17 ++++++++++-- surfa/io/framed.py | 50 ++++++++++++++++++++++++++++++---- surfa/io/fsio.py | 13 +++++++++ surfa/io/utils.py | 65 ++++++++++++++++++++++++++++++++++++++++++++ 4 files changed, 137 insertions(+), 8 deletions(-) diff --git a/surfa/core/framed.py b/surfa/core/framed.py index 89400ba..a0b4d98 100644 --- a/surfa/core/framed.py +++ b/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): """ diff --git a/surfa/io/framed.py b/surfa/io/framed.py index 2564f0b..19bd311 100644 --- a/surfa/io/framed.py +++ b/surfa/io/framed.py @@ -9,6 +9,7 @@ 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 @@ -16,6 +17,8 @@ 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)) diff --git a/surfa/io/fsio.py b/surfa/io/fsio.py index 23d6eb7..fe598c2 100644 --- a/surfa/io/fsio.py +++ b/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 ] diff --git a/surfa/io/utils.py b/surfa/io/utils.py index 1cf30ab..78d42e1 100644 --- a/surfa/io/utils.py +++ b/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')) +