From 9c474f1e7c702efc9623defb48f5fe9a6219549f Mon Sep 17 00:00:00 2001 From: Oscar Esteban Date: Thu, 17 Aug 2023 09:23:10 +0200 Subject: [PATCH 1/5] fix: clean up ``scipy.signal.cubic`` deprecation Resolves: #370. --- sdcflows/tests/test_transform.py | 11 ++++++++++- sdcflows/transform.py | 27 +++++++++++++++------------ 2 files changed, 25 insertions(+), 13 deletions(-) diff --git a/sdcflows/tests/test_transform.py b/sdcflows/tests/test_transform.py index 865d3c4cad..72851876be 100644 --- a/sdcflows/tests/test_transform.py +++ b/sdcflows/tests/test_transform.py @@ -348,7 +348,16 @@ def test_grid_bspline_weights(): nb.Nifti1Image(np.zeros(target_shape), target_aff), nb.Nifti1Image(np.zeros(ctrl_shape), ctrl_aff), ).tocsr() - assert weights.shape == (64, 1000) + assert weights.shape == (np.prod(np.array(ctrl_shape) + 2), np.prod(target_shape)) + + # Calculate the legacy mask to index the weights below + legacy_mask = np.pad( + np.ones(ctrl_shape, dtype=bool), + ((1, 1), (1, 1), (1, 1)), + ).reshape(-1) + + # Drop scipy's padding + weights = weights[legacy_mask] # Empirically determined numbers intended to indicate that something # significant has changed. If it turns out we've been doing this wrong, # these numbers will probably change. diff --git a/sdcflows/transform.py b/sdcflows/transform.py index 68949a767c..3df8c4a91f 100644 --- a/sdcflows/transform.py +++ b/sdcflows/transform.py @@ -56,8 +56,8 @@ import numpy as np from warnings import warn from scipy import ndimage as ndi -from scipy.signal import cubic -from scipy.sparse import vstack as sparse_vstack, kron, lil_array +from scipy.interpolate import BSpline +from scipy.sparse import vstack as sparse_vstack, kron import nibabel as nb import nitransforms as nt @@ -325,7 +325,12 @@ def fit( for level in coeffs: wmat = grid_bspline_weights(target_reference, level) weights.append(wmat) - coeffs_data.append(level.get_fdata(dtype="float32").reshape(-1)) + + # Coefficients must be zero-padded because Scipy's BSpline pads the knots. + coeffs_data.append(np.pad( + level.get_fdata(dtype="float32"), + ((1, 1), (1, 1), (1, 1)), + ).reshape(-1)) # Reconstruct the fieldmap (in Hz) from coefficients fmap = np.zeros(projected_reference.shape[:3], dtype="float32") @@ -732,18 +737,16 @@ def grid_bspline_weights(target_nii, ctrl_nii, dtype="float32"): coords[axis] = np.arange(sample_shape[axis], dtype=dtype) # Calculate the index component of samples w.r.t. B-Spline knots along current axis + # Size of locations is L locs = nb.affines.apply_affine(target_to_grid, coords.T)[:, axis] - knots = np.arange(knots_shape[axis], dtype=dtype) - - distance = np.abs(locs[np.newaxis, ...] - knots[..., np.newaxis]) - within_support = distance < 2.0 - d_vals, d_idxs = np.unique(distance[within_support], return_inverse=True) - bs_w = cubic(d_vals) - colloc_ax = lil_array((knots_shape[axis], sample_shape[axis]), dtype=dtype) - colloc_ax[within_support] = bs_w[d_idxs] + # Size of knots is K + 6 so that all locations are fully covered by basis + knots = np.arange(-3, knots_shape[axis] + 3, dtype=dtype) - wd.append(colloc_ax) + # Our original design matrices had size (K, L) + # However, BSpline.design_matrix() generates a size of (L, K + 2), + # hence the trasposition (and zero-padding of 1 at every face when using these) + wd.append(BSpline.design_matrix(locs, knots, 3).T) # Calculate the tensor product of the three design matrices return kron(kron(wd[0], wd[1]), wd[2]).astype(dtype) From 25743844e6b8ec5e30b4066217927ff1852baef5 Mon Sep 17 00:00:00 2001 From: Oscar Esteban Date: Mon, 21 Aug 2023 11:33:58 +0200 Subject: [PATCH 2/5] Apply suggestions from code review Co-authored-by: Chris Markiewicz --- sdcflows/transform.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/sdcflows/transform.py b/sdcflows/transform.py index 3df8c4a91f..be5ffa99e5 100644 --- a/sdcflows/transform.py +++ b/sdcflows/transform.py @@ -745,7 +745,7 @@ def grid_bspline_weights(target_nii, ctrl_nii, dtype="float32"): # Our original design matrices had size (K, L) # However, BSpline.design_matrix() generates a size of (L, K + 2), - # hence the trasposition (and zero-padding of 1 at every face when using these) + # hence the transposition (and zero-padding of 1 at every face when using these) wd.append(BSpline.design_matrix(locs, knots, 3).T) # Calculate the tensor product of the three design matrices From b6a735bb3bb1b90d6b128f2df4d6aa2c43bc8d6a Mon Sep 17 00:00:00 2001 From: Chris Markiewicz Date: Fri, 22 Sep 2023 16:53:26 -0400 Subject: [PATCH 3/5] RF: Sample spline basis less densely Use a construction that will potentially allow us to calculate derivatives more easily. --- sdcflows/tests/test_transform.py | 11 +---------- sdcflows/transform.py | 26 ++++++++++++++------------ 2 files changed, 15 insertions(+), 22 deletions(-) diff --git a/sdcflows/tests/test_transform.py b/sdcflows/tests/test_transform.py index 72851876be..865d3c4cad 100644 --- a/sdcflows/tests/test_transform.py +++ b/sdcflows/tests/test_transform.py @@ -348,16 +348,7 @@ def test_grid_bspline_weights(): nb.Nifti1Image(np.zeros(target_shape), target_aff), nb.Nifti1Image(np.zeros(ctrl_shape), ctrl_aff), ).tocsr() - assert weights.shape == (np.prod(np.array(ctrl_shape) + 2), np.prod(target_shape)) - - # Calculate the legacy mask to index the weights below - legacy_mask = np.pad( - np.ones(ctrl_shape, dtype=bool), - ((1, 1), (1, 1), (1, 1)), - ).reshape(-1) - - # Drop scipy's padding - weights = weights[legacy_mask] + assert weights.shape == (64, 1000) # Empirically determined numbers intended to indicate that something # significant has changed. If it turns out we've been doing this wrong, # these numbers will probably change. diff --git a/sdcflows/transform.py b/sdcflows/transform.py index be5ffa99e5..1ae71fd564 100644 --- a/sdcflows/transform.py +++ b/sdcflows/transform.py @@ -57,7 +57,7 @@ from warnings import warn from scipy import ndimage as ndi from scipy.interpolate import BSpline -from scipy.sparse import vstack as sparse_vstack, kron +from scipy.sparse import vstack as sparse_vstack, kron, lil_array import nibabel as nb import nitransforms as nt @@ -309,7 +309,6 @@ def fit( atol=1e-3, ) - weights = [] if approx: from sdcflows.utils.tools import deoblique_and_zooms @@ -321,16 +320,12 @@ def fit( ) # Generate tensor-product B-Spline weights + weights = [] coeffs_data = [] for level in coeffs: wmat = grid_bspline_weights(target_reference, level) weights.append(wmat) - - # Coefficients must be zero-padded because Scipy's BSpline pads the knots. - coeffs_data.append(np.pad( - level.get_fdata(dtype="float32"), - ((1, 1), (1, 1), (1, 1)), - ).reshape(-1)) + coeffs_data.append(level.get_fdata(dtype="float32").reshape(-1)) # Reconstruct the fieldmap (in Hz) from coefficients fmap = np.zeros(projected_reference.shape[:3], dtype="float32") @@ -743,10 +738,17 @@ def grid_bspline_weights(target_nii, ctrl_nii, dtype="float32"): # Size of knots is K + 6 so that all locations are fully covered by basis knots = np.arange(-3, knots_shape[axis] + 3, dtype=dtype) - # Our original design matrices had size (K, L) - # However, BSpline.design_matrix() generates a size of (L, K + 2), - # hence the transposition (and zero-padding of 1 at every face when using these) - wd.append(BSpline.design_matrix(locs, knots, 3).T) + bspl = BSpline(knots, np.eye(len(knots) - 3 - 1), 3) + + # Construct a sparse design matrix (L, K) + distance = np.abs(locs[..., np.newaxis] - knots[np.newaxis, 3:-3]) + within_support = distance < 2.0 + + colloc_ax = lil_array(distance.shape, dtype=dtype) + colloc_ax[within_support] = bspl(locs)[:, 1:-1][within_support] + + # Transpose to (K, L) and convert to CSR for efficient multiplication + wd.append(colloc_ax.T.tocsr()) # Calculate the tensor product of the three design matrices return kron(kron(wd[0], wd[1]), wd[2]).astype(dtype) From 58ecbe794628493aa671840dc1bc16e3991682c7 Mon Sep 17 00:00:00 2001 From: Chris Markiewicz Date: Sun, 24 Sep 2023 10:58:49 -0400 Subject: [PATCH 4/5] RF: Simplify spline fit/apply by keeping in (L, K) shape --- sdcflows/interfaces/bspline.py | 14 +++++++------- sdcflows/tests/test_transform.py | 6 +++--- sdcflows/transform.py | 27 ++++++++++++--------------- 3 files changed, 22 insertions(+), 25 deletions(-) diff --git a/sdcflows/interfaces/bspline.py b/sdcflows/interfaces/bspline.py index 767e055927..9219ddf7fd 100644 --- a/sdcflows/interfaces/bspline.py +++ b/sdcflows/interfaces/bspline.py @@ -130,7 +130,7 @@ class BSplineApprox(SimpleInterface): def _run_interface(self, runtime): from sklearn import linear_model as lm - from scipy.sparse import vstack as sparse_vstack + from scipy.sparse import hstack as sparse_hstack # Output name baseline out_name = fname_presuffix( @@ -197,9 +197,9 @@ def _run_interface(self, runtime): data -= center # Calculate collocation matrix from (possibly resized) image and knot grids - colmat = sparse_vstack( - grid_bspline_weights(fmapnii, grid) for grid in bs_grids - ).T.tocsr() + colmat = sparse_hstack( + [grid_bspline_weights(fmapnii, grid) for grid in bs_grids] + ).tocsr() bs_grids_str = ["x".join(str(s) for s in grid.shape) for grid in bs_grids] bs_grids_str[-1] = f"and {bs_grids_str[-1]}" @@ -254,9 +254,9 @@ def _run_interface(self, runtime): mask = np.asanyarray(masknii.dataobj) > 1e-4 else: mask = np.ones_like(fmapnii.dataobj, dtype=bool) - colmat = sparse_vstack( - grid_bspline_weights(fmapnii, grid) for grid in bs_grids - ).T.tocsr() + colmat = sparse_hstack( + [grid_bspline_weights(fmapnii, grid) for grid in bs_grids] + ).tocsr() regressors = colmat[mask.reshape(-1), :] interp_data = np.zeros_like(data) diff --git a/sdcflows/tests/test_transform.py b/sdcflows/tests/test_transform.py index 865d3c4cad..eb8bd09bb2 100644 --- a/sdcflows/tests/test_transform.py +++ b/sdcflows/tests/test_transform.py @@ -348,11 +348,11 @@ def test_grid_bspline_weights(): nb.Nifti1Image(np.zeros(target_shape), target_aff), nb.Nifti1Image(np.zeros(ctrl_shape), ctrl_aff), ).tocsr() - assert weights.shape == (64, 1000) + assert weights.shape == (1000, 64) # Empirically determined numbers intended to indicate that something # significant has changed. If it turns out we've been doing this wrong, # these numbers will probably change. assert np.isclose(weights[0, 0], 0.00089725334) assert np.isclose(weights[-1, -1], 0.18919244) - assert np.isclose(weights.sum(axis=1).max(), 129.3907) - assert np.isclose(weights.sum(axis=1).min(), 0.0052327816) + assert np.isclose(weights.sum(axis=0).max(), 129.3907) + assert np.isclose(weights.sum(axis=0).min(), 0.0052327816) diff --git a/sdcflows/transform.py b/sdcflows/transform.py index 1ae71fd564..2be476ae21 100644 --- a/sdcflows/transform.py +++ b/sdcflows/transform.py @@ -57,7 +57,7 @@ from warnings import warn from scipy import ndimage as ndi from scipy.interpolate import BSpline -from scipy.sparse import vstack as sparse_vstack, kron, lil_array +from scipy.sparse import hstack as sparse_hstack, kron, lil_array import nibabel as nb import nitransforms as nt @@ -320,18 +320,15 @@ def fit( ) # Generate tensor-product B-Spline weights - weights = [] - coeffs_data = [] - for level in coeffs: - wmat = grid_bspline_weights(target_reference, level) - weights.append(wmat) - coeffs_data.append(level.get_fdata(dtype="float32").reshape(-1)) + colmat = sparse_hstack( + [grid_bspline_weights(target_reference, level) for level in coeffs] + ).tocsr() + coefficients = np.hstack( + [level.get_fdata(dtype="float32").reshape(-1) for level in coeffs] + ) # Reconstruct the fieldmap (in Hz) from coefficients - fmap = np.zeros(projected_reference.shape[:3], dtype="float32") - fmap = (np.squeeze(np.hstack(coeffs_data).T) @ sparse_vstack(weights)).reshape( - fmap.shape - ) + fmap = np.reshape(colmat @ coefficients, projected_reference.shape[:3]) # Generate a NIfTI object hdr = target_reference.header.copy() @@ -703,7 +700,7 @@ def grid_bspline_weights(target_nii, ctrl_nii, dtype="float32"): Returns ------- - weights : :obj:`numpy.ndarray` (:math:`K \times N`) + weights : :obj:`numpy.ndarray` (:math:`N \times K`) A sparse matrix of interpolating weights :math:`\Psi^3(\mathbf{k}, \mathbf{s})` for the *N* voxels of the target EPI, for each of the total *K* knots. This sparse matrix can be directly used as design matrix for the fitting @@ -747,11 +744,11 @@ def grid_bspline_weights(target_nii, ctrl_nii, dtype="float32"): colloc_ax = lil_array(distance.shape, dtype=dtype) colloc_ax[within_support] = bspl(locs)[:, 1:-1][within_support] - # Transpose to (K, L) and convert to CSR for efficient multiplication - wd.append(colloc_ax.T.tocsr()) + # Convert to CSR for efficient multiplication + wd.append(colloc_ax.tocsr()) # Calculate the tensor product of the three design matrices - return kron(kron(wd[0], wd[1]), wd[2]).astype(dtype) + return kron(kron(wd[0], wd[1]), wd[2]) def _move_coeff(in_coeff, fmap_ref, transform, fmap_target=None): From 1d715a9f1b933c433328b37e76fb55369a69c0af Mon Sep 17 00:00:00 2001 From: Chris Markiewicz Date: Wed, 27 Sep 2023 10:44:30 -0400 Subject: [PATCH 5/5] Update sdcflows/transform.py --- sdcflows/transform.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/sdcflows/transform.py b/sdcflows/transform.py index 2be476ae21..77ac1687e1 100644 --- a/sdcflows/transform.py +++ b/sdcflows/transform.py @@ -321,7 +321,7 @@ def fit( # Generate tensor-product B-Spline weights colmat = sparse_hstack( - [grid_bspline_weights(target_reference, level) for level in coeffs] + [grid_bspline_weights(projected_reference, level) for level in coeffs] ).tocsr() coefficients = np.hstack( [level.get_fdata(dtype="float32").reshape(-1) for level in coeffs]