Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

share a 2D nufft across Stacked Operators #87

Merged
merged 6 commits into from
Apr 22, 2024
Merged

share a 2D nufft across Stacked Operators #87

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
6 commits
Select commit Hold shift + click to select a range
3ed0924
feat: allow reuse of 2D operator in stacked nufft.
paquiteau Mar 18, 2024
b044ff0
use shared operator for cufinufft as well.
paquiteau Mar 19, 2024
cbcc266
add test for stacked reuse
paquiteau Apr 4, 2024
04c9a55
Merge branch 'master' into stacked-shared
paquiteau Apr 5, 2024
630bfbd
fix: E501
paquiteau Apr 5, 2024
1ac1e26
fix: ruff.
paquiteau Apr 22, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
133 changes: 91 additions & 42 deletions src/mrinufft/operators/stacked.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -6,7 +6,12 @@
import scipy as sp

from mrinufft._utils import proper_trajectory, power_method, get_array_module, auto_cast
from mrinufft.operators.base import FourierOperatorBase, check_backend, get_operator
from mrinufft.operators.base import (
FourierOperatorBase,
check_backend,
get_operator,
with_numpy_cupy,
)
from mrinufft.operators.interfaces.utils import (
is_cuda_array,
is_host_array,
Expand Down Expand Up @@ -36,8 +41,13 @@ class MRIStackedNUFFT(FourierOperatorBase):
Shape of the image.
z_index: array-like
Cartesian z index of masked plan.
backend: str
backend: str or FourierOperatorBase
Backend to use.
If str, a NUFFT operator is initialized with str being a registered backend.
If FourierOperatorBase, operator is checked for compatibility and used as is
notably one should have:
``n_coils = self.n_coils*len(z_index), squeeze_dims=True, smaps=None``

smaps: array-like
Sensitivity maps.
n_coils: int
Expand Down Expand Up @@ -73,22 +83,48 @@ def __init__(
**kwargs,
):
super().__init__()
samples2d, z_index_ = self._init_samples(samples, z_index, shape)
self.shape = shape
self.samples = samples2d.reshape(-1, 2)
self.z_index = z_index_
self.n_coils = n_coils
self.n_batchs = n_batchs
self.squeeze_dims = squeeze_dims
self.smaps = smaps
self.operator = get_operator(backend)(
self.samples,
shape[:-1],
n_coils=self.n_coils * len(self.z_index),
smaps=None,
squeeze_dims=True,
**kwargs,
)
if isinstance(backend, str):
samples2d, z_index_ = self._init_samples(samples, z_index, shape)
self.samples = samples2d.reshape(-1, 2)
self.z_index = z_index_
self.operator = get_operator(backend)(
self.samples,
shape[:-1],
n_coils=self.n_coils * len(self.z_index),
smaps=None,
squeeze_dims=True,
**kwargs,
)
elif isinstance(backend, FourierOperatorBase):
# get all the interesting values from the operator
if backend.shape != shape[:-1]:
raise ValueError("Backend operator should have compatible shape")

samples2d, z_index_ = self._init_samples(backend.samples, z_index, shape)
self.samples = samples2d.reshape(-1, 2)
self.z_index = z_index_

if backend.n_coils != self.n_coils * (len(z_index_)):
raise ValueError(
"The backend operator should have ``n_coils * len(z_index)``"
" specified for its coil dimension."
)
if backend.uses_sense:
raise ValueError("Backend operator should not uses smaps.")
if not backend.squeeze_dims:
raise ValueError("Backend operator should have ``squeeze_dims=True``")
self.operator = backend

else:
raise ValueError(
"backend should either be a 2D nufft operator,"
" or a str specifying which nufft library to use."
)

@staticmethod
def _init_samples(samples, z_index, shape):
Expand Down Expand Up @@ -145,6 +181,7 @@ def _ifftz(data):
sp.fft.ifft(sp.fft.ifftshift(data, axes=-1), axis=-1, norm="ortho"), axes=-1
) / np.sqrt(2)

@with_numpy_cupy
def op(self, data, ksp=None):
"""Forward operator."""
if self.uses_sense:
Expand Down Expand Up @@ -189,6 +226,7 @@ def _op_calibless(self, data, ksp=None):
ksp = ksp.reshape((B, C, NZ * NS))
return ksp

