From 8ca012cd9d96ebb4488eb291a3401d3584189171 Mon Sep 17 00:00:00 2001 From: chaithyagr Date: Fri, 26 Apr 2024 14:45:12 +0200 Subject: [PATCH 01/16] Added --- .../operators/interfaces/nudft_numpy.py | 40 +++-- tests/test_autodiff.py | 145 ++++++++++++++++++ tests/test_ndft.py | 59 +++---- 3 files changed, 205 insertions(+), 39 deletions(-) create mode 100644 tests/test_autodiff.py diff --git a/src/mrinufft/operators/interfaces/nudft_numpy.py b/src/mrinufft/operators/interfaces/nudft_numpy.py index 31a9fa18a..9d26d02f5 100644 --- a/src/mrinufft/operators/interfaces/nudft_numpy.py +++ b/src/mrinufft/operators/interfaces/nudft_numpy.py @@ -8,6 +8,7 @@ from ..base import FourierOperatorCPU +<<<<<<< Updated upstream def get_fourier_matrix(ktraj, shape): """Get the NDFT Fourier Matrix.""" n = np.prod(shape) @@ -17,47 +18,64 @@ def get_fourier_matrix(ktraj, shape): grid_r = np.reshape(np.meshgrid(*r, indexing="ij"), (ndim, np.prod(shape))) traj_grid = ktraj @ grid_r matrix = np.exp(-2j * np.pi * traj_grid) +======= +def get_fourier_matrix(ktraj, shape, dtype=np.complex64, normalize=False): + """Get the NDFT Fourier Matrix.""" + n = np.prod(shape) + ndim = len(shape) + matrix = np.zeros((len(ktraj), n), dtype=dtype) + r = [np.linspace(-s/2, s/2-1, s) for s in shape] + grid_r = np.reshape(np.meshgrid(*r, indexing="ij"), (ndim, np.prod(shape))) + traj_grid = ktraj @ grid_r + matrix = np.exp(-2j * np.pi * traj_grid, dtype=dtype) + if normalize: + matrix /= (np.sqrt(np.prod(shape)) * np.power(np.sqrt(2), len(shape))) +>>>>>>> Stashed changes return matrix -def implicit_type2_ndft(ktraj, image, shape): +def implicit_type2_ndft(ktraj, image, shape, normalize=False): """Compute the NDFT using the implicit type 2 (image -> kspace) algorithm.""" - r = [np.arange(s) for s in shape] + r = [np.linspace(-s/2, s/2-1, s) for s in shape] grid_r = np.reshape( np.meshgrid(*r, indexing="ij"), (len(shape), np.prod(image.shape)) ) res = np.zeros(len(ktraj), dtype=image.dtype) for j in range(np.prod(image.shape)): res += image[j] * np.exp(-2j * np.pi * ktraj @ grid_r[:, j]) + if normalize: + matrix /= (np.sqrt(np.prod(shape)) * np.power(np.sqrt(2), len(shape))) return res -def implicit_type1_ndft(ktraj, coeffs, shape): +def implicit_type1_ndft(ktraj, coeffs, shape, normalize=False): """Compute the NDFT using the implicit type 1 (kspace -> image) algorithm.""" - r = [np.arange(s) for s in shape] + r = [np.linspace(-s/2, s/2-1, s) for s in shape] grid_r = np.reshape(np.meshgrid(*r, indexing="ij"), (len(shape), np.prod(shape))) res = np.zeros(np.prod(shape), dtype=coeffs.dtype) for i in range(len(ktraj)): res += coeffs[i] * np.exp(2j * np.pi * ktraj[i] @ grid_r) + if normalize: + matrix /= (np.sqrt(np.prod(shape)) * np.power(np.sqrt(2), len(shape))) return res -def get_implicit_matrix(ktraj, shape): +def get_implicit_matrix(ktraj, shape, normalize=False): """Get the NDFT Fourier Matrix as implicit operator. This is more memory efficient than the explicit matrix. """ return sp.sparse.linalg.LinearOperator( (len(ktraj), np.prod(shape)), - matvec=lambda x: implicit_type2_ndft(ktraj, x, shape), - rmatvec=lambda x: implicit_type1_ndft(ktraj, x, shape), + matvec=lambda x: implicit_type2_ndft(ktraj, x, shape, normalize), + rmatvec=lambda x: implicit_type1_ndft(ktraj, x, shape, normalize), ) class RawNDFT: """Implementation of the NUDFT using numpy.""" - def __init__(self, samples, shape, explicit_matrix=True): + def __init__(self, samples, shape, explicit_matrix=True, normalize=False): self.samples = samples self.shape = shape self.n_samples = len(samples) @@ -65,13 +83,13 @@ def __init__(self, samples, shape, explicit_matrix=True): if explicit_matrix: try: self._fourier_matrix = sp.sparse.linalg.aslinearoperator( - get_fourier_matrix(self.samples, self.shape) + get_fourier_matrix(self.samples, self.shape, normalize=normalize) ) except MemoryError: warnings.warn("Not enough memory, using an implicit definition anyway") - self._fourier_matrix = get_implicit_matrix(self.samples, self.shape) + self._fourier_matrix = get_implicit_matrix(self.samples, self.shape, normalize) else: - self._fourier_matrix = get_implicit_matrix(self.samples, self.shape) + self._fourier_matrix = get_implicit_matrix(self.samples, self.shape, normalize) def op(self, coeffs, image): """Compute the forward NUDFT.""" diff --git a/tests/test_autodiff.py b/tests/test_autodiff.py new file mode 100644 index 000000000..01de2f91c --- /dev/null +++ b/tests/test_autodiff.py @@ -0,0 +1,145 @@ +"""Test the autodiff functionnality.""" + +import numpy as np +from mrinufft.operators.interfaces.nudft_numpy import get_fourier_matrix +import pytest +from pytest_cases import parametrize_with_cases, parametrize, fixture +from case_trajectories import CasesTrajectories +from mrinufft.operators import get_operator + + +from helpers import ( + kspace_from_op, + image_from_op, + to_interface, +) + + +TORCH_AVAILABLE = True +try: + import torch + import torch.testing as tt +except ImportError: + TORCH_AVAILABLE = False + + +@fixture(scope="module") +@parametrize(backend=["cufinufft", "finufft"]) +@parametrize(autograd=["data"]) +@parametrize_with_cases( + "kspace_loc, shape", + cases=[ + CasesTrajectories.case_grid2D, + CasesTrajectories.case_nyquist_radial2D, + ], # 2D cases only for reduced memory footprint. +) +def operator(kspace_loc, shape, backend, autograd): + """Create NUFFT operator with autodiff capabilities.""" + kspace_loc = kspace_loc.astype(np.float32) + + nufft = get_operator(backend_name=backend, autograd=autograd)( + samples=kspace_loc, + shape=shape, + smaps=None, + ) + + return nufft + + +@fixture(scope="module") +def ndft_matrix(operator): + """Get the NDFT matrix from the operator.""" + return get_fourier_matrix(operator.samples, operator.shape, normalize=True) + + +@pytest.mark.parametrize("interface", ["torch-gpu", "torch-cpu"]) +@pytest.mark.skipif(not TORCH_AVAILABLE, reason="Pytorch is not installed") +def test_adjoint_and_grad(operator, ndft_matrix, interface): + """Test the adjoint and gradient of the operator.""" + if operator.backend == "finufft" and "gpu" in interface: + pytest.skip("GPU not supported for finufft backend") + ndft_matrix_torch = to_interface(ndft_matrix, interface=interface) + ksp_data = to_interface(kspace_from_op(operator), interface=interface) + img_data = to_interface(image_from_op(operator), interface=interface) + ksp_data.requires_grad = True + + with torch.autograd.set_detect_anomaly(True): + adj_data = operator.adj_op(ksp_data).reshape(img_data.shape) + adj_data_ndft = (ndft_matrix_torch.conj().T @ ksp_data.flatten()).reshape( + adj_data.shape + ) + loss_nufft = torch.mean(torch.abs(adj_data) ** 2) + loss_ndft = torch.mean(torch.abs(adj_data_ndft) ** 2) + + # Check if nufft and ndft are close in the backprop + grad_ndft_kdata = torch.autograd.grad(loss_ndft, ksp_data, retain_graph=True)[0] + grad_nufft_kdata = torch.autograd.grad(loss_nufft, ksp_data, retain_graph=True)[0] + tt.assert_close(grad_ndft_kdata, grad_nufft_kdata, rtol=1, atol=1) + + +@pytest.mark.parametrize("interface", ["torch-gpu", "torch-cpu"]) +@pytest.mark.skipif(not TORCH_AVAILABLE, reason="Pytorch is not installed") +def test_adjoint_and_gradauto(operator, ndft_matrix, interface): + """Test the adjoint and gradient of the operator using autograd gradcheck.""" + if operator.backend == "finufft" and "gpu" in interface: + pytest.skip("GPU not supported for finufft backend") + + ksp_data = to_interface(kspace_from_op(operator), interface=interface) + ksp_data = torch.ones(ksp_data.shape, requires_grad=True, dtype=ksp_data.dtype) + print(ksp_data.shape) + # todo: tighten the tolerance + assert torch.autograd.gradcheck( + operator.adjoint, + ksp_data, + eps=1e-10, + rtol=1, + atol=1, + nondet_tol=1, + raise_exception=True, + ) + + +@pytest.mark.parametrize("interface", ["torch-gpu", "torch-cpu"]) +@pytest.mark.skipif(not TORCH_AVAILABLE, reason="Pytorch is not installed") +def test_forward_and_grad(operator, ndft_matrix, interface): + """Test the adjoint and gradient of the operator.""" + if operator.backend == "finufft" and "gpu" in interface: + pytest.skip("GPU not supported for finufft backend") + + ndft_matrix_torch = to_interface(ndft_matrix, interface=interface) + ksp_data_ref = to_interface(kspace_from_op(operator), interface=interface) + img_data = to_interface(image_from_op(operator), interface=interface) + img_data.requires_grad = True + + with torch.autograd.set_detect_anomaly(True): + ksp_data = operator.op(img_data).reshape(ksp_data_ref.shape) + ksp_data_ndft = (ndft_matrix_torch @ img_data.flatten()).reshape(ksp_data.shape) + loss_nufft = torch.mean(torch.abs(ksp_data - ksp_data_ref) ** 2) + loss_ndft = torch.mean(torch.abs(ksp_data_ndft - ksp_data_ref) ** 2) + + # Check if nufft and ndft are close in the backprop + grad_ndft_kdata = torch.autograd.grad(loss_ndft, img_data, retain_graph=True)[0] + grad_nufft_kdata = torch.autograd.grad(loss_nufft, img_data, retain_graph=True)[0] + assert torch.allclose(grad_ndft_kdata, grad_nufft_kdata, atol=6e-3) + + +@pytest.mark.parametrize("interface", ["torch-gpu", "torch-cpu"]) +@pytest.mark.skipif(not TORCH_AVAILABLE, reason="Pytorch is not installed") +def test_forward_and_gradauto(operator, ndft_matrix, interface): + """Test the forward and gradient of the operator using autograd gradcheck.""" + if operator.backend == "finufft" and "gpu" in interface: + pytest.skip("GPU not supported for finufft backend") + + img_data = to_interface(image_from_op(operator), interface=interface) + img_data = torch.ones(img_data.shape, requires_grad=True, dtype=img_data.dtype) + print(img_data.shape) + # todo: tighten the tolerance + assert torch.autograd.gradcheck( + operator.adjoint, + img_data, + eps=1e-10, + rtol=1, + atol=1, + nondet_tol=1, + raise_exception=True, + ) diff --git a/tests/test_ndft.py b/tests/test_ndft.py index 7ae595a89..5bcb3e8c9 100644 --- a/tests/test_ndft.py +++ b/tests/test_ndft.py @@ -13,32 +13,9 @@ from case_trajectories import CasesTrajectories, case_grid1D from helpers import assert_almost_allclose +from mrinufft import get_operator -@parametrize_with_cases( - "kspace_grid, shape", - cases=[ - case_grid1D, - CasesTrajectories.case_grid2D, - ], # 3D is ignored (too much possibility for the reordering) -) -def test_ndft_grid_matrix(kspace_grid, shape): - """Test that the ndft matrix is a good matrix for doing fft.""" - ndft_matrix = get_fourier_matrix(kspace_grid, shape) - # Create a random image - fft_matrix = [None] * len(shape) - for i in range(len(shape)): - fft_matrix[i] = sp.fft.fft(np.eye(shape[i])) - fft_mat = fft_matrix[0] - if len(shape) == 2: - fft_mat = fft_matrix[0].flatten()[:, None] @ fft_matrix[1].flatten()[None, :] - fft_mat = ( - fft_mat.reshape(shape * 2) - .transpose(2, 0, 1, 3) - .reshape((np.prod(shape),) * 2) - ) - assert np.allclose(ndft_matrix, fft_mat) - @parametrize_with_cases( "kspace, shape", @@ -56,7 +33,7 @@ def test_ndft_implicit2(kspace, shape): linop_coef = implicit_type2_ndft(kspace, random_image.flatten(), shape) matrix_coef = matrix @ random_image.flatten() - assert np.allclose(linop_coef, matrix_coef) + assert_almost_allclose(linop_coef, matrix_coef, atol=1e-4, rtol=1e-4, mismatch=5) @parametrize_with_cases( @@ -76,7 +53,32 @@ def test_ndft_implicit1(kspace, shape): linop_coef = implicit_type1_ndft(kspace, random_kspace.flatten(), shape) matrix_coef = matrix.conj().T @ random_kspace.flatten() - assert np.allclose(linop_coef, matrix_coef) + assert_almost_allclose(linop_coef, matrix_coef, atol=1e-4, rtol=1e-4, mismatch=5) + +@parametrize_with_cases( + "kspace, shape", + cases=[ + CasesTrajectories.case_random2D, + CasesTrajectories.case_grid2D, + CasesTrajectories.case_grid3D, + ], +) +def test_ndft_nufft(kspace, shape): + "Test that NDFT matches NUFFT" + ndft_op = RawNDFT(kspace, shape, normalize=True) + random_kspace = 1j * np.random.randn(len(kspace)) + random_kspace += np.random.randn(len(kspace)) + random_image = np.random.randn(*shape) + 1j * np.random.randn(*shape) + operator = get_operator("pynfft")(kspace, shape) # FIXME: @PAC, we need to get ref here + nufft_k = operator.op(random_image) + nufft_i = operator.adj_op(random_kspace) + + ndft_k = np.empty(ndft_op.n_samples, dtype=random_image.dtype) + ndft_i = np.empty(shape, dtype=random_kspace.dtype) + ndft_op.op(ndft_k, random_image) + ndft_op.adj_op(random_kspace, ndft_i) + assert_almost_allclose(nufft_k, ndft_k, atol=1e-4, rtol=1e-4, mismatch=5) + assert_almost_allclose(nufft_i, ndft_i, atol=1e-4, rtol=1e-4, mismatch=5) @parametrize_with_cases( @@ -98,6 +100,7 @@ def test_ndft_fft(kspace_grid, shape): kspace = kspace.reshape(img.shape) if len(shape) >= 2: kspace = kspace.swapaxes(0, 1) - kspace_fft = sp.fft.fftn(img) + kspace_fft = sp.fft.fftn(sp.fft.fftshift(img)) + + assert_almost_allclose(kspace, kspace_fft, atol=1e-4, rtol=1e-4, mismatch=5) - assert_almost_allclose(kspace, kspace_fft, atol=1e-5, rtol=1e-5, mismatch=5) From 093edfbe3def12e22c4aec3657f70377c238525b Mon Sep 17 00:00:00 2001 From: chaithyagr Date: Fri, 26 Apr 2024 14:46:50 +0200 Subject: [PATCH 02/16] Fix --- src/mrinufft/operators/interfaces/nudft_numpy.py | 12 ------------ 1 file changed, 12 deletions(-) diff --git a/src/mrinufft/operators/interfaces/nudft_numpy.py b/src/mrinufft/operators/interfaces/nudft_numpy.py index 9d26d02f5..68690cc17 100644 --- a/src/mrinufft/operators/interfaces/nudft_numpy.py +++ b/src/mrinufft/operators/interfaces/nudft_numpy.py @@ -8,17 +8,6 @@ from ..base import FourierOperatorCPU -<<<<<<< Updated upstream -def get_fourier_matrix(ktraj, shape): - """Get the NDFT Fourier Matrix.""" - n = np.prod(shape) - ndim = len(shape) - matrix = np.zeros((len(ktraj), n), dtype=complex) - r = [np.arange(shape[i]) for i in range(ndim)] - grid_r = np.reshape(np.meshgrid(*r, indexing="ij"), (ndim, np.prod(shape))) - traj_grid = ktraj @ grid_r - matrix = np.exp(-2j * np.pi * traj_grid) -======= def get_fourier_matrix(ktraj, shape, dtype=np.complex64, normalize=False): """Get the NDFT Fourier Matrix.""" n = np.prod(shape) @@ -30,7 +19,6 @@ def get_fourier_matrix(ktraj, shape, dtype=np.complex64, normalize=False): matrix = np.exp(-2j * np.pi * traj_grid, dtype=dtype) if normalize: matrix /= (np.sqrt(np.prod(shape)) * np.power(np.sqrt(2), len(shape))) ->>>>>>> Stashed changes return matrix From 74c1ecd6d0a07dcc141156aa506be9d1967d07e5 Mon Sep 17 00:00:00 2001 From: chaithyagr Date: Fri, 26 Apr 2024 14:48:12 +0200 Subject: [PATCH 03/16] Remove bymistake add --- tests/test_autodiff.py | 145 ----------------------------------------- 1 file changed, 145 deletions(-) delete mode 100644 tests/test_autodiff.py diff --git a/tests/test_autodiff.py b/tests/test_autodiff.py deleted file mode 100644 index 01de2f91c..000000000 --- a/tests/test_autodiff.py +++ /dev/null @@ -1,145 +0,0 @@ -"""Test the autodiff functionnality.""" - -import numpy as np -from mrinufft.operators.interfaces.nudft_numpy import get_fourier_matrix -import pytest -from pytest_cases import parametrize_with_cases, parametrize, fixture -from case_trajectories import CasesTrajectories -from mrinufft.operators import get_operator - - -from helpers import ( - kspace_from_op, - image_from_op, - to_interface, -) - - -TORCH_AVAILABLE = True -try: - import torch - import torch.testing as tt -except ImportError: - TORCH_AVAILABLE = False - - -@fixture(scope="module") -@parametrize(backend=["cufinufft", "finufft"]) -@parametrize(autograd=["data"]) -@parametrize_with_cases( - "kspace_loc, shape", - cases=[ - CasesTrajectories.case_grid2D, - CasesTrajectories.case_nyquist_radial2D, - ], # 2D cases only for reduced memory footprint. -) -def operator(kspace_loc, shape, backend, autograd): - """Create NUFFT operator with autodiff capabilities.""" - kspace_loc = kspace_loc.astype(np.float32) - - nufft = get_operator(backend_name=backend, autograd=autograd)( - samples=kspace_loc, - shape=shape, - smaps=None, - ) - - return nufft - - -@fixture(scope="module") -def ndft_matrix(operator): - """Get the NDFT matrix from the operator.""" - return get_fourier_matrix(operator.samples, operator.shape, normalize=True) - - -@pytest.mark.parametrize("interface", ["torch-gpu", "torch-cpu"]) -@pytest.mark.skipif(not TORCH_AVAILABLE, reason="Pytorch is not installed") -def test_adjoint_and_grad(operator, ndft_matrix, interface): - """Test the adjoint and gradient of the operator.""" - if operator.backend == "finufft" and "gpu" in interface: - pytest.skip("GPU not supported for finufft backend") - ndft_matrix_torch = to_interface(ndft_matrix, interface=interface) - ksp_data = to_interface(kspace_from_op(operator), interface=interface) - img_data = to_interface(image_from_op(operator), interface=interface) - ksp_data.requires_grad = True - - with torch.autograd.set_detect_anomaly(True): - adj_data = operator.adj_op(ksp_data).reshape(img_data.shape) - adj_data_ndft = (ndft_matrix_torch.conj().T @ ksp_data.flatten()).reshape( - adj_data.shape - ) - loss_nufft = torch.mean(torch.abs(adj_data) ** 2) - loss_ndft = torch.mean(torch.abs(adj_data_ndft) ** 2) - - # Check if nufft and ndft are close in the backprop - grad_ndft_kdata = torch.autograd.grad(loss_ndft, ksp_data, retain_graph=True)[0] - grad_nufft_kdata = torch.autograd.grad(loss_nufft, ksp_data, retain_graph=True)[0] - tt.assert_close(grad_ndft_kdata, grad_nufft_kdata, rtol=1, atol=1) - - -@pytest.mark.parametrize("interface", ["torch-gpu", "torch-cpu"]) -@pytest.mark.skipif(not TORCH_AVAILABLE, reason="Pytorch is not installed") -def test_adjoint_and_gradauto(operator, ndft_matrix, interface): - """Test the adjoint and gradient of the operator using autograd gradcheck.""" - if operator.backend == "finufft" and "gpu" in interface: - pytest.skip("GPU not supported for finufft backend") - - ksp_data = to_interface(kspace_from_op(operator), interface=interface) - ksp_data = torch.ones(ksp_data.shape, requires_grad=True, dtype=ksp_data.dtype) - print(ksp_data.shape) - # todo: tighten the tolerance - assert torch.autograd.gradcheck( - operator.adjoint, - ksp_data, - eps=1e-10, - rtol=1, - atol=1, - nondet_tol=1, - raise_exception=True, - ) - - -@pytest.mark.parametrize("interface", ["torch-gpu", "torch-cpu"]) -@pytest.mark.skipif(not TORCH_AVAILABLE, reason="Pytorch is not installed") -def test_forward_and_grad(operator, ndft_matrix, interface): - """Test the adjoint and gradient of the operator.""" - if operator.backend == "finufft" and "gpu" in interface: - pytest.skip("GPU not supported for finufft backend") - - ndft_matrix_torch = to_interface(ndft_matrix, interface=interface) - ksp_data_ref = to_interface(kspace_from_op(operator), interface=interface) - img_data = to_interface(image_from_op(operator), interface=interface) - img_data.requires_grad = True - - with torch.autograd.set_detect_anomaly(True): - ksp_data = operator.op(img_data).reshape(ksp_data_ref.shape) - ksp_data_ndft = (ndft_matrix_torch @ img_data.flatten()).reshape(ksp_data.shape) - loss_nufft = torch.mean(torch.abs(ksp_data - ksp_data_ref) ** 2) - loss_ndft = torch.mean(torch.abs(ksp_data_ndft - ksp_data_ref) ** 2) - - # Check if nufft and ndft are close in the backprop - grad_ndft_kdata = torch.autograd.grad(loss_ndft, img_data, retain_graph=True)[0] - grad_nufft_kdata = torch.autograd.grad(loss_nufft, img_data, retain_graph=True)[0] - assert torch.allclose(grad_ndft_kdata, grad_nufft_kdata, atol=6e-3) - - -@pytest.mark.parametrize("interface", ["torch-gpu", "torch-cpu"]) -@pytest.mark.skipif(not TORCH_AVAILABLE, reason="Pytorch is not installed") -def test_forward_and_gradauto(operator, ndft_matrix, interface): - """Test the forward and gradient of the operator using autograd gradcheck.""" - if operator.backend == "finufft" and "gpu" in interface: - pytest.skip("GPU not supported for finufft backend") - - img_data = to_interface(image_from_op(operator), interface=interface) - img_data = torch.ones(img_data.shape, requires_grad=True, dtype=img_data.dtype) - print(img_data.shape) - # todo: tighten the tolerance - assert torch.autograd.gradcheck( - operator.adjoint, - img_data, - eps=1e-10, - rtol=1, - atol=1, - nondet_tol=1, - raise_exception=True, - ) From 0250aa8d3753bad0191cdc5f42cd1c112f589f44 Mon Sep 17 00:00:00 2001 From: chaithyagr Date: Fri, 26 Apr 2024 15:38:27 +0200 Subject: [PATCH 04/16] Fix --- tests/test_ndft.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/test_ndft.py b/tests/test_ndft.py index 5bcb3e8c9..3d972ff56 100644 --- a/tests/test_ndft.py +++ b/tests/test_ndft.py @@ -64,7 +64,7 @@ def test_ndft_implicit1(kspace, shape): ], ) def test_ndft_nufft(kspace, shape): - "Test that NDFT matches NUFFT" + """Test that NDFT matches NUFFT""" ndft_op = RawNDFT(kspace, shape, normalize=True) random_kspace = 1j * np.random.randn(len(kspace)) random_kspace += np.random.randn(len(kspace)) From 060a8bdd125d2140b82dfaa6a492ec78682a80bf Mon Sep 17 00:00:00 2001 From: chaithyagr Date: Fri, 26 Apr 2024 15:39:16 +0200 Subject: [PATCH 05/16] Fixed lint --- .../operators/interfaces/nudft_numpy.py | 20 +++++++++++-------- tests/test_ndft.py | 9 +++++---- 2 files changed, 17 insertions(+), 12 deletions(-) diff --git a/src/mrinufft/operators/interfaces/nudft_numpy.py b/src/mrinufft/operators/interfaces/nudft_numpy.py index 68690cc17..3e8e81aa3 100644 --- a/src/mrinufft/operators/interfaces/nudft_numpy.py +++ b/src/mrinufft/operators/interfaces/nudft_numpy.py @@ -13,18 +13,18 @@ def get_fourier_matrix(ktraj, shape, dtype=np.complex64, normalize=False): n = np.prod(shape) ndim = len(shape) matrix = np.zeros((len(ktraj), n), dtype=dtype) - r = [np.linspace(-s/2, s/2-1, s) for s in shape] + r = [np.linspace(-s / 2, s / 2 - 1, s) for s in shape] grid_r = np.reshape(np.meshgrid(*r, indexing="ij"), (ndim, np.prod(shape))) traj_grid = ktraj @ grid_r matrix = np.exp(-2j * np.pi * traj_grid, dtype=dtype) if normalize: - matrix /= (np.sqrt(np.prod(shape)) * np.power(np.sqrt(2), len(shape))) + matrix /= np.sqrt(np.prod(shape)) * np.power(np.sqrt(2), len(shape)) return matrix def implicit_type2_ndft(ktraj, image, shape, normalize=False): """Compute the NDFT using the implicit type 2 (image -> kspace) algorithm.""" - r = [np.linspace(-s/2, s/2-1, s) for s in shape] + r = [np.linspace(-s / 2, s / 2 - 1, s) for s in shape] grid_r = np.reshape( np.meshgrid(*r, indexing="ij"), (len(shape), np.prod(image.shape)) ) @@ -32,19 +32,19 @@ def implicit_type2_ndft(ktraj, image, shape, normalize=False): for j in range(np.prod(image.shape)): res += image[j] * np.exp(-2j * np.pi * ktraj @ grid_r[:, j]) if normalize: - matrix /= (np.sqrt(np.prod(shape)) * np.power(np.sqrt(2), len(shape))) + matrix /= np.sqrt(np.prod(shape)) * np.power(np.sqrt(2), len(shape)) return res def implicit_type1_ndft(ktraj, coeffs, shape, normalize=False): """Compute the NDFT using the implicit type 1 (kspace -> image) algorithm.""" - r = [np.linspace(-s/2, s/2-1, s) for s in shape] + r = [np.linspace(-s / 2, s / 2 - 1, s) for s in shape] grid_r = np.reshape(np.meshgrid(*r, indexing="ij"), (len(shape), np.prod(shape))) res = np.zeros(np.prod(shape), dtype=coeffs.dtype) for i in range(len(ktraj)): res += coeffs[i] * np.exp(2j * np.pi * ktraj[i] @ grid_r) if normalize: - matrix /= (np.sqrt(np.prod(shape)) * np.power(np.sqrt(2), len(shape))) + matrix /= np.sqrt(np.prod(shape)) * np.power(np.sqrt(2), len(shape)) return res @@ -75,9 +75,13 @@ def __init__(self, samples, shape, explicit_matrix=True, normalize=False): ) except MemoryError: warnings.warn("Not enough memory, using an implicit definition anyway") - self._fourier_matrix = get_implicit_matrix(self.samples, self.shape, normalize) + self._fourier_matrix = get_implicit_matrix( + self.samples, self.shape, normalize + ) else: - self._fourier_matrix = get_implicit_matrix(self.samples, self.shape, normalize) + self._fourier_matrix = get_implicit_matrix( + self.samples, self.shape, normalize + ) def op(self, coeffs, image): """Compute the forward NUDFT.""" diff --git a/tests/test_ndft.py b/tests/test_ndft.py index 3d972ff56..7f90d14ea 100644 --- a/tests/test_ndft.py +++ b/tests/test_ndft.py @@ -16,7 +16,6 @@ from mrinufft import get_operator - @parametrize_with_cases( "kspace, shape", cases=[ @@ -55,6 +54,7 @@ def test_ndft_implicit1(kspace, shape): assert_almost_allclose(linop_coef, matrix_coef, atol=1e-4, rtol=1e-4, mismatch=5) + @parametrize_with_cases( "kspace, shape", cases=[ @@ -69,10 +69,12 @@ def test_ndft_nufft(kspace, shape): random_kspace = 1j * np.random.randn(len(kspace)) random_kspace += np.random.randn(len(kspace)) random_image = np.random.randn(*shape) + 1j * np.random.randn(*shape) - operator = get_operator("pynfft")(kspace, shape) # FIXME: @PAC, we need to get ref here + operator = get_operator("pynfft")( + kspace, shape + ) # FIXME: @PAC, we need to get ref here nufft_k = operator.op(random_image) nufft_i = operator.adj_op(random_kspace) - + ndft_k = np.empty(ndft_op.n_samples, dtype=random_image.dtype) ndft_i = np.empty(shape, dtype=random_kspace.dtype) ndft_op.op(ndft_k, random_image) @@ -103,4 +105,3 @@ def test_ndft_fft(kspace_grid, shape): kspace_fft = sp.fft.fftn(sp.fft.fftshift(img)) assert_almost_allclose(kspace, kspace_fft, atol=1e-4, rtol=1e-4, mismatch=5) - From aecb844c74ae53ae67deb852204bc9e647ac28fd Mon Sep 17 00:00:00 2001 From: chaithyagr Date: Fri, 26 Apr 2024 15:40:50 +0200 Subject: [PATCH 06/16] Lint --- src/mrinufft/operators/interfaces/nudft_numpy.py | 4 ++-- tests/test_ndft.py | 2 +- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/src/mrinufft/operators/interfaces/nudft_numpy.py b/src/mrinufft/operators/interfaces/nudft_numpy.py index 3e8e81aa3..bcc6c0338 100644 --- a/src/mrinufft/operators/interfaces/nudft_numpy.py +++ b/src/mrinufft/operators/interfaces/nudft_numpy.py @@ -32,7 +32,7 @@ def implicit_type2_ndft(ktraj, image, shape, normalize=False): for j in range(np.prod(image.shape)): res += image[j] * np.exp(-2j * np.pi * ktraj @ grid_r[:, j]) if normalize: - matrix /= np.sqrt(np.prod(shape)) * np.power(np.sqrt(2), len(shape)) + res /= np.sqrt(np.prod(shape)) * np.power(np.sqrt(2), len(shape)) return res @@ -44,7 +44,7 @@ def implicit_type1_ndft(ktraj, coeffs, shape, normalize=False): for i in range(len(ktraj)): res += coeffs[i] * np.exp(2j * np.pi * ktraj[i] @ grid_r) if normalize: - matrix /= np.sqrt(np.prod(shape)) * np.power(np.sqrt(2), len(shape)) + res /= np.sqrt(np.prod(shape)) * np.power(np.sqrt(2), len(shape)) return res diff --git a/tests/test_ndft.py b/tests/test_ndft.py index 7f90d14ea..fa66b8b26 100644 --- a/tests/test_ndft.py +++ b/tests/test_ndft.py @@ -64,7 +64,7 @@ def test_ndft_implicit1(kspace, shape): ], ) def test_ndft_nufft(kspace, shape): - """Test that NDFT matches NUFFT""" + """Test that NDFT matches NUFFT.""" ndft_op = RawNDFT(kspace, shape, normalize=True) random_kspace = 1j * np.random.randn(len(kspace)) random_kspace += np.random.randn(len(kspace)) From 3130bc1c5f443294a2f71dcae30178bb8357d392 Mon Sep 17 00:00:00 2001 From: chaithyagr Date: Fri, 26 Apr 2024 17:15:00 +0200 Subject: [PATCH 07/16] Added refbackend --- tests/test_ndft.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/tests/test_ndft.py b/tests/test_ndft.py index fa66b8b26..57aedfa6d 100644 --- a/tests/test_ndft.py +++ b/tests/test_ndft.py @@ -63,18 +63,18 @@ def test_ndft_implicit1(kspace, shape): CasesTrajectories.case_grid3D, ], ) -def test_ndft_nufft(kspace, shape): +def test_ndft_nufft(kspace, shape, request): """Test that NDFT matches NUFFT.""" ndft_op = RawNDFT(kspace, shape, normalize=True) random_kspace = 1j * np.random.randn(len(kspace)) random_kspace += np.random.randn(len(kspace)) random_image = np.random.randn(*shape) + 1j * np.random.randn(*shape) - operator = get_operator("pynfft")( + operator = get_operator(request.config.getoption("ref"))( kspace, shape - ) # FIXME: @PAC, we need to get ref here + ) nufft_k = operator.op(random_image) nufft_i = operator.adj_op(random_kspace) - + ndft_k = np.empty(ndft_op.n_samples, dtype=random_image.dtype) ndft_i = np.empty(shape, dtype=random_kspace.dtype) ndft_op.op(ndft_k, random_image) From bc014b8973e3b355f17365a9eb933cc57b92fb4b Mon Sep 17 00:00:00 2001 From: chaithyagr Date: Fri, 26 Apr 2024 17:17:48 +0200 Subject: [PATCH 08/16] Fix NDFT --- tests/test_ndft.py | 6 ++---- 1 file changed, 2 insertions(+), 4 deletions(-) diff --git a/tests/test_ndft.py b/tests/test_ndft.py index 57aedfa6d..7a157d349 100644 --- a/tests/test_ndft.py +++ b/tests/test_ndft.py @@ -69,12 +69,10 @@ def test_ndft_nufft(kspace, shape, request): random_kspace = 1j * np.random.randn(len(kspace)) random_kspace += np.random.randn(len(kspace)) random_image = np.random.randn(*shape) + 1j * np.random.randn(*shape) - operator = get_operator(request.config.getoption("ref"))( - kspace, shape - ) + operator = get_operator(request.config.getoption("ref"))(kspace, shape) nufft_k = operator.op(random_image) nufft_i = operator.adj_op(random_kspace) - + ndft_k = np.empty(ndft_op.n_samples, dtype=random_image.dtype) ndft_i = np.empty(shape, dtype=random_kspace.dtype) ndft_op.op(ndft_k, random_image) From 0cc73c41cf743ea19ffa053f0cd43b54f43f192e Mon Sep 17 00:00:00 2001 From: Pierre-antoine Comby Date: Mon, 29 Apr 2024 10:48:25 +0200 Subject: [PATCH 09/16] feat: use finufft as ref backend. --- tests/conftest.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/conftest.py b/tests/conftest.py index 69598fdb4..4e89f0ed5 100644 --- a/tests/conftest.py +++ b/tests/conftest.py @@ -15,7 +15,7 @@ def pytest_addoption(parser): ) parser.addoption( "--ref", - default="pynfft", + default="finufft", help="Reference backend on which the tests are performed.", ) From 21e090f21803e9e57cc721010661d30449ce1a0b Mon Sep 17 00:00:00 2001 From: Pierre-antoine Comby Date: Mon, 29 Apr 2024 10:49:37 +0200 Subject: [PATCH 10/16] feat(tests): move ndft vs nufft tests to own file. --- tests/operators/test_operator_ref.py | 74 ++++++++++++++++++++++++++++ tests/test_ndft.py | 26 ---------- 2 files changed, 74 insertions(+), 26 deletions(-) create mode 100644 tests/operators/test_operator_ref.py diff --git a/tests/operators/test_operator_ref.py b/tests/operators/test_operator_ref.py new file mode 100644 index 000000000..b51e1633b --- /dev/null +++ b/tests/operators/test_operator_ref.py @@ -0,0 +1,74 @@ +"""Tests for the reference backend.""" + +from pytest_cases import parametrize_with_cases, fixture +from case_trajectories import CasesTrajectories + +from mrinufft import get_operator +from mrinufft.operators.interfaces.nudft_numpy import MRInumpy +from helpers import assert_almost_allclose, kspace_from_op, image_from_op + + +@fixture(scope="session", autouse=True) +def ref_backend(request): + """Get the reference backend from the CLI.""" + return request.config.getoption("ref") + + +@fixture(scope="module") +@parametrize_with_cases( + "kspace, shape", + cases=[ + CasesTrajectories.case_random2D, + CasesTrajectories.case_grid2D, + CasesTrajectories.case_grid3D, + ], +) +def ref_operator(request, ref_backend, kspace, shape): + """Generate a NFFT operator, matching the property of the first operator.""" + return get_operator(ref_backend)(kspace, shape) + + +@fixture(scope="module") +def ndft_operator(ref_operator): + """Get a NDFT operator matching the reference operator.""" + return MRInumpy(ref_operator.samples, ref_operator.shape) + + +@fixture(scope="module") +def image_data(ref_operator): + """Generate a random image. Remains constant for the module.""" + return image_from_op(ref_operator) + + +@fixture(scope="module") +def kspace_data(ref_operator): + """Generate a random kspace. Remains constant for the module.""" + return kspace_from_op(ref_operator) + + +def test_ref_nufft_forward(ref_operator, ndft_operator, image_data): + """Test that the reference nufft matches the NDFT.""" + nufft_ksp = ref_operator.op(image_data) + ndft_ksp = ndft_operator.op(image_data) + + assert_almost_allclose( + nufft_ksp, + ndft_ksp, + atol=1e-4, + rtol=1e-4, + mismatch=5, + ) + + +def test_ref_nufft_adjoint(ref_operator, ndft_operator, kspace_data): + """Test that the reference nufft matches the NDFT adjoint.""" + nufft_img = ref_operator.adj_op(kspace_data) + ndft_img = ndft_operator.adj_op(kspace_data) + + assert_almost_allclose( + nufft_img, + ndft_img, + atol=1e-4, + rtol=1e-4, + mismatch=5, + ) diff --git a/tests/test_ndft.py b/tests/test_ndft.py index 7a157d349..cd21622ea 100644 --- a/tests/test_ndft.py +++ b/tests/test_ndft.py @@ -55,32 +55,6 @@ def test_ndft_implicit1(kspace, shape): assert_almost_allclose(linop_coef, matrix_coef, atol=1e-4, rtol=1e-4, mismatch=5) -@parametrize_with_cases( - "kspace, shape", - cases=[ - CasesTrajectories.case_random2D, - CasesTrajectories.case_grid2D, - CasesTrajectories.case_grid3D, - ], -) -def test_ndft_nufft(kspace, shape, request): - """Test that NDFT matches NUFFT.""" - ndft_op = RawNDFT(kspace, shape, normalize=True) - random_kspace = 1j * np.random.randn(len(kspace)) - random_kspace += np.random.randn(len(kspace)) - random_image = np.random.randn(*shape) + 1j * np.random.randn(*shape) - operator = get_operator(request.config.getoption("ref"))(kspace, shape) - nufft_k = operator.op(random_image) - nufft_i = operator.adj_op(random_kspace) - - ndft_k = np.empty(ndft_op.n_samples, dtype=random_image.dtype) - ndft_i = np.empty(shape, dtype=random_kspace.dtype) - ndft_op.op(ndft_k, random_image) - ndft_op.adj_op(random_kspace, ndft_i) - assert_almost_allclose(nufft_k, ndft_k, atol=1e-4, rtol=1e-4, mismatch=5) - assert_almost_allclose(nufft_i, ndft_i, atol=1e-4, rtol=1e-4, mismatch=5) - - @parametrize_with_cases( "kspace_grid, shape", cases=[ From df10d86f10489f50f10bf8042e24dc6c57d9b038 Mon Sep 17 00:00:00 2001 From: chaithyagr Date: Mon, 3 Feb 2025 13:05:29 +0100 Subject: [PATCH 11/16] Added a bunch of extra codes --- src/mrinufft/io/siemens.py | 8 +++++++- src/mrinufft/io/utils.py | 24 ++++++++++++++++++++++++ 2 files changed, 31 insertions(+), 1 deletion(-) diff --git a/src/mrinufft/io/siemens.py b/src/mrinufft/io/siemens.py index b336ca457..ab56b1daa 100644 --- a/src/mrinufft/io/siemens.py +++ b/src/mrinufft/io/siemens.py @@ -75,7 +75,12 @@ def read_siemens_rawdat( "n_slices": int(twixObj.image.NSli), "n_average": int(twixObj.image.NAve), "orientation": siemens_quat_to_rot_mat(twixObj.image.slicePos[0][-4:]), + "acs": None, } + if "refscan" in twixObj.keys(): + twixObj.refscan.squeeze = True + acs = twixObj.refscan[""].astype(np.float32) + hdr["acs"] = acs.swapaxes(0, 1) if slice_num is not None and hdr["n_slices"] < slice_num: raise ValueError("The slice number is out of bounds.") if contrast_num is not None and hdr["n_contrasts"] < contrast_num: @@ -97,7 +102,8 @@ def read_siemens_rawdat( data = data.reshape( hdr["n_coils"], - hdr["n_shots"] * hdr["n_adc_samples"], + hdr["n_shots"], + hdr["n_adc_samples"], hdr["n_slices"] if slice_num is None else 1, hdr["n_contrasts"] if contrast_num is None else 1, hdr["n_average"] if hdr["n_average"] > 1 and not doAverage else 1, diff --git a/src/mrinufft/io/utils.py b/src/mrinufft/io/utils.py index 7c6433b71..552dccbcf 100644 --- a/src/mrinufft/io/utils.py +++ b/src/mrinufft/io/utils.py @@ -60,3 +60,27 @@ def siemens_quat_to_rot_mat(quat): R[2] = -R[2] R[-1, -1] = 1 return R + +def remove_extra_kspace_samples(kspace_data, num_samples_per_shot): + """ + Remove extra samples from k-space data. + This function is useful when the k-space data has extra samples + mainly as ADC samples at only upto + + Parameters + ---------- + kspace_data : np.ndarray + The k-space data ordered as NCha X NShot X NSamples. + num_samples_per_shot : int + The number of samples per shot in trajectory + + Returns + ------- + np.ndarray + The k-space data with extra samples removed. + """ + n_samples = kspace_data.shape[-1] + n_extra_samples = n_samples - num_samples_per_shot + if n_extra_samples > 0: + kspace_data = kspace_data[..., :-n_extra_samples] + return kspace_data From f280aae5f92e445880772648a3cbab92b8991de3 Mon Sep 17 00:00:00 2001 From: chaithyagr Date: Mon, 3 Feb 2025 13:06:09 +0100 Subject: [PATCH 12/16] PEP fixes --- src/mrinufft/io/utils.py | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/src/mrinufft/io/utils.py b/src/mrinufft/io/utils.py index 552dccbcf..2366bfa65 100644 --- a/src/mrinufft/io/utils.py +++ b/src/mrinufft/io/utils.py @@ -61,9 +61,10 @@ def siemens_quat_to_rot_mat(quat): R[-1, -1] = 1 return R + def remove_extra_kspace_samples(kspace_data, num_samples_per_shot): - """ - Remove extra samples from k-space data. + """Remove extra samples from k-space data. + This function is useful when the k-space data has extra samples mainly as ADC samples at only upto From 9e8010e2c1ba537282233ae4c274f47dea1fd65f Mon Sep 17 00:00:00 2001 From: Chaithya G R Date: Tue, 4 Feb 2025 21:02:01 +0100 Subject: [PATCH 13/16] Update siemens.py --- src/mrinufft/io/siemens.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/mrinufft/io/siemens.py b/src/mrinufft/io/siemens.py index ab56b1daa..b828caa09 100644 --- a/src/mrinufft/io/siemens.py +++ b/src/mrinufft/io/siemens.py @@ -79,7 +79,7 @@ def read_siemens_rawdat( } if "refscan" in twixObj.keys(): twixObj.refscan.squeeze = True - acs = twixObj.refscan[""].astype(np.float32) + acs = twixObj.refscan[""].astype(np.complex64) hdr["acs"] = acs.swapaxes(0, 1) if slice_num is not None and hdr["n_slices"] < slice_num: raise ValueError("The slice number is out of bounds.") From 53cad71c1e99c9d3648539572b248f493ee0013f Mon Sep 17 00:00:00 2001 From: chaithyagr Date: Wed, 5 Feb 2025 10:12:45 +0100 Subject: [PATCH 14/16] Added fixes --- src/mrinufft/io/siemens.py | 1 + src/mrinufft/io/utils.py | 4 +++- 2 files changed, 4 insertions(+), 1 deletion(-) diff --git a/src/mrinufft/io/siemens.py b/src/mrinufft/io/siemens.py index b828caa09..0b1d3af65 100644 --- a/src/mrinufft/io/siemens.py +++ b/src/mrinufft/io/siemens.py @@ -42,6 +42,7 @@ def read_siemens_rawdat( Imported data formatted as n_coils X n_samples X n_slices X n_contrasts hdr: dict Extra information about the data parsed from the twix file + This header also contains the ACS data as "acs" if it was found in raw data. Raises ------ diff --git a/src/mrinufft/io/utils.py b/src/mrinufft/io/utils.py index 2366bfa65..428a81a00 100644 --- a/src/mrinufft/io/utils.py +++ b/src/mrinufft/io/utils.py @@ -66,7 +66,9 @@ def remove_extra_kspace_samples(kspace_data, num_samples_per_shot): """Remove extra samples from k-space data. This function is useful when the k-space data has extra samples - mainly as ADC samples at only upto + mainly as ADC samples at only at specific number of samples. + This sometimes leads to a situation where we will have more ADC samples + than what is expected. Parameters ---------- From 12c2c6a111a901a7c7a426f984e9273063a7cc12 Mon Sep 17 00:00:00 2001 From: chaithyagr Date: Wed, 5 Feb 2025 10:17:13 +0100 Subject: [PATCH 15/16] add [docs] --- src/mrinufft/io/utils.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/mrinufft/io/utils.py b/src/mrinufft/io/utils.py index 428a81a00..dd2230fd5 100644 --- a/src/mrinufft/io/utils.py +++ b/src/mrinufft/io/utils.py @@ -67,8 +67,8 @@ def remove_extra_kspace_samples(kspace_data, num_samples_per_shot): This function is useful when the k-space data has extra samples mainly as ADC samples at only at specific number of samples. - This sometimes leads to a situation where we will have more ADC samples - than what is expected. + This sometimes leads to a situation where we will have more ADC samples + than what is expected. Parameters ---------- From 3f87e64f68788242db3920ef32758eec30b45723 Mon Sep 17 00:00:00 2001 From: chaithyagr Date: Wed, 5 Feb 2025 13:32:44 +0100 Subject: [PATCH 16/16] Fixes and updates on the locatuions --- src/mrinufft/io/nsp.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/mrinufft/io/nsp.py b/src/mrinufft/io/nsp.py index b0594e2e3..dfe86152e 100644 --- a/src/mrinufft/io/nsp.py +++ b/src/mrinufft/io/nsp.py @@ -484,7 +484,7 @@ def read_arbgrad_rawdat( if "ARBGRAD_VE11C" in data_type: hdr["type"] = "ARBGRAD_GRE" hdr["shifts"] = () - for s in [7, 6, 8]: + for s in [6, 7, 8]: shift = twixObj.search_header_for_val( "Phoenix", ("sWiPMemBlock", "adFree", str(s)) )