From db1b250dbbfa928785a2d469f246369e130f3624 Mon Sep 17 00:00:00 2001 From: Oscar Esteban Date: Thu, 16 Nov 2023 17:49:46 +0100 Subject: [PATCH 1/7] fix: postpone coordinate mapping on linear array transforms Resolves: #173. --- nitransforms/linear.py | 68 +++++++++++++++++++++++++----------------- 1 file changed, 41 insertions(+), 27 deletions(-) diff --git a/nitransforms/linear.py b/nitransforms/linear.py index 9c430d3b..0709d50b 100644 --- a/nitransforms/linear.py +++ b/nitransforms/linear.py @@ -436,6 +436,7 @@ def apply( The data imaged after resampling to reference space. """ + if reference is not None and isinstance(reference, (str, Path)): reference = _nbload(str(reference)) @@ -446,40 +447,53 @@ def apply( if isinstance(spatialimage, (str, Path)): spatialimage = _nbload(str(spatialimage)) - data = np.squeeze(np.asanyarray(spatialimage.dataobj)) - output_dtype = output_dtype or data.dtype + # Avoid opening the data array just yet + input_dtype = spatialimage.header.get_data_dtype() + output_dtype = output_dtype or input_dtype - ycoords = self.map(_ref.ndcoords.T) - targets = ImageGrid(spatialimage).index( # data should be an image - _as_homogeneous(np.vstack(ycoords), dim=_ref.ndim) - ) + # Prepare physical coordinates of input (grid, points) + xcoords = _ref.ndcoords.astype("f4") - if data.ndim == 4: - if len(self) != data.shape[-1]: + # Invert target's (moving) affine once + ras2vox = ~Affine(spatialimage.affine) + + if spatialimage.ndim == 4: + if len(self) != spatialimage.shape[-1]: raise ValueError( "Attempting to apply %d transforms on a file with " - "%d timepoints" % (len(self), data.shape[-1]) + "%d timepoints" % (len(self), spatialimage.shape[-1]) ) - targets = targets.reshape((len(self), -1, targets.shape[-1])) - resampled = np.stack( - [ - ndi.map_coordinates( - data[..., t], - targets[t, ..., : _ref.ndim].T, - output=output_dtype, - order=order, - mode=mode, - cval=cval, - prefilter=prefilter, - ) - for t in range(data.shape[-1]) - ], - axis=0, + + # Order F ensures individual volumes are contiguous in memory + # Also matches NIfTI, making final save more efficient + resampled = np.zeros( + (xcoords.T.shape[0], ) + spatialimage.shape[-1:], dtype=output_dtype, order="F" ) - elif data.ndim in (2, 3): + + for t in range(spatialimage.shape[-1]): + # Map the input coordinates on to timepoint t of the target (moving) + ycoords = Affine(self.matrix[t]).map(xcoords.T)[..., : _ref.ndim] + + # Calculate corresponding voxel coordinates + yvoxels = ras2vox.map(ycoords)[..., : _ref.ndim] + + # Interpolate + resampled[..., t] = ndi.map_coordinates( + spatialimage.dataobj[..., t].astype(input_dtype, copy=False), + yvoxels.T, + output=output_dtype, + order=order, + mode=mode, + cval=cval, + prefilter=prefilter, + ) + elif spatialimage.ndim in (2, 3): + ycoords = self.map(xcoords.T)[..., : _ref.ndim] + yvoxels = ras2vox.map(ycoords)[..., : _ref.ndim] + resampled = ndi.map_coordinates( - data, - targets[..., : _ref.ndim].T, + spatialimage.dataobj.astype(input_dtype, copy=False), + yvoxels.T, output=output_dtype, order=order, mode=mode, From c0e7f2758a0a013d7f5377d25059456365465232 Mon Sep 17 00:00:00 2001 From: Oscar Esteban Date: Thu, 16 Nov 2023 21:50:39 +0100 Subject: [PATCH 2/7] Update nitransforms/linear.py --- nitransforms/linear.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/nitransforms/linear.py b/nitransforms/linear.py index 0709d50b..3995bd28 100644 --- a/nitransforms/linear.py +++ b/nitransforms/linear.py @@ -448,7 +448,7 @@ def apply( spatialimage = _nbload(str(spatialimage)) # Avoid opening the data array just yet - input_dtype = spatialimage.header.get_data_dtype() + input_dtype = nb.arrayproxy.get_obj_dtype(spatialimage.dataobj) output_dtype = output_dtype or input_dtype # Prepare physical coordinates of input (grid, points) From 2c36d088e7d8dae658a064211fdd394caf214ed9 Mon Sep 17 00:00:00 2001 From: Oscar Esteban Date: Thu, 16 Nov 2023 21:53:01 +0100 Subject: [PATCH 3/7] fix: missing import --- nitransforms/linear.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/nitransforms/linear.py b/nitransforms/linear.py index 3995bd28..dfeca36d 100644 --- a/nitransforms/linear.py +++ b/nitransforms/linear.py @@ -14,6 +14,7 @@ from nibabel.loadsave import load as _nbload from nibabel.affines import from_matvec +from nibabel.arrayproxy import get_obj_dtype from nitransforms.base import ( ImageGrid, @@ -448,7 +449,7 @@ def apply( spatialimage = _nbload(str(spatialimage)) # Avoid opening the data array just yet - input_dtype = nb.arrayproxy.get_obj_dtype(spatialimage.dataobj) + input_dtype = get_obj_dtype(spatialimage.dataobj) output_dtype = output_dtype or input_dtype # Prepare physical coordinates of input (grid, points) From 9772710419968ef362e1fb486b29c31ae00cbd1f Mon Sep 17 00:00:00 2001 From: Oscar Esteban Date: Thu, 16 Nov 2023 23:15:44 +0100 Subject: [PATCH 4/7] fix: shape and order of resampled array --- nitransforms/linear.py | 68 +++++++++++++++++++++--------------------- 1 file changed, 34 insertions(+), 34 deletions(-) diff --git a/nitransforms/linear.py b/nitransforms/linear.py index dfeca36d..8511124c 100644 --- a/nitransforms/linear.py +++ b/nitransforms/linear.py @@ -317,6 +317,11 @@ def __init__(self, transforms, reference=None): ) self._inverse = np.linalg.inv(self._matrix) + def __iter__(self): + """Enable iterating over the series of transforms.""" + for _m in self.matrix: + yield Affine(_m, reference=self._reference) + def __getitem__(self, i): """Enable indexed access to the series of matrices.""" return Affine(self.matrix[i, ...], reference=self._reference) @@ -458,42 +463,37 @@ def apply( # Invert target's (moving) affine once ras2vox = ~Affine(spatialimage.affine) - if spatialimage.ndim == 4: - if len(self) != spatialimage.shape[-1]: - raise ValueError( - "Attempting to apply %d transforms on a file with " - "%d timepoints" % (len(self), spatialimage.shape[-1]) - ) - - # Order F ensures individual volumes are contiguous in memory - # Also matches NIfTI, making final save more efficient - resampled = np.zeros( - (xcoords.T.shape[0], ) + spatialimage.shape[-1:], dtype=output_dtype, order="F" + if spatialimage.ndim == 4 and (len(self) != spatialimage.shape[-1]): + raise ValueError( + "Attempting to apply %d transforms on a file with " + "%d timepoints" % (len(self), spatialimage.shape[-1]) ) - for t in range(spatialimage.shape[-1]): - # Map the input coordinates on to timepoint t of the target (moving) - ycoords = Affine(self.matrix[t]).map(xcoords.T)[..., : _ref.ndim] - - # Calculate corresponding voxel coordinates - yvoxels = ras2vox.map(ycoords)[..., : _ref.ndim] - - # Interpolate - resampled[..., t] = ndi.map_coordinates( - spatialimage.dataobj[..., t].astype(input_dtype, copy=False), - yvoxels.T, - output=output_dtype, - order=order, - mode=mode, - cval=cval, - prefilter=prefilter, - ) - elif spatialimage.ndim in (2, 3): - ycoords = self.map(xcoords.T)[..., : _ref.ndim] + # Order F ensures individual volumes are contiguous in memory + # Also matches NIfTI, making final save more efficient + resampled = np.zeros( + (xcoords.T.shape[0], len(self)), dtype=output_dtype, order="F" + ) + + dataobj = ( + np.asanyarray(spatialimage.dataobj, dtype=input_dtype) + if spatialimage.ndim in (2, 3) + else None + ) + + for t, xfm_t in enumerate(self): + # Map the input coordinates on to timepoint t of the target (moving) + ycoords = xfm_t.map(xcoords.T)[..., : _ref.ndim] + + # Calculate corresponding voxel coordinates yvoxels = ras2vox.map(ycoords)[..., : _ref.ndim] - resampled = ndi.map_coordinates( - spatialimage.dataobj.astype(input_dtype, copy=False), + # Interpolate + resampled[..., t] = ndi.map_coordinates( + ( + dataobj if dataobj is not None + else np.asanyarray(spatialimage.dataobj[..., t], dtype=input_dtype) + ), yvoxels.T, output=output_dtype, order=order, @@ -503,9 +503,9 @@ def apply( ) if isinstance(_ref, ImageGrid): # If reference is grid, reshape - newdata = resampled.reshape((len(self), *_ref.shape)) + newdata = resampled.reshape(_ref.shape + (len(self), )) moved = spatialimage.__class__( - np.moveaxis(newdata, 0, -1), _ref.affine, spatialimage.header + newdata, _ref.affine, spatialimage.header ) moved.header.set_data_dtype(output_dtype) return moved From 40e3c13c66054e66c960dacc8d5da101986d574d Mon Sep 17 00:00:00 2001 From: Oscar Esteban Date: Thu, 16 Nov 2023 23:18:31 +0100 Subject: [PATCH 5/7] enh: one less iterated operation (transpose) --- nitransforms/linear.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/nitransforms/linear.py b/nitransforms/linear.py index 8511124c..12409885 100644 --- a/nitransforms/linear.py +++ b/nitransforms/linear.py @@ -458,7 +458,7 @@ def apply( output_dtype = output_dtype or input_dtype # Prepare physical coordinates of input (grid, points) - xcoords = _ref.ndcoords.astype("f4") + xcoords = _ref.ndcoords.astype("f4").T # Invert target's (moving) affine once ras2vox = ~Affine(spatialimage.affine) @@ -472,7 +472,7 @@ def apply( # Order F ensures individual volumes are contiguous in memory # Also matches NIfTI, making final save more efficient resampled = np.zeros( - (xcoords.T.shape[0], len(self)), dtype=output_dtype, order="F" + (xcoords.shape[0], len(self)), dtype=output_dtype, order="F" ) dataobj = ( @@ -483,7 +483,7 @@ def apply( for t, xfm_t in enumerate(self): # Map the input coordinates on to timepoint t of the target (moving) - ycoords = xfm_t.map(xcoords.T)[..., : _ref.ndim] + ycoords = xfm_t.map(xcoords)[..., : _ref.ndim] # Calculate corresponding voxel coordinates yvoxels = ras2vox.map(ycoords)[..., : _ref.ndim] From a6fe3ebd8e0bc5f1010375e8532d38fb99943086 Mon Sep 17 00:00:00 2001 From: Oscar Esteban Date: Thu, 16 Nov 2023 23:26:29 +0100 Subject: [PATCH 6/7] sty: run black --- nitransforms/linear.py | 14 ++++++-------- 1 file changed, 6 insertions(+), 8 deletions(-) diff --git a/nitransforms/linear.py b/nitransforms/linear.py index 12409885..fcee7432 100644 --- a/nitransforms/linear.py +++ b/nitransforms/linear.py @@ -217,14 +217,13 @@ def from_filename(cls, filename, fmt=None, reference=None, moving=None): is_array = cls != Affine errors = [] for potential_fmt in fmtlist: - if (potential_fmt == "itk" and Path(filename).suffix == ".mat"): + if potential_fmt == "itk" and Path(filename).suffix == ".mat": is_array = False cls = Affine try: struct = get_linear_factory( - potential_fmt, - is_array=is_array + potential_fmt, is_array=is_array ).from_filename(filename) except (TransformFileError, FileNotFoundError) as err: errors.append((potential_fmt, err)) @@ -491,7 +490,8 @@ def apply( # Interpolate resampled[..., t] = ndi.map_coordinates( ( - dataobj if dataobj is not None + dataobj + if dataobj is not None else np.asanyarray(spatialimage.dataobj[..., t], dtype=input_dtype) ), yvoxels.T, @@ -503,10 +503,8 @@ def apply( ) if isinstance(_ref, ImageGrid): # If reference is grid, reshape - newdata = resampled.reshape(_ref.shape + (len(self), )) - moved = spatialimage.__class__( - newdata, _ref.affine, spatialimage.header - ) + newdata = resampled.reshape(_ref.shape + (len(self),)) + moved = spatialimage.__class__(newdata, _ref.affine, spatialimage.header) moved.header.set_data_dtype(output_dtype) return moved From d148e85e333f5cb6c737879e9d33fd573c8ae9b4 Mon Sep 17 00:00:00 2001 From: Oscar Esteban Date: Thu, 16 Nov 2023 23:55:39 +0100 Subject: [PATCH 7/7] fix: access ``__getitem__`` directly Co-authored-by: Chris Markiewicz --- nitransforms/linear.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/nitransforms/linear.py b/nitransforms/linear.py index fcee7432..eb4a95d7 100644 --- a/nitransforms/linear.py +++ b/nitransforms/linear.py @@ -492,7 +492,7 @@ def apply( ( dataobj if dataobj is not None - else np.asanyarray(spatialimage.dataobj[..., t], dtype=input_dtype) + else spatialimage.dataobj[..., t].astype(input_dtype, copy=False) ), yvoxels.T, output=output_dtype,