@with_numpy_cupy
def adj_op(self, coeffs, img=None):
"""Adjoint operator."""
if self.uses_sense:
Expand Down Expand Up @@ -285,6 +323,7 @@ def __init__(
squeeze_dims=False,
smaps_cached=False,
density=False,
backend="cufinufft",
**kwargs,
):
if not (CUPY_AVAILABLE and check_backend("cufinufft")):
Expand All @@ -293,25 +332,49 @@ def __init__(
if (n_batchs * n_coils) % n_trans != 0:
raise ValueError("n_batchs * n_coils should be a multiple of n_transf")

samples2d, z_index_ = self._init_samples(samples, z_index, shape)
self.shape = shape
self.samples = samples2d.reshape(-1, 2)
self.z_index = z_index_
self.n_coils = n_coils
self.n_batchs = n_batchs
self.n_trans = n_trans
self.squeeze_dims = squeeze_dims
if isinstance(backend, str):
samples2d, z_index_ = self._init_samples(samples, z_index, shape)
self.samples = samples2d.reshape(-1, 2)
self.z_index = z_index_
self.operator = get_operator(backend)(
self.samples,
shape[:-1],
n_coils=self.n_trans * len(self.z_index),
n_trans=len(self.z_index),
smaps=None,
squeeze_dims=True,
**kwargs,
)
elif isinstance(backend, FourierOperatorBase):
# get all the interesting values from the operator
if backend.shape != shape[:-1]:
raise ValueError("Backend operator should have compatible shape")

samples2d, z_index_ = self._init_samples(backend.samples, z_index, shape)
self.samples = samples2d.reshape(-1, 2)
self.z_index = z_index_

if backend.n_coils != self.n_trans * len(z_index_):
raise ValueError(
"The backend operator should have ``n_coils * len(z_index)``"
" specified for its coil dimension."
)
if backend.uses_sense:
raise ValueError("Backend operator should not uses smaps.")
if not backend.squeeze_dims:
raise ValueError("Backend operator should have ``squeeze_dims=True``")
self.operator = backend
else:
raise ValueError(
"backend should either be a 2D nufft operator,"
" or a str specifying which nufft library to use."
)

self.operator = get_operator("cufinufft")(
self.samples,
shape[:-1],
n_coils=n_trans * len(self.z_index),
n_trans=len(self.z_index),
smaps=None,
squeeze_dims=True,
density=density,
**kwargs,
)
# Smaps support
self.smaps = smaps
self.smaps_cached = False
Expand All @@ -320,7 +383,6 @@ def __init__(
raise ValueError(
"Smaps should be either a C-ordered ndarray, " "or a GPUArray."
)
self.smaps_cached = False
if smaps_cached:
warnings.warn(
f"{sizeof_fmt(smaps.size * np.dtype(self.cpx_dtype).itemsize)}"
Expand Down Expand Up @@ -367,12 +429,10 @@ def _ifftz(data):
axes=-1,
)

@with_numpy_cupy
def op(self, data, ksp=None):
"""Forward operator."""
# Dispatch to special case.
xp = get_array_module(data)
if xp.__name__ == "torch" and data.is_cpu:
data = data.numpy()
data = auto_cast(data, self.cpx_dtype)

if self.uses_sense and is_cuda_array(data):
Expand All @@ -385,10 +445,6 @@ def op(self, data, ksp=None):
op_func = self._op_calibless_host
ret = op_func(data, ksp)

if xp.__name__ == "torch" and is_cuda_array(ret):
ret = xp.as_tensor(ret, device=data.device)
elif xp.__name__ == "torch":
ret = xp.from_numpy(ret)
return self._safe_squeeze(ret)

def _op_sense_host(self, data, ksp=None):
Expand Down Expand Up @@ -527,12 +583,10 @@ def _op_calibless_device(self, data, ksp=None):
ksp = ksp.reshape((B, C, NZ * NS))
return ksp

@with_numpy_cupy
def adj_op(self, coeffs, img=None):
"""Adjoint operator."""
# Dispatch to special case.
xp = get_array_module(coeffs)
if xp.__name__ == "torch" and coeffs.is_cpu:
coeffs = coeffs.numpy()
coeffs = auto_cast(coeffs, self.cpx_dtype)

if self.uses_sense and is_cuda_array(coeffs):
Expand All @@ -546,11 +600,6 @@ def adj_op(self, coeffs, img=None):

ret = adj_op_func(coeffs, img)

if xp.__name__ == "torch" and is_cuda_array(ret):
ret = xp.as_tensor(ret, device=coeffs.device)
elif xp.__name__ == "torch":
ret = xp.from_numpy(ret)

return self._safe_squeeze(ret)

def _adj_op_sense_host(self, coeffs, img_d=None):
Expand Down
16 changes: 16 additions & 0 deletions tests/operators/test_stacked.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -134,3 +134,19 @@ def test_stacked2traj3d():

npt.assert_allclose(traj2d, traj2d)
npt.assert_allclose(z_index, z_index2)


def test_stack_reuse(operator, stacked_op):
"""Test the reuse of the stacked operator."""
nufft_2d = stacked_op.operator

reuse_op = MRIStackedNUFFT(
backend=nufft_2d,
shape=stacked_op.shape,
samples=stacked_op.samples,
z_index=stacked_op.z_index,
n_coils=stacked_op.n_coils,
n_batchs=stacked_op.n_batchs,
smaps=stacked_op.smaps,
)
assert reuse_op.operator is nufft_2d
Loading