From 49bf4c36306e682b3b2889afcdd78b0bf4f9b4a8 Mon Sep 17 00:00:00 2001 From: Julien Marabotto Date: Tue, 7 May 2024 12:19:30 +0200 Subject: [PATCH 01/12] X5.py created --- nitransforms/io/__init__.py | 3 ++- nitransforms/io/x5.py | 53 +++++++++++++++++++++++++++++++++++++ 2 files changed, 55 insertions(+), 1 deletion(-) create mode 100644 nitransforms/io/x5.py diff --git a/nitransforms/io/__init__.py b/nitransforms/io/__init__.py index c38d11c2..f456a6d1 100644 --- a/nitransforms/io/__init__.py +++ b/nitransforms/io/__init__.py @@ -1,7 +1,7 @@ # emacs: -*- mode: python-mode; py-indent-offset: 4; indent-tabs-mode: nil -*- # vi: set ft=python sts=4 ts=4 sw=4 et: """Read and write transforms.""" -from nitransforms.io import afni, fsl, itk, lta +from nitransforms.io import afni, fsl, itk, lta, x5 from nitransforms.io.base import TransformIOError, TransformFileError __all__ = [ @@ -22,6 +22,7 @@ "fs": (lta, "FSLinearTransform"), "fsl": (fsl, "FSLLinearTransform"), "afni": (afni, "AFNILinearTransform"), + "x5": (x5, "X5LinearTransform") } diff --git a/nitransforms/io/x5.py b/nitransforms/io/x5.py new file mode 100644 index 00000000..b38ef9ab --- /dev/null +++ b/nitransforms/io/x5.py @@ -0,0 +1,53 @@ +"""Read/write x5 transforms.""" +import warnings +import numpy as np +from scipy.io import loadmat as _read_mat, savemat as _save_mat +from h5py import File as H5File +from nibabel import Nifti1Header, Nifti1Image +from nibabel.affines import from_matvec +from nitransforms.io.base import ( + BaseLinearTransformList, + DisplacementsField, + LinearParameters, + TransformIOError, + TransformFileError, +) + +LPS = np.diag([-1, -1, 1, 1]) + +class X5LinearTransform(LinearParameters): + """A string-based structure for X5 linear transforms.""" + + template_dtype = np.dtype( + [ + ("type", "i4"), + ("index", "i4"), + ("parameters", "f8", (4, 4)), + ("offset", "f4", 3), # Center of rotation + ] + ) + dtype = template_dtype + + def __init__(self, parameters=None, offset=None): + return + + def __str__(self): + return + + def to_filename(self, filename): + """Store this transform to a file with the appropriate format.""" + sa = self.structarr + affine = np.array( + np.hstack( + (sa["parameters"][:3, :3].reshape(-1), sa["parameters"][:3, 3]) + )[..., np.newaxis], + dtype="f8", + ) + return + + @classmethod + def from_filename(cls, filename): + """Read the struct from a X5 file given its path.""" + if str(filename).endswith(".h5"): + with H5File(str(filename)) as f: + return cls.from_h5obj(f) From 74090ac7a7ea492b5ebba848f67996dcc587ad82 Mon Sep 17 00:00:00 2001 From: Julien Marabotto Date: Fri, 10 May 2024 09:32:23 +0200 Subject: [PATCH 02/12] enh: update implementation of X5 --- nitransforms/io/x5.py | 16 ++++++++++++++++ 1 file changed, 16 insertions(+) diff --git a/nitransforms/io/x5.py b/nitransforms/io/x5.py index b38ef9ab..11403c92 100644 --- a/nitransforms/io/x5.py +++ b/nitransforms/io/x5.py @@ -51,3 +51,19 @@ def from_filename(cls, filename): if str(filename).endswith(".h5"): with H5File(str(filename)) as f: return cls.from_h5obj(f) + +class X5LinearTransformArray(BaseLinearTransformList): + """A string-based structure for series of X5 linear transforms.""" + + _inner_type = X5LinearTransform + + @property + def xforms(self): + """Get the list of internal ITKLinearTransforms.""" + return self._xforms + + @xforms.setter + def xforms(self, value): + self._xforms = list(value) + for i, val in enumerate(self.xforms): + val["index"] = i \ No newline at end of file From ba8a389156728c7b300ff0f6b59cfc9a63f50775 Mon Sep 17 00:00:00 2001 From: Julien Marabotto Date: Fri, 10 May 2024 14:49:03 +0200 Subject: [PATCH 03/12] enh: update implementation of X5 --- nitransforms/io/x5.py | 15 ++------------- 1 file changed, 2 insertions(+), 13 deletions(-) diff --git a/nitransforms/io/x5.py b/nitransforms/io/x5.py index 11403c92..264df3d7 100644 --- a/nitransforms/io/x5.py +++ b/nitransforms/io/x5.py @@ -35,14 +35,9 @@ def __str__(self): return def to_filename(self, filename): - """Store this transform to a file with the appropriate format.""" + '''store this transform to a file with the X5 format''' sa = self.structarr - affine = np.array( - np.hstack( - (sa["parameters"][:3, :3].reshape(-1), sa["parameters"][:3, 3]) - )[..., np.newaxis], - dtype="f8", - ) + affine = '''some affine that will return a 4x4 array''' return @classmethod @@ -61,9 +56,3 @@ class X5LinearTransformArray(BaseLinearTransformList): def xforms(self): """Get the list of internal ITKLinearTransforms.""" return self._xforms - - @xforms.setter - def xforms(self, value): - self._xforms = list(value) - for i, val in enumerate(self.xforms): - val["index"] = i \ No newline at end of file From e4e3f44c0951704d2a1a5653aedd9ad6dab64958 Mon Sep 17 00:00:00 2001 From: Julien Marabotto Date: Fri, 17 May 2024 16:17:51 +0200 Subject: [PATCH 04/12] enh: x5 implementation --- nitransforms/base.py | 45 +++++++++++++++++++++++++++++++++++++ nitransforms/io/__init__.py | 2 +- nitransforms/io/x5.py | 41 ++++++++++++++++----------------- nitransforms/linear.py | 27 ++++++++++++++++++---- 4 files changed, 88 insertions(+), 27 deletions(-) diff --git a/nitransforms/base.py b/nitransforms/base.py index 96f00edb..a9e8448d 100644 --- a/nitransforms/base.py +++ b/nitransforms/base.py @@ -350,7 +350,52 @@ def to_filename(self, filename, fmt="X5"): def _to_hdf5(self, x5_root): """Serialize this object into the x5 file format.""" + transform_group = x5_root.create_group("TransformGroup") + + """Group '0' containing Affine transform""" + transform_0 = transform_group.create_group("0") + + transform_0.attrs["Type"] = "Affine" + transform_0.create_dataset("Transform", data=self._matrix) + transform_0.create_dataset("Inverse", data=np.linalg.inv(self._matrix)) + + metadata = {"key": "value"} + transform_0.attrs["Metadata"] = str(metadata) + + """sub-group 'Domain' contained within group '0' """ + domain_group = transform_0.create_group("Domain") + domain_group.attrs["Grid"] = self.grid + domain_group.create_dataset("Size", data=_as_homogeneous(self._reference.shape)) + domain_group.create_dataset("Mapping", data=self.map) + raise NotImplementedError + + def read_x5(x5_root): + variables = {} + with h5py.File(x5_root, "r") as f: + f.visititems(lambda name, x5_root: _from_hdf5(name, x5_root, variables)) + + _transform = variables["TransformGroup/0/Transform"] + _inverse = variables["TransformGroup/0/Inverse"] + _size = variables["TransformGroup/0/Domain/Size"] + _map = variables["TransformGroup/0/Domain/Mapping"] + + return _transform, _inverse, _size, _map + + def _from_hdf5(name, x5_root, storage): + if isinstance(x5_root, h5py.Dataset): + storage[name] = { + 'type': 'dataset', + 'attrs': dict(x5_root.attrs), + 'shape': x5_root.shape, + 'data': x5_root[()] # Read the data + } + elif isinstance(x5_root, h5py.Group): + storage[name] = { + 'type': 'group', + 'attrs': dict(x5_root.attrs), + 'members': {} + } def _as_homogeneous(xyz, dtype="float32", dim=3): diff --git a/nitransforms/io/__init__.py b/nitransforms/io/__init__.py index f456a6d1..3627fd34 100644 --- a/nitransforms/io/__init__.py +++ b/nitransforms/io/__init__.py @@ -22,7 +22,7 @@ "fs": (lta, "FSLinearTransform"), "fsl": (fsl, "FSLLinearTransform"), "afni": (afni, "AFNILinearTransform"), - "x5": (x5, "X5LinearTransform") + "x5": (x5, "X5Transform") } diff --git a/nitransforms/io/x5.py b/nitransforms/io/x5.py index 264df3d7..72b08a0f 100644 --- a/nitransforms/io/x5.py +++ b/nitransforms/io/x5.py @@ -13,20 +13,10 @@ TransformFileError, ) -LPS = np.diag([-1, -1, 1, 1]) - -class X5LinearTransform(LinearParameters): +class X5Transform: """A string-based structure for X5 linear transforms.""" - template_dtype = np.dtype( - [ - ("type", "i4"), - ("index", "i4"), - ("parameters", "f8", (4, 4)), - ("offset", "f4", 3), # Center of rotation - ] - ) - dtype = template_dtype + _transform = None def __init__(self, parameters=None, offset=None): return @@ -34,25 +24,32 @@ def __init__(self, parameters=None, offset=None): def __str__(self): return - def to_filename(self, filename): - '''store this transform to a file with the X5 format''' - sa = self.structarr - affine = '''some affine that will return a 4x4 array''' - return - @classmethod def from_filename(cls, filename): """Read the struct from a X5 file given its path.""" if str(filename).endswith(".h5"): - with H5File(str(filename)) as f: - return cls.from_h5obj(f) + with H5File(str(filename), 'r') as hdf: + return cls.from_h5obj(hdf) + + @classmethod + def from_h5obj(cls, h5obj): + """Read the transformations in an X5 file.""" + xfm_list = list(h5obj.keys()) + + xfm = xfm_list["Transform"] + inv = xfm_list["Inverse"] + coords = xfm_list["Size"] + map = xfm_list["Mapping"] + + return xfm, inv, coords, map + class X5LinearTransformArray(BaseLinearTransformList): """A string-based structure for series of X5 linear transforms.""" - _inner_type = X5LinearTransform + _inner_type = X5Transform @property def xforms(self): - """Get the list of internal ITKLinearTransforms.""" + """Get the list of internal X5LinearTransforms.""" return self._xforms diff --git a/nitransforms/linear.py b/nitransforms/linear.py index af14f396..f7eb5f1b 100644 --- a/nitransforms/linear.py +++ b/nitransforms/linear.py @@ -180,15 +180,34 @@ def map(self, x, inverse=False): affine = self._inverse return affine.dot(coords).T[..., :-1] - def _to_hdf5(self, x5_root): + def _to_hdf5(self, x, x5_root): """Serialize this object into the x5 file format.""" - xform = x5_root.create_dataset("Transform", data=[self._matrix]) - xform.attrs["Type"] = "affine" - x5_root.create_dataset("Inverse", data=[(~self).matrix]) + transgrp = x5_root.create_group("TransformGroup") + affine = self._x5group_affine(transgrp) + coords = self._x5group_domain(x, affine) if self._reference: self.reference._to_hdf5(x5_root.create_group("Reference")) + return #nothing? + + def _x5group_affine(self, TransformGroup): + """Create group "0" for affine in x5_root/TransformGroup/ according to x5 file format""" + aff = TransformGroup.create_group("0") + aff.attrs["Type"] = "affine" #Should have shape {scalar} + aff.attrs["Metadata"] = 'metadata' #This is a draft for metadata. Should have shape {scalar} + aff.create_dataset("Transform", data=[self._matrix]) #Should have shape {3,4} + aff.create_dataset("Inverse", data=[(~self).matrix]) #Should have shape {4,3} + return aff + + def _x5group_domain(self, x, transform): + """Create group "Domain" in x5_root/TransformGroup/0/ according to x5 file format""" + coords = transform.create_group("Domain") + coords.attrs["Grid"] = "grid" #How do I interpet this 'grid'? Should have shape {scalar} + coords.create_dataset("Size", data=_as_homogeneous(x, dim=self._matrix.shape[0] - 1).T) #Should have shape {3} + coords.create_dataset("Mapping", data=[self.map(self, x)]) #Should have shape {4,4} + return coords + def to_filename(self, filename, fmt="X5", moving=None): """Store the transform in the requested output format.""" writer = get_linear_factory(fmt, is_array=False) From f51ca65cec9cd6d90374fd817207f4bedb5a3352 Mon Sep 17 00:00:00 2001 From: Julien Marabotto Date: Thu, 23 May 2024 10:59:09 +0200 Subject: [PATCH 05/12] enh: small changes to implementation of . x5 --- nitransforms/base.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/nitransforms/base.py b/nitransforms/base.py index a9e8448d..ef81c521 100644 --- a/nitransforms/base.py +++ b/nitransforms/base.py @@ -370,10 +370,10 @@ def _to_hdf5(self, x5_root): raise NotImplementedError - def read_x5(x5_root): + def read_x5(self, x5_root): variables = {} with h5py.File(x5_root, "r") as f: - f.visititems(lambda name, x5_root: _from_hdf5(name, x5_root, variables)) + f.visititems(lambda filename, x5_root: self._from_hdf5(filename, x5_root, variables)) _transform = variables["TransformGroup/0/Transform"] _inverse = variables["TransformGroup/0/Inverse"] @@ -382,7 +382,7 @@ def read_x5(x5_root): return _transform, _inverse, _size, _map - def _from_hdf5(name, x5_root, storage): + def _from_hdf5(self, name, x5_root, storage): if isinstance(x5_root, h5py.Dataset): storage[name] = { 'type': 'dataset', From 7ba96c3bb559556caab6c5b9fea423a4dc027068 Mon Sep 17 00:00:00 2001 From: sgiavasis Date: Thu, 18 Jul 2024 00:58:43 -0400 Subject: [PATCH 06/12] Set up I/O functions to be expandable, deferring some design decisions to the future. Started the beginning of some xfm_utils for quicker testing and validation. --- nitransforms/base.py | 105 ++++++++++++++------------------- nitransforms/cli.py | 43 +++++++++++++- nitransforms/io/base.py | 126 ++++++++++++++++++++++++++++++++++++++++ 3 files changed, 209 insertions(+), 65 deletions(-) diff --git a/nitransforms/base.py b/nitransforms/base.py index a9e8448d..321b88a3 100644 --- a/nitransforms/base.py +++ b/nitransforms/base.py @@ -178,14 +178,48 @@ def __ne__(self, other): class TransformBase: """Abstract image class to represent transforms.""" - __slots__ = ("_reference", "_ndim",) - - def __init__(self, reference=None): + __slots__ = ("_reference", "_ndim", "_affine", "_shape", "_header", + "_grid", "_mapping", "_hdf5_dct", "_x5_dct") + + x5_struct = { + 'TransformGroup/0': { + 'Type': None, + 'Transform': None, + 'Metadata': None, + 'Inverse': None + }, + 'TransformGroup/0/Domain': { + 'Grid': None, + 'Size': None, + 'Mapping': None + }, + 'TransformGroup/1': {}, + 'TransformChain': {} + } + + def __init__(self, x5=None, hdf5=None, nifti=None, shape=None, affine=None, + header=None, reference=None): """Instantiate a transform.""" + self._reference = None if reference: self.reference = reference + if nifti is not None: + self._x5_dct = self.init_x5_structure(nifti) + elif hdf5: + self.update_x5_structure(hdf5) + elif x5: + self.update_x5_structure(x5) + + self._shape = shape + self._affine = affine + self._header = header + + # TO-DO + self._grid = None + self._mapping = None + def __call__(self, x, inverse=False): """Apply y = f(x).""" return self.map(x, inverse=inverse) @@ -222,6 +256,12 @@ def ndim(self): """Access the dimensions of the reference space.""" raise TypeError("TransformBase has no dimensions") + def init_x5_structure(self, xfm_data=None): + self.x5_struct['TransformGroup/0/Transform'] = xfm_data + + def update_x5_structure(self, hdf5_struct=None): + self.x5_struct.update(hdf5_struct) + def apply( self, spatialimage, @@ -338,65 +378,6 @@ def map(self, x, inverse=False): """ return x - def to_filename(self, filename, fmt="X5"): - """Store the transform in BIDS-Transforms HDF5 file format (.x5).""" - with h5py.File(filename, "w") as out_file: - out_file.attrs["Format"] = "X5" - out_file.attrs["Version"] = np.uint16(1) - root = out_file.create_group("/0") - self._to_hdf5(root) - - return filename - - def _to_hdf5(self, x5_root): - """Serialize this object into the x5 file format.""" - transform_group = x5_root.create_group("TransformGroup") - - """Group '0' containing Affine transform""" - transform_0 = transform_group.create_group("0") - - transform_0.attrs["Type"] = "Affine" - transform_0.create_dataset("Transform", data=self._matrix) - transform_0.create_dataset("Inverse", data=np.linalg.inv(self._matrix)) - - metadata = {"key": "value"} - transform_0.attrs["Metadata"] = str(metadata) - - """sub-group 'Domain' contained within group '0' """ - domain_group = transform_0.create_group("Domain") - domain_group.attrs["Grid"] = self.grid - domain_group.create_dataset("Size", data=_as_homogeneous(self._reference.shape)) - domain_group.create_dataset("Mapping", data=self.map) - - raise NotImplementedError - - def read_x5(x5_root): - variables = {} - with h5py.File(x5_root, "r") as f: - f.visititems(lambda name, x5_root: _from_hdf5(name, x5_root, variables)) - - _transform = variables["TransformGroup/0/Transform"] - _inverse = variables["TransformGroup/0/Inverse"] - _size = variables["TransformGroup/0/Domain/Size"] - _map = variables["TransformGroup/0/Domain/Mapping"] - - return _transform, _inverse, _size, _map - - def _from_hdf5(name, x5_root, storage): - if isinstance(x5_root, h5py.Dataset): - storage[name] = { - 'type': 'dataset', - 'attrs': dict(x5_root.attrs), - 'shape': x5_root.shape, - 'data': x5_root[()] # Read the data - } - elif isinstance(x5_root, h5py.Group): - storage[name] = { - 'type': 'group', - 'attrs': dict(x5_root.attrs), - 'members': {} - } - def _as_homogeneous(xyz, dtype="float32", dim=3): """ diff --git a/nitransforms/cli.py b/nitransforms/cli.py index 63b8bed4..bcae1541 100644 --- a/nitransforms/cli.py +++ b/nitransforms/cli.py @@ -2,10 +2,13 @@ import os from textwrap import dedent +from nitransforms.io.base import xfm_loader +from nitransforms.linear import load as linload +from nitransforms.nonlinear import load as nlinload -from .linear import load as linload -from .nonlinear import load as nlinload +from nitransforms.base import TransformBase +import pprint def cli_apply(pargs): """ @@ -45,9 +48,27 @@ def cli_apply(pargs): cval=pargs.cval, prefilter=pargs.prefilter, ) - moved.to_filename(pargs.out or f"nt_{os.path.basename(pargs.moving)}") + #moved.to_filename(pargs.out or f"nt_{os.path.basename(pargs.moving)}") +def cli_xfm_util(pargs): + """ + """ + + xfm_data = xfm_loader(pargs.transform) + xfm_x5 = TransformBase(**xfm_data) + + if pargs.info: + pprint.pprint(xfm_x5.x5_struct) + print(f"Shape:\n{xfm_x5._shape}") + print(f"Affine:\n{xfm_x5._affine}") + + if pargs.x5: + filename = f"{os.path.basename(pargs.transform).split('.')[0]}.x5" + xfm_x5.to_filename(filename) + print(f"Writing out {filename}") + + def get_parser(): desc = dedent( """ @@ -56,6 +77,7 @@ def get_parser(): Commands: apply Apply a transformation to an image + xfm_util Assorted transform utilities For command specific information, use 'nt -h'. """ @@ -120,6 +142,17 @@ def _add_subparser(name, description): help="Determines if the image's data array is prefiltered with a spline filter before " "interpolation (default: True)", ) + + xfm_util = _add_subparser("xfm_util", cli_xfm_util.__doc__) + xfm_util.set_defaults(func=cli_xfm_util) + xfm_util.add_argument("transform", help="The transform file") + xfm_util.add_argument("--info", + action="store_true", + help="Get information about the transform") + xfm_util.add_argument("--x5", + action="store_true", + help="Convert transform to .x5 file format.") + return parser, subparsers @@ -133,3 +166,7 @@ def main(pargs=None): subparser = subparsers.choices[pargs.command] subparser.print_help() raise (e) + + +if __name__ == "__main__": + main() diff --git a/nitransforms/io/base.py b/nitransforms/io/base.py index 6d1a7c8e..f50186dd 100644 --- a/nitransforms/io/base.py +++ b/nitransforms/io/base.py @@ -1,11 +1,137 @@ """Read/write linear transforms.""" from pathlib import Path import numpy as np +import nibabel as nb from nibabel import load as loadimg +import h5py + from ..patched import LabeledWrapStruct +def get_xfm_filetype(xfm_file): + path = Path(xfm_file) + ext = path.suffix + if ext == '.gz' and path.name.endswith('.nii.gz'): + return 'nifti' + + file_types = { + '.nii': 'nifti', + '.h5': 'hdf5', + '.x5': 'x5', + '.txt': 'txt', + '.mat': 'txt' + } + return file_types.get(ext, 'unknown') + +def gather_fields(x5=None, hdf5=None, nifti=None, shape=None, affine=None, header=None): + xfm_fields = { + "x5": x5, + "hdf5": hdf5, + "nifti": nifti, + "header": header, + "shape": shape, + "affine": affine + } + return xfm_fields + +def load_nifti(nifti_file): + nifti_xfm = nb.load(nifti_file) + xfm_data = nifti_xfm.get_fdata() + shape = nifti_xfm.shape + affine = nifti_xfm.affine + header = getattr(nifti_xfm, "header", None) + return gather_fields(nifti=xfm_data, shape=shape, affine=affine, header=header) + +def load_hdf5(hdf5_file): + storage = {} + + def get_hdf5_items(name, x5_root): + if isinstance(x5_root, h5py.Dataset): + storage[name] = { + 'type': 'dataset', + 'attrs': dict(x5_root.attrs), + 'shape': x5_root.shape, + 'data': x5_root[()] + } + elif isinstance(x5_root, h5py.Group): + storage[name] = { + 'type': 'group', + 'attrs': dict(x5_root.attrs), + 'members': {} + } + + with h5py.File(hdf5_file, 'r') as f: + f.visititems(get_hdf5_items) + if storage: + hdf5_storage = {'hdf5': storage} + return hdf5_storage + +def load_x5(x5_file): + load_hdf5(x5_file) + +def load_mat(mat_file): + affine_matrix = np.loadtxt(mat_file) + affine = nb.affines.from_matvec(affine_matrix[:,:3], affine_matrix[:,3]) + return gather_fields(affine=affine) + +def xfm_loader(xfm_file): + loaders = { + 'nifti': load_nifti, + 'hdf5': load_hdf5, + 'x5': load_x5, + 'txt': load_mat, + 'mat': load_mat + } + xfm_filetype = get_xfm_filetype(xfm_file) + loader = loaders.get(xfm_filetype) + if loader is None: + raise ValueError(f"Unsupported file type: {xfm_filetype}") + return loader(xfm_file) + +def to_filename(self, filename, fmt="X5"): + """Store the transform in BIDS-Transforms HDF5 file format (.x5).""" + with h5py.File(filename, "w") as out_file: + out_file.attrs["Format"] = "X5" + out_file.attrs["Version"] = np.uint16(1) + root = out_file.create_group("/0") + self._to_hdf5(root) + + return filename + +def _to_hdf5(self, x5_root): + """Serialize this object into the x5 file format.""" + transform_group = x5_root.create_group("TransformGroup") + + """Group '0' containing Affine transform""" + transform_0 = transform_group.create_group("0") + + transform_0.attrs["Type"] = "Affine" + transform_0.create_dataset("Transform", data=self._affine) + transform_0.create_dataset("Inverse", data=np.linalg.inv(self._affine)) + + metadata = {"key": "value"} + transform_0.attrs["Metadata"] = str(metadata) + + """sub-group 'Domain' contained within group '0' """ + domain_group = transform_0.create_group("Domain") + #domain_group.attrs["Grid"] = self._grid + #domain_group.create_dataset("Size", data=_as_homogeneous(self._reference.shape)) + #domain_group.create_dataset("Mapping", data=self.mapping) + +def _from_x5(self, x5_root): + variables = {} + + x5_root.visititems(lambda name, x5_root: loader(name, x5_root, variables)) + + _transform = variables["TransformGroup/0/Transform"] + _inverse = variables["TransformGroup/0/Inverse"] + _size = variables["TransformGroup/0/Domain/Size"] + _mapping = variables["TransformGroup/0/Domain/Mapping"] + + return _transform, _inverse, _size, _map + + class TransformIOError(IOError): """General I/O exception while reading/writing transforms.""" From 0e213cb259030a92af81eb54dc0bd792bed4abf2 Mon Sep 17 00:00:00 2001 From: Julien Marabotto Date: Tue, 7 May 2024 12:19:30 +0200 Subject: [PATCH 07/12] X5.py created --- nitransforms/io/__init__.py | 3 ++- nitransforms/io/x5.py | 53 +++++++++++++++++++++++++++++++++++++ 2 files changed, 55 insertions(+), 1 deletion(-) create mode 100644 nitransforms/io/x5.py diff --git a/nitransforms/io/__init__.py b/nitransforms/io/__init__.py index c38d11c2..f456a6d1 100644 --- a/nitransforms/io/__init__.py +++ b/nitransforms/io/__init__.py @@ -1,7 +1,7 @@ # emacs: -*- mode: python-mode; py-indent-offset: 4; indent-tabs-mode: nil -*- # vi: set ft=python sts=4 ts=4 sw=4 et: """Read and write transforms.""" -from nitransforms.io import afni, fsl, itk, lta +from nitransforms.io import afni, fsl, itk, lta, x5 from nitransforms.io.base import TransformIOError, TransformFileError __all__ = [ @@ -22,6 +22,7 @@ "fs": (lta, "FSLinearTransform"), "fsl": (fsl, "FSLLinearTransform"), "afni": (afni, "AFNILinearTransform"), + "x5": (x5, "X5LinearTransform") } diff --git a/nitransforms/io/x5.py b/nitransforms/io/x5.py new file mode 100644 index 00000000..b38ef9ab --- /dev/null +++ b/nitransforms/io/x5.py @@ -0,0 +1,53 @@ +"""Read/write x5 transforms.""" +import warnings +import numpy as np +from scipy.io import loadmat as _read_mat, savemat as _save_mat +from h5py import File as H5File +from nibabel import Nifti1Header, Nifti1Image +from nibabel.affines import from_matvec +from nitransforms.io.base import ( + BaseLinearTransformList, + DisplacementsField, + LinearParameters, + TransformIOError, + TransformFileError, +) + +LPS = np.diag([-1, -1, 1, 1]) + +class X5LinearTransform(LinearParameters): + """A string-based structure for X5 linear transforms.""" + + template_dtype = np.dtype( + [ + ("type", "i4"), + ("index", "i4"), + ("parameters", "f8", (4, 4)), + ("offset", "f4", 3), # Center of rotation + ] + ) + dtype = template_dtype + + def __init__(self, parameters=None, offset=None): + return + + def __str__(self): + return + + def to_filename(self, filename): + """Store this transform to a file with the appropriate format.""" + sa = self.structarr + affine = np.array( + np.hstack( + (sa["parameters"][:3, :3].reshape(-1), sa["parameters"][:3, 3]) + )[..., np.newaxis], + dtype="f8", + ) + return + + @classmethod + def from_filename(cls, filename): + """Read the struct from a X5 file given its path.""" + if str(filename).endswith(".h5"): + with H5File(str(filename)) as f: + return cls.from_h5obj(f) From 3d263bb36d74e21dc461987f4cf055af22f0dab8 Mon Sep 17 00:00:00 2001 From: Julien Marabotto Date: Fri, 10 May 2024 09:32:23 +0200 Subject: [PATCH 08/12] enh: update implementation of X5 --- nitransforms/io/x5.py | 16 ++++++++++++++++ 1 file changed, 16 insertions(+) diff --git a/nitransforms/io/x5.py b/nitransforms/io/x5.py index b38ef9ab..11403c92 100644 --- a/nitransforms/io/x5.py +++ b/nitransforms/io/x5.py @@ -51,3 +51,19 @@ def from_filename(cls, filename): if str(filename).endswith(".h5"): with H5File(str(filename)) as f: return cls.from_h5obj(f) + +class X5LinearTransformArray(BaseLinearTransformList): + """A string-based structure for series of X5 linear transforms.""" + + _inner_type = X5LinearTransform + + @property + def xforms(self): + """Get the list of internal ITKLinearTransforms.""" + return self._xforms + + @xforms.setter + def xforms(self, value): + self._xforms = list(value) + for i, val in enumerate(self.xforms): + val["index"] = i \ No newline at end of file From a102f1e2bee031c55708dc6dc4e2598c94353a39 Mon Sep 17 00:00:00 2001 From: Julien Marabotto Date: Fri, 10 May 2024 14:49:03 +0200 Subject: [PATCH 09/12] enh: update implementation of X5 --- nitransforms/io/x5.py | 15 ++------------- 1 file changed, 2 insertions(+), 13 deletions(-) diff --git a/nitransforms/io/x5.py b/nitransforms/io/x5.py index 11403c92..264df3d7 100644 --- a/nitransforms/io/x5.py +++ b/nitransforms/io/x5.py @@ -35,14 +35,9 @@ def __str__(self): return def to_filename(self, filename): - """Store this transform to a file with the appropriate format.""" + '''store this transform to a file with the X5 format''' sa = self.structarr - affine = np.array( - np.hstack( - (sa["parameters"][:3, :3].reshape(-1), sa["parameters"][:3, 3]) - )[..., np.newaxis], - dtype="f8", - ) + affine = '''some affine that will return a 4x4 array''' return @classmethod @@ -61,9 +56,3 @@ class X5LinearTransformArray(BaseLinearTransformList): def xforms(self): """Get the list of internal ITKLinearTransforms.""" return self._xforms - - @xforms.setter - def xforms(self, value): - self._xforms = list(value) - for i, val in enumerate(self.xforms): - val["index"] = i \ No newline at end of file From e465332c8ec71b1a05aaa4cbb928b8a287586712 Mon Sep 17 00:00:00 2001 From: Julien Marabotto Date: Fri, 17 May 2024 16:17:51 +0200 Subject: [PATCH 10/12] enh: x5 implementation --- nitransforms/base.py | 63 ++++++++++++------ nitransforms/cli.py | 46 +++++++++++-- nitransforms/io/__init__.py | 2 +- nitransforms/io/base.py | 126 ++++++++++++++++++++++++++++++++++++ nitransforms/io/x5.py | 41 ++++++------ nitransforms/linear.py | 27 ++++++-- 6 files changed, 253 insertions(+), 52 deletions(-) diff --git a/nitransforms/base.py b/nitransforms/base.py index 81ed1a5e..ac6f5c22 100644 --- a/nitransforms/base.py +++ b/nitransforms/base.py @@ -248,17 +248,48 @@ def __ne__(self, other): class TransformBase: """Abstract image class to represent transforms.""" - __slots__ = ( - "_reference", - "_ndim", - ) - - def __init__(self, reference=None): + __slots__ = ("_reference", "_ndim", "_affine", "_shape", "_header", + "_grid", "_mapping", "_hdf5_dct", "_x5_dct") + + x5_struct = { + 'TransformGroup/0': { + 'Type': None, + 'Transform': None, + 'Metadata': None, + 'Inverse': None + }, + 'TransformGroup/0/Domain': { + 'Grid': None, + 'Size': None, + 'Mapping': None + }, + 'TransformGroup/1': {}, + 'TransformChain': {} + } + + def __init__(self, x5=None, hdf5=None, nifti=None, shape=None, affine=None, + header=None, reference=None): """Instantiate a transform.""" + self._reference = None if reference: self.reference = reference + if nifti is not None: + self._x5_dct = self.init_x5_structure(nifti) + elif hdf5: + self.update_x5_structure(hdf5) + elif x5: + self.update_x5_structure(x5) + + self._shape = shape + self._affine = affine + self._header = header + + # TO-DO + self._grid = None + self._mapping = None + def __call__(self, x, inverse=False): """Apply y = f(x).""" return self.map(x, inverse=inverse) @@ -295,6 +326,12 @@ def ndim(self): """Access the dimensions of the reference space.""" raise TypeError("TransformBase has no dimensions") + def init_x5_structure(self, xfm_data=None): + self.x5_struct['TransformGroup/0/Transform'] = xfm_data + + def update_x5_structure(self, hdf5_struct=None): + self.x5_struct.update(hdf5_struct) + def map(self, x, inverse=False): r""" Apply :math:`y = f(x)`. @@ -316,20 +353,6 @@ def map(self, x, inverse=False): """ return x - def to_filename(self, filename, fmt="X5"): - """Store the transform in BIDS-Transforms HDF5 file format (.x5).""" - with h5py.File(filename, "w") as out_file: - out_file.attrs["Format"] = "X5" - out_file.attrs["Version"] = np.uint16(1) - root = out_file.create_group("/0") - self._to_hdf5(root) - - return filename - - def _to_hdf5(self, x5_root): - """Serialize this object into the x5 file format.""" - raise NotImplementedError - def apply(self, *args, **kwargs): """Apply the transform to a dataset. diff --git a/nitransforms/cli.py b/nitransforms/cli.py index 8f8f5ce0..3ade8486 100644 --- a/nitransforms/cli.py +++ b/nitransforms/cli.py @@ -2,11 +2,13 @@ import os from textwrap import dedent +from nitransforms.base import TransformBase +from nitransforms.io.base import xfm_loader +from nitransforms.linear import load as linload +from nitransforms.nonlinear import load as nlinload +from nitransforms.resampling import apply -from .linear import load as linload -from .nonlinear import load as nlinload -from .resampling import apply - +import pprint def cli_apply(pargs): """ @@ -47,8 +49,26 @@ def cli_apply(pargs): cval=pargs.cval, prefilter=pargs.prefilter, ) - moved.to_filename(pargs.out or f"nt_{os.path.basename(pargs.moving)}") + #moved.to_filename(pargs.out or f"nt_{os.path.basename(pargs.moving)}") + + +def cli_xfm_util(pargs): + """ + """ + + xfm_data = xfm_loader(pargs.transform) + xfm_x5 = TransformBase(**xfm_data) + + if pargs.info: + pprint.pprint(xfm_x5.x5_struct) + print(f"Shape:\n{xfm_x5._shape}") + print(f"Affine:\n{xfm_x5._affine}") + if pargs.x5: + filename = f"{os.path.basename(pargs.transform).split('.')[0]}.x5" + xfm_x5.to_filename(filename) + print(f"Writing out {filename}") + def get_parser(): desc = dedent( @@ -58,6 +78,7 @@ def get_parser(): Commands: apply Apply a transformation to an image + xfm_util Assorted transform utilities For command specific information, use 'nt -h'. """ @@ -122,6 +143,17 @@ def _add_subparser(name, description): help="Determines if the image's data array is prefiltered with a spline filter before " "interpolation (default: True)", ) + + xfm_util = _add_subparser("xfm_util", cli_xfm_util.__doc__) + xfm_util.set_defaults(func=cli_xfm_util) + xfm_util.add_argument("transform", help="The transform file") + xfm_util.add_argument("--info", + action="store_true", + help="Get information about the transform") + xfm_util.add_argument("--x5", + action="store_true", + help="Convert transform to .x5 file format.") + return parser, subparsers @@ -135,3 +167,7 @@ def main(pargs=None): subparser = subparsers.choices[pargs.command] subparser.print_help() raise (e) + + +if __name__ == "__main__": + main() diff --git a/nitransforms/io/__init__.py b/nitransforms/io/__init__.py index f456a6d1..3627fd34 100644 --- a/nitransforms/io/__init__.py +++ b/nitransforms/io/__init__.py @@ -22,7 +22,7 @@ "fs": (lta, "FSLinearTransform"), "fsl": (fsl, "FSLLinearTransform"), "afni": (afni, "AFNILinearTransform"), - "x5": (x5, "X5LinearTransform") + "x5": (x5, "X5Transform") } diff --git a/nitransforms/io/base.py b/nitransforms/io/base.py index 3c923426..78c28827 100644 --- a/nitransforms/io/base.py +++ b/nitransforms/io/base.py @@ -1,11 +1,137 @@ """Read/write linear transforms.""" from pathlib import Path import numpy as np +import nibabel as nb from nibabel import load as loadimg +import h5py + from ..patched import LabeledWrapStruct +def get_xfm_filetype(xfm_file): + path = Path(xfm_file) + ext = path.suffix + if ext == '.gz' and path.name.endswith('.nii.gz'): + return 'nifti' + + file_types = { + '.nii': 'nifti', + '.h5': 'hdf5', + '.x5': 'x5', + '.txt': 'txt', + '.mat': 'txt' + } + return file_types.get(ext, 'unknown') + +def gather_fields(x5=None, hdf5=None, nifti=None, shape=None, affine=None, header=None): + xfm_fields = { + "x5": x5, + "hdf5": hdf5, + "nifti": nifti, + "header": header, + "shape": shape, + "affine": affine + } + return xfm_fields + +def load_nifti(nifti_file): + nifti_xfm = nb.load(nifti_file) + xfm_data = nifti_xfm.get_fdata() + shape = nifti_xfm.shape + affine = nifti_xfm.affine + header = getattr(nifti_xfm, "header", None) + return gather_fields(nifti=xfm_data, shape=shape, affine=affine, header=header) + +def load_hdf5(hdf5_file): + storage = {} + + def get_hdf5_items(name, x5_root): + if isinstance(x5_root, h5py.Dataset): + storage[name] = { + 'type': 'dataset', + 'attrs': dict(x5_root.attrs), + 'shape': x5_root.shape, + 'data': x5_root[()] + } + elif isinstance(x5_root, h5py.Group): + storage[name] = { + 'type': 'group', + 'attrs': dict(x5_root.attrs), + 'members': {} + } + + with h5py.File(hdf5_file, 'r') as f: + f.visititems(get_hdf5_items) + if storage: + hdf5_storage = {'hdf5': storage} + return hdf5_storage + +def load_x5(x5_file): + load_hdf5(x5_file) + +def load_mat(mat_file): + affine_matrix = np.loadtxt(mat_file) + affine = nb.affines.from_matvec(affine_matrix[:,:3], affine_matrix[:,3]) + return gather_fields(affine=affine) + +def xfm_loader(xfm_file): + loaders = { + 'nifti': load_nifti, + 'hdf5': load_hdf5, + 'x5': load_x5, + 'txt': load_mat, + 'mat': load_mat + } + xfm_filetype = get_xfm_filetype(xfm_file) + loader = loaders.get(xfm_filetype) + if loader is None: + raise ValueError(f"Unsupported file type: {xfm_filetype}") + return loader(xfm_file) + +def to_filename(self, filename, fmt="X5"): + """Store the transform in BIDS-Transforms HDF5 file format (.x5).""" + with h5py.File(filename, "w") as out_file: + out_file.attrs["Format"] = "X5" + out_file.attrs["Version"] = np.uint16(1) + root = out_file.create_group("/0") + self._to_hdf5(root) + + return filename + +def _to_hdf5(self, x5_root): + """Serialize this object into the x5 file format.""" + transform_group = x5_root.create_group("TransformGroup") + + """Group '0' containing Affine transform""" + transform_0 = transform_group.create_group("0") + + transform_0.attrs["Type"] = "Affine" + transform_0.create_dataset("Transform", data=self._affine) + transform_0.create_dataset("Inverse", data=np.linalg.inv(self._affine)) + + metadata = {"key": "value"} + transform_0.attrs["Metadata"] = str(metadata) + + """sub-group 'Domain' contained within group '0' """ + domain_group = transform_0.create_group("Domain") + #domain_group.attrs["Grid"] = self._grid + #domain_group.create_dataset("Size", data=_as_homogeneous(self._reference.shape)) + #domain_group.create_dataset("Mapping", data=self.mapping) + +def _from_x5(self, x5_root): + variables = {} + + x5_root.visititems(lambda name, x5_root: loader(name, x5_root, variables)) + + _transform = variables["TransformGroup/0/Transform"] + _inverse = variables["TransformGroup/0/Inverse"] + _size = variables["TransformGroup/0/Domain/Size"] + _mapping = variables["TransformGroup/0/Domain/Mapping"] + + return _transform, _inverse, _size, _map + + class TransformIOError(IOError): """General I/O exception while reading/writing transforms.""" diff --git a/nitransforms/io/x5.py b/nitransforms/io/x5.py index 264df3d7..72b08a0f 100644 --- a/nitransforms/io/x5.py +++ b/nitransforms/io/x5.py @@ -13,20 +13,10 @@ TransformFileError, ) -LPS = np.diag([-1, -1, 1, 1]) - -class X5LinearTransform(LinearParameters): +class X5Transform: """A string-based structure for X5 linear transforms.""" - template_dtype = np.dtype( - [ - ("type", "i4"), - ("index", "i4"), - ("parameters", "f8", (4, 4)), - ("offset", "f4", 3), # Center of rotation - ] - ) - dtype = template_dtype + _transform = None def __init__(self, parameters=None, offset=None): return @@ -34,25 +24,32 @@ def __init__(self, parameters=None, offset=None): def __str__(self): return - def to_filename(self, filename): - '''store this transform to a file with the X5 format''' - sa = self.structarr - affine = '''some affine that will return a 4x4 array''' - return - @classmethod def from_filename(cls, filename): """Read the struct from a X5 file given its path.""" if str(filename).endswith(".h5"): - with H5File(str(filename)) as f: - return cls.from_h5obj(f) + with H5File(str(filename), 'r') as hdf: + return cls.from_h5obj(hdf) + + @classmethod + def from_h5obj(cls, h5obj): + """Read the transformations in an X5 file.""" + xfm_list = list(h5obj.keys()) + + xfm = xfm_list["Transform"] + inv = xfm_list["Inverse"] + coords = xfm_list["Size"] + map = xfm_list["Mapping"] + + return xfm, inv, coords, map + class X5LinearTransformArray(BaseLinearTransformList): """A string-based structure for series of X5 linear transforms.""" - _inner_type = X5LinearTransform + _inner_type = X5Transform @property def xforms(self): - """Get the list of internal ITKLinearTransforms.""" + """Get the list of internal X5LinearTransforms.""" return self._xforms diff --git a/nitransforms/linear.py b/nitransforms/linear.py index 71df6a16..d8597dd1 100644 --- a/nitransforms/linear.py +++ b/nitransforms/linear.py @@ -180,15 +180,34 @@ def map(self, x, inverse=False): affine = self._inverse return affine.dot(coords).T[..., :-1] - def _to_hdf5(self, x5_root): + def _to_hdf5(self, x, x5_root): """Serialize this object into the x5 file format.""" - xform = x5_root.create_dataset("Transform", data=[self._matrix]) - xform.attrs["Type"] = "affine" - x5_root.create_dataset("Inverse", data=[(~self).matrix]) + transgrp = x5_root.create_group("TransformGroup") + affine = self._x5group_affine(transgrp) + coords = self._x5group_domain(x, affine) if self._reference: self.reference._to_hdf5(x5_root.create_group("Reference")) + return #nothing? + + def _x5group_affine(self, TransformGroup): + """Create group "0" for affine in x5_root/TransformGroup/ according to x5 file format""" + aff = TransformGroup.create_group("0") + aff.attrs["Type"] = "affine" #Should have shape {scalar} + aff.attrs["Metadata"] = 'metadata' #This is a draft for metadata. Should have shape {scalar} + aff.create_dataset("Transform", data=[self._matrix]) #Should have shape {3,4} + aff.create_dataset("Inverse", data=[(~self).matrix]) #Should have shape {4,3} + return aff + + def _x5group_domain(self, x, transform): + """Create group "Domain" in x5_root/TransformGroup/0/ according to x5 file format""" + coords = transform.create_group("Domain") + coords.attrs["Grid"] = "grid" #How do I interpet this 'grid'? Should have shape {scalar} + coords.create_dataset("Size", data=_as_homogeneous(x, dim=self._matrix.shape[0] - 1).T) #Should have shape {3} + coords.create_dataset("Mapping", data=[self.map(self, x)]) #Should have shape {4,4} + return coords + def to_filename(self, filename, fmt="X5", moving=None): """Store the transform in the requested output format.""" writer = get_linear_factory(fmt, is_array=False) From dd399a75772bcb0430f811b5792629bf96a4107a Mon Sep 17 00:00:00 2001 From: Julien Marabotto Date: Thu, 23 May 2024 10:59:09 +0200 Subject: [PATCH 11/12] enh: small changes to implementation of . x5 --- nitransforms/base.py | 49 ++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 49 insertions(+) diff --git a/nitransforms/base.py b/nitransforms/base.py index ac6f5c22..5a7fa005 100644 --- a/nitransforms/base.py +++ b/nitransforms/base.py @@ -366,6 +366,55 @@ def apply(self, *args, **kwargs): return apply(self, *args, **kwargs) + def _to_hdf5(self, x5_root): + """Serialize this object into the x5 file format.""" + transform_group = x5_root.create_group("TransformGroup") + + """Group '0' containing Affine transform""" + transform_0 = transform_group.create_group("0") + + transform_0.attrs["Type"] = "Affine" + transform_0.create_dataset("Transform", data=self._matrix) + transform_0.create_dataset("Inverse", data=np.linalg.inv(self._matrix)) + + metadata = {"key": "value"} + transform_0.attrs["Metadata"] = str(metadata) + + """sub-group 'Domain' contained within group '0' """ + domain_group = transform_0.create_group("Domain") + domain_group.attrs["Grid"] = self.grid + domain_group.create_dataset("Size", data=_as_homogeneous(self._reference.shape)) + domain_group.create_dataset("Mapping", data=self.map) + + raise NotImplementedError + + def read_x5(self, x5_root): + variables = {} + with h5py.File(x5_root, "r") as f: + f.visititems(lambda filename, x5_root: self._from_hdf5(filename, x5_root, variables)) + + _transform = variables["TransformGroup/0/Transform"] + _inverse = variables["TransformGroup/0/Inverse"] + _size = variables["TransformGroup/0/Domain/Size"] + _map = variables["TransformGroup/0/Domain/Mapping"] + + return _transform, _inverse, _size, _map + + def _from_hdf5(self, name, x5_root, storage): + if isinstance(x5_root, h5py.Dataset): + storage[name] = { + 'type': 'dataset', + 'attrs': dict(x5_root.attrs), + 'shape': x5_root.shape, + 'data': x5_root[()] # Read the data + } + elif isinstance(x5_root, h5py.Group): + storage[name] = { + 'type': 'group', + 'attrs': dict(x5_root.attrs), + 'members': {} + } + def _as_homogeneous(xyz, dtype="float32", dim=3): """ From 19bc3531568294014f3e25554ba17db580eec202 Mon Sep 17 00:00:00 2001 From: Oscar Esteban Date: Fri, 2 Aug 2024 10:39:31 +0200 Subject: [PATCH 12/12] sty: run ruff --- nitransforms/base.py | 92 +++++++++++++++++++++---------------- nitransforms/cli.py | 24 +++++----- nitransforms/io/__init__.py | 3 +- nitransforms/io/base.py | 74 ++++++++++++++++------------- nitransforms/io/x5.py | 15 ++---- nitransforms/linear.py | 25 ++++++---- 6 files changed, 130 insertions(+), 103 deletions(-) diff --git a/nitransforms/base.py b/nitransforms/base.py index 5a7fa005..c433898d 100644 --- a/nitransforms/base.py +++ b/nitransforms/base.py @@ -7,6 +7,7 @@ # ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ## """Common interface for transforms.""" + from pathlib import Path import numpy as np import h5py @@ -146,13 +147,13 @@ def from_arrays(cls, coordinates, triangles): darrays = [ nb.gifti.GiftiDataArray( coordinates.astype(np.float32), - intent=nb.nifti1.intent_codes['NIFTI_INTENT_POINTSET'], - datatype=nb.nifti1.data_type_codes['NIFTI_TYPE_FLOAT32'], + intent=nb.nifti1.intent_codes["NIFTI_INTENT_POINTSET"], + datatype=nb.nifti1.data_type_codes["NIFTI_TYPE_FLOAT32"], ), nb.gifti.GiftiDataArray( triangles.astype(np.int32), - intent=nb.nifti1.intent_codes['NIFTI_INTENT_TRIANGLE'], - datatype=nb.nifti1.data_type_codes['NIFTI_TYPE_INT32'], + intent=nb.nifti1.intent_codes["NIFTI_INTENT_TRIANGLE"], + datatype=nb.nifti1.data_type_codes["NIFTI_TYPE_INT32"], ), ] gii = nb.gifti.GiftiImage(darrays=darrays) @@ -248,27 +249,40 @@ def __ne__(self, other): class TransformBase: """Abstract image class to represent transforms.""" - __slots__ = ("_reference", "_ndim", "_affine", "_shape", "_header", - "_grid", "_mapping", "_hdf5_dct", "_x5_dct") + __slots__ = ( + "_reference", + "_ndim", + "_affine", + "_shape", + "_header", + "_grid", + "_mapping", + "_hdf5_dct", + "_x5_dct", + ) x5_struct = { - 'TransformGroup/0': { - 'Type': None, - 'Transform': None, - 'Metadata': None, - 'Inverse': None - }, - 'TransformGroup/0/Domain': { - 'Grid': None, - 'Size': None, - 'Mapping': None + "TransformGroup/0": { + "Type": None, + "Transform": None, + "Metadata": None, + "Inverse": None, }, - 'TransformGroup/1': {}, - 'TransformChain': {} + "TransformGroup/0/Domain": {"Grid": None, "Size": None, "Mapping": None}, + "TransformGroup/1": {}, + "TransformChain": {}, } - def __init__(self, x5=None, hdf5=None, nifti=None, shape=None, affine=None, - header=None, reference=None): + def __init__( + self, + x5=None, + hdf5=None, + nifti=None, + shape=None, + affine=None, + header=None, + reference=None, + ): """Instantiate a transform.""" self._reference = None @@ -281,7 +295,7 @@ def __init__(self, x5=None, hdf5=None, nifti=None, shape=None, affine=None, self.update_x5_structure(hdf5) elif x5: self.update_x5_structure(x5) - + self._shape = shape self._affine = affine self._header = header @@ -327,8 +341,8 @@ def ndim(self): raise TypeError("TransformBase has no dimensions") def init_x5_structure(self, xfm_data=None): - self.x5_struct['TransformGroup/0/Transform'] = xfm_data - + self.x5_struct["TransformGroup/0/Transform"] = xfm_data + def update_x5_structure(self, hdf5_struct=None): self.x5_struct.update(hdf5_struct) @@ -358,9 +372,7 @@ def apply(self, *args, **kwargs): Deprecated. Please use ``nitransforms.resampling.apply`` instead. """ - message = ( - "The `apply` method is deprecated. Please use `nitransforms.resampling.apply` instead." - ) + message = "The `apply` method is deprecated. Please use `nitransforms.resampling.apply` instead." warnings.warn(message, DeprecationWarning, stacklevel=2) from .resampling import apply @@ -372,7 +384,7 @@ def _to_hdf5(self, x5_root): """Group '0' containing Affine transform""" transform_0 = transform_group.create_group("0") - + transform_0.attrs["Type"] = "Affine" transform_0.create_dataset("Transform", data=self._matrix) transform_0.create_dataset("Inverse", data=np.linalg.inv(self._matrix)) @@ -385,13 +397,15 @@ def _to_hdf5(self, x5_root): domain_group.attrs["Grid"] = self.grid domain_group.create_dataset("Size", data=_as_homogeneous(self._reference.shape)) domain_group.create_dataset("Mapping", data=self.map) - + raise NotImplementedError - + def read_x5(self, x5_root): variables = {} with h5py.File(x5_root, "r") as f: - f.visititems(lambda filename, x5_root: self._from_hdf5(filename, x5_root, variables)) + f.visititems( + lambda filename, x5_root: self._from_hdf5(filename, x5_root, variables) + ) _transform = variables["TransformGroup/0/Transform"] _inverse = variables["TransformGroup/0/Inverse"] @@ -399,20 +413,20 @@ def read_x5(self, x5_root): _map = variables["TransformGroup/0/Domain/Mapping"] return _transform, _inverse, _size, _map - + def _from_hdf5(self, name, x5_root, storage): if isinstance(x5_root, h5py.Dataset): storage[name] = { - 'type': 'dataset', - 'attrs': dict(x5_root.attrs), - 'shape': x5_root.shape, - 'data': x5_root[()] # Read the data - } + "type": "dataset", + "attrs": dict(x5_root.attrs), + "shape": x5_root.shape, + "data": x5_root[()], # Read the data + } elif isinstance(x5_root, h5py.Group): storage[name] = { - 'type': 'group', - 'attrs': dict(x5_root.attrs), - 'members': {} + "type": "group", + "attrs": dict(x5_root.attrs), + "members": {}, } diff --git a/nitransforms/cli.py b/nitransforms/cli.py index 3ade8486..1bd21d8f 100644 --- a/nitransforms/cli.py +++ b/nitransforms/cli.py @@ -10,6 +10,7 @@ import pprint + def cli_apply(pargs): """ Apply a transformation to an image, resampling on the reference. @@ -34,8 +35,8 @@ def cli_apply(pargs): xfm = ( nlinload(pargs.transform, fmt=fmt) - if pargs.nonlinear else - linload(pargs.transform, fmt=fmt) + if pargs.nonlinear + else linload(pargs.transform, fmt=fmt) ) # ensure a reference is set @@ -49,12 +50,11 @@ def cli_apply(pargs): cval=pargs.cval, prefilter=pargs.prefilter, ) - #moved.to_filename(pargs.out or f"nt_{os.path.basename(pargs.moving)}") + # moved.to_filename(pargs.out or f"nt_{os.path.basename(pargs.moving)}") def cli_xfm_util(pargs): - """ - """ + """ """ xfm_data = xfm_loader(pargs.transform) xfm_x5 = TransformBase(**xfm_data) @@ -68,7 +68,7 @@ def cli_xfm_util(pargs): filename = f"{os.path.basename(pargs.transform).split('.')[0]}.x5" xfm_x5.to_filename(filename) print(f"Writing out {filename}") - + def get_parser(): desc = dedent( @@ -147,12 +147,12 @@ def _add_subparser(name, description): xfm_util = _add_subparser("xfm_util", cli_xfm_util.__doc__) xfm_util.set_defaults(func=cli_xfm_util) xfm_util.add_argument("transform", help="The transform file") - xfm_util.add_argument("--info", - action="store_true", - help="Get information about the transform") - xfm_util.add_argument("--x5", - action="store_true", - help="Convert transform to .x5 file format.") + xfm_util.add_argument( + "--info", action="store_true", help="Get information about the transform" + ) + xfm_util.add_argument( + "--x5", action="store_true", help="Convert transform to .x5 file format." + ) return parser, subparsers diff --git a/nitransforms/io/__init__.py b/nitransforms/io/__init__.py index 3627fd34..a85e377a 100644 --- a/nitransforms/io/__init__.py +++ b/nitransforms/io/__init__.py @@ -1,6 +1,7 @@ # emacs: -*- mode: python-mode; py-indent-offset: 4; indent-tabs-mode: nil -*- # vi: set ft=python sts=4 ts=4 sw=4 et: """Read and write transforms.""" + from nitransforms.io import afni, fsl, itk, lta, x5 from nitransforms.io.base import TransformIOError, TransformFileError @@ -22,7 +23,7 @@ "fs": (lta, "FSLinearTransform"), "fsl": (fsl, "FSLLinearTransform"), "afni": (afni, "AFNILinearTransform"), - "x5": (x5, "X5Transform") + "x5": (x5, "X5Transform"), } diff --git a/nitransforms/io/base.py b/nitransforms/io/base.py index 78c28827..2ce991bd 100644 --- a/nitransforms/io/base.py +++ b/nitransforms/io/base.py @@ -1,4 +1,5 @@ """Read/write linear transforms.""" + from pathlib import Path import numpy as np import nibabel as nb @@ -12,17 +13,18 @@ def get_xfm_filetype(xfm_file): path = Path(xfm_file) ext = path.suffix - if ext == '.gz' and path.name.endswith('.nii.gz'): - return 'nifti' - + if ext == ".gz" and path.name.endswith(".nii.gz"): + return "nifti" + file_types = { - '.nii': 'nifti', - '.h5': 'hdf5', - '.x5': 'x5', - '.txt': 'txt', - '.mat': 'txt' + ".nii": "nifti", + ".h5": "hdf5", + ".x5": "x5", + ".txt": "txt", + ".mat": "txt", } - return file_types.get(ext, 'unknown') + return file_types.get(ext, "unknown") + def gather_fields(x5=None, hdf5=None, nifti=None, shape=None, affine=None, header=None): xfm_fields = { @@ -31,10 +33,11 @@ def gather_fields(x5=None, hdf5=None, nifti=None, shape=None, affine=None, heade "nifti": nifti, "header": header, "shape": shape, - "affine": affine + "affine": affine, } return xfm_fields + def load_nifti(nifti_file): nifti_xfm = nb.load(nifti_file) xfm_data = nifti_xfm.get_fdata() @@ -43,45 +46,49 @@ def load_nifti(nifti_file): header = getattr(nifti_xfm, "header", None) return gather_fields(nifti=xfm_data, shape=shape, affine=affine, header=header) + def load_hdf5(hdf5_file): storage = {} def get_hdf5_items(name, x5_root): if isinstance(x5_root, h5py.Dataset): storage[name] = { - 'type': 'dataset', - 'attrs': dict(x5_root.attrs), - 'shape': x5_root.shape, - 'data': x5_root[()] - } + "type": "dataset", + "attrs": dict(x5_root.attrs), + "shape": x5_root.shape, + "data": x5_root[()], + } elif isinstance(x5_root, h5py.Group): storage[name] = { - 'type': 'group', - 'attrs': dict(x5_root.attrs), - 'members': {} + "type": "group", + "attrs": dict(x5_root.attrs), + "members": {}, } - with h5py.File(hdf5_file, 'r') as f: + with h5py.File(hdf5_file, "r") as f: f.visititems(get_hdf5_items) if storage: - hdf5_storage = {'hdf5': storage} + hdf5_storage = {"hdf5": storage} return hdf5_storage + def load_x5(x5_file): load_hdf5(x5_file) + def load_mat(mat_file): affine_matrix = np.loadtxt(mat_file) - affine = nb.affines.from_matvec(affine_matrix[:,:3], affine_matrix[:,3]) + affine = nb.affines.from_matvec(affine_matrix[:, :3], affine_matrix[:, 3]) return gather_fields(affine=affine) + def xfm_loader(xfm_file): loaders = { - 'nifti': load_nifti, - 'hdf5': load_hdf5, - 'x5': load_x5, - 'txt': load_mat, - 'mat': load_mat + "nifti": load_nifti, + "hdf5": load_hdf5, + "x5": load_x5, + "txt": load_mat, + "mat": load_mat, } xfm_filetype = get_xfm_filetype(xfm_file) loader = loaders.get(xfm_filetype) @@ -89,6 +96,7 @@ def xfm_loader(xfm_file): raise ValueError(f"Unsupported file type: {xfm_filetype}") return loader(xfm_file) + def to_filename(self, filename, fmt="X5"): """Store the transform in BIDS-Transforms HDF5 file format (.x5).""" with h5py.File(filename, "w") as out_file: @@ -99,13 +107,14 @@ def to_filename(self, filename, fmt="X5"): return filename + def _to_hdf5(self, x5_root): """Serialize this object into the x5 file format.""" transform_group = x5_root.create_group("TransformGroup") """Group '0' containing Affine transform""" transform_0 = transform_group.create_group("0") - + transform_0.attrs["Type"] = "Affine" transform_0.create_dataset("Transform", data=self._affine) transform_0.create_dataset("Inverse", data=np.linalg.inv(self._affine)) @@ -115,10 +124,11 @@ def _to_hdf5(self, x5_root): """sub-group 'Domain' contained within group '0' """ domain_group = transform_0.create_group("Domain") - #domain_group.attrs["Grid"] = self._grid - #domain_group.create_dataset("Size", data=_as_homogeneous(self._reference.shape)) - #domain_group.create_dataset("Mapping", data=self.mapping) - + # domain_group.attrs["Grid"] = self._grid + # domain_group.create_dataset("Size", data=_as_homogeneous(self._reference.shape)) + # domain_group.create_dataset("Mapping", data=self.mapping) + + def _from_x5(self, x5_root): variables = {} @@ -130,7 +140,7 @@ def _from_x5(self, x5_root): _mapping = variables["TransformGroup/0/Domain/Mapping"] return _transform, _inverse, _size, _map - + class TransformIOError(IOError): """General I/O exception while reading/writing transforms.""" diff --git a/nitransforms/io/x5.py b/nitransforms/io/x5.py index 72b08a0f..24331db7 100644 --- a/nitransforms/io/x5.py +++ b/nitransforms/io/x5.py @@ -1,18 +1,11 @@ """Read/write x5 transforms.""" -import warnings -import numpy as np -from scipy.io import loadmat as _read_mat, savemat as _save_mat + from h5py import File as H5File -from nibabel import Nifti1Header, Nifti1Image -from nibabel.affines import from_matvec from nitransforms.io.base import ( BaseLinearTransformList, - DisplacementsField, - LinearParameters, - TransformIOError, - TransformFileError, ) + class X5Transform: """A string-based structure for X5 linear transforms.""" @@ -28,9 +21,9 @@ def __str__(self): def from_filename(cls, filename): """Read the struct from a X5 file given its path.""" if str(filename).endswith(".h5"): - with H5File(str(filename), 'r') as hdf: + with H5File(str(filename), "r") as hdf: return cls.from_h5obj(hdf) - + @classmethod def from_h5obj(cls, h5obj): """Read the transformations in an X5 file.""" diff --git a/nitransforms/linear.py b/nitransforms/linear.py index d8597dd1..a296f360 100644 --- a/nitransforms/linear.py +++ b/nitransforms/linear.py @@ -7,6 +7,7 @@ # ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ## """Linear transforms.""" + import warnings import numpy as np from pathlib import Path @@ -189,23 +190,31 @@ def _to_hdf5(self, x, x5_root): if self._reference: self.reference._to_hdf5(x5_root.create_group("Reference")) - return #nothing? + return # nothing? def _x5group_affine(self, TransformGroup): """Create group "0" for affine in x5_root/TransformGroup/ according to x5 file format""" aff = TransformGroup.create_group("0") - aff.attrs["Type"] = "affine" #Should have shape {scalar} - aff.attrs["Metadata"] = 'metadata' #This is a draft for metadata. Should have shape {scalar} - aff.create_dataset("Transform", data=[self._matrix]) #Should have shape {3,4} - aff.create_dataset("Inverse", data=[(~self).matrix]) #Should have shape {4,3} + aff.attrs["Type"] = "affine" # Should have shape {scalar} + aff.attrs["Metadata"] = ( + "metadata" # This is a draft for metadata. Should have shape {scalar} + ) + aff.create_dataset("Transform", data=[self._matrix]) # Should have shape {3,4} + aff.create_dataset("Inverse", data=[(~self).matrix]) # Should have shape {4,3} return aff def _x5group_domain(self, x, transform): """Create group "Domain" in x5_root/TransformGroup/0/ according to x5 file format""" coords = transform.create_group("Domain") - coords.attrs["Grid"] = "grid" #How do I interpet this 'grid'? Should have shape {scalar} - coords.create_dataset("Size", data=_as_homogeneous(x, dim=self._matrix.shape[0] - 1).T) #Should have shape {3} - coords.create_dataset("Mapping", data=[self.map(self, x)]) #Should have shape {4,4} + coords.attrs["Grid"] = ( + "grid" # How do I interpet this 'grid'? Should have shape {scalar} + ) + coords.create_dataset( + "Size", data=_as_homogeneous(x, dim=self._matrix.shape[0] - 1).T + ) # Should have shape {3} + coords.create_dataset( + "Mapping", data=[self.map(self, x)] + ) # Should have shape {4,4} return coords def to_filename(self, filename, fmt="X5", moving=None):