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.

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

FourierOp Norm #542

Draft
wants to merge 7 commits into
base: main
Choose a base branch
from
Draft
Show file tree
Hide file tree
Changes from all commits
Commits
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
4 changes: 4 additions & 0 deletions src/mrpro/data/SpatialDimension.py
Original file line number Diff line number Diff line change
Expand Up @@ -414,3 +414,7 @@ def shape(self) -> tuple[int, ...]:
return self.x.shape
else:
raise ValueError('Inconsistent shapes')

def prod(self: SpatialDimension[T_co]) -> T_co:
"""Calculate the product x*y*z."""
return self.x * self.y * self.z
7 changes: 7 additions & 0 deletions src/mrpro/operators/DensityCompensationOp.py
Original file line number Diff line number Diff line change
@@ -1,6 +1,7 @@
"""Class for Density Compensation Operator."""

import torch
from typing_extensions import Self

from mrpro.data.DcfData import DcfData
from mrpro.operators.EinsumOp import EinsumOp
Expand All @@ -25,3 +26,9 @@ def __init__(self, dcf: DcfData | torch.Tensor) -> None:
else:
dcf_tensor = dcf
super().__init__(dcf_tensor, '... k2 k1 k0 ,... coil k2 k1 k0 ->... coil k2 k1 k0')

def __pow__(self, exponent: float) -> Self:
"""Raise the operator to a power."""
if not exponent > 0:
raise NotImplementedError('Can only raise to positive powers')
return type(self)(dcf=self.matrix**exponent)
4 changes: 2 additions & 2 deletions src/mrpro/operators/FastFourierOp.py
Original file line number Diff line number Diff line change
Expand Up @@ -132,10 +132,10 @@ def adjoint(self, y: torch.Tensor) -> tuple[torch.Tensor,]:
-------
IFFT of y
"""
# FFT
return self._pad_op.adjoint(
(x,) = self._pad_op.adjoint(
torch.fft.fftshift(
torch.fft.ifftn(torch.fft.ifftshift(y, dim=self._dim), dim=self._dim, norm='ortho'),
dim=self._dim,
),
)
return (x,)
13 changes: 8 additions & 5 deletions src/mrpro/operators/FourierOp.py
Original file line number Diff line number Diff line change
Expand Up @@ -126,6 +126,7 @@ def get_traj(traj: KTrajectory, dims: Sequence[int]):
self._fwd_nufft_op = None
self._adj_nufft_op = None
self._kshape = traj.broadcasted_shape
self._nufft_scale = torch.tensor(get_spatial_dims(encoding_matrix, self._nufft_dims)).prod().sqrt().reciprocal()

@classmethod
def from_kdata(cls, kdata: KData, recon_shape: SpatialDimension[int] | None = None) -> Self:
Expand Down Expand Up @@ -164,7 +165,6 @@ def forward(self, x: torch.Tensor) -> tuple[torch.Tensor,]:
keep_dims = [-4, *self._nufft_dims] # -4 is always coil
permute = [i for i in range(-x.ndim, 0) if i not in keep_dims] + keep_dims
unpermute = np.argsort(permute)

x = x.permute(*permute)
permuted_x_shape = x.shape
x = x.flatten(end_dim=-len(keep_dims) - 1)
Expand All @@ -174,16 +174,18 @@ def forward(self, x: torch.Tensor) -> tuple[torch.Tensor,]:
omega = omega.broadcast_to(*permuted_x_shape[: -len(keep_dims)], *omega.shape[-len(keep_dims) :])
omega = omega.flatten(end_dim=-len(keep_dims) - 1).flatten(start_dim=-len(keep_dims) + 1)

x = self._fwd_nufft_op(x, omega, norm='ortho')
x = self._fwd_nufft_op(x, omega, norm=None)

shape_nufft_dims = [self._kshape[i] for i in self._nufft_dims]
x = x.reshape(*permuted_x_shape[: -len(keep_dims)], -1, *shape_nufft_dims) # -1 is coils
x = x.permute(*unpermute)

# manual scaling instead of ortho to be independent of oversampling ratio
x = x * self._nufft_scale

if self._fast_fourier_op is not None and self._cart_sampling_op is not None:
# FFT
(x,) = self._cart_sampling_op(self._fast_fourier_op(x)[0])

return (x,)

def adjoint(self, x: torch.Tensor) -> tuple[torch.Tensor,]:
Expand All @@ -209,7 +211,6 @@ def adjoint(self, x: torch.Tensor) -> tuple[torch.Tensor,]:
keep_dims = [-4, *self._nufft_dims] # -4 is coil
permute = [i for i in range(-x.ndim, 0) if i not in keep_dims] + keep_dims
unpermute = np.argsort(permute)

x = x.permute(*permute)
permuted_x_shape = x.shape
x = x.flatten(end_dim=-len(keep_dims) - 1).flatten(start_dim=-len(keep_dims) + 1)
Expand All @@ -218,11 +219,13 @@ def adjoint(self, x: torch.Tensor) -> tuple[torch.Tensor,]:
omega = omega.broadcast_to(*permuted_x_shape[: -len(keep_dims)], *omega.shape[-len(keep_dims) :])
omega = omega.flatten(end_dim=-len(keep_dims) - 1).flatten(start_dim=-len(keep_dims) + 1)

x = self._adj_nufft_op(x, omega, norm='ortho')
x = self._adj_nufft_op(x, omega, norm=None)

x = x.reshape(*permuted_x_shape[: -len(keep_dims)], *x.shape[-len(keep_dims) :])
x = x.permute(*unpermute)

# manual scaling instead of ortho to be independent of oversampling ratio
x = x * self._nufft_scale
return (x,)

@property
Expand Down
92 changes: 53 additions & 39 deletions tests/conftest.py
Original file line number Diff line number Diff line change
Expand Up @@ -240,7 +240,7 @@ def create_traj(k_shape, nkx, nky, nkz, type_kx, type_ky, type_kz):
k_list = []
for spacing, nk in zip([type_kz, type_ky, type_kx], [nkz, nky, nkx], strict=True):
if spacing == 'non-uniform':
k = random_generator.float32_tensor(size=nk, low=-1, high=1) * max(nk)
k = random_generator.float32_tensor(size=nk, low=-max(nk) // 2, high=max(nk) // 2)
elif spacing == 'uniform':
k = create_uniform_traj(nk, k_shape=k_shape)
elif spacing == 'zero':
Expand Down Expand Up @@ -293,10 +293,10 @@ def ismrmrd_cart(ellipse_phantom, tmp_path_factory):
'zero', # type_k2
),
( # (2) 2d non-Cartesian mri with 2 coils
(1, 2, 1, 96, 128), # im_shape
(1, 2, 1, 16, 192), # k_shape
(1, 1, 16, 192), # nkx
(1, 1, 16, 192), # nky
(1, 2, 1, 128, 128), # im_shape
(1, 2, 1, 16, 128), # k_shape
(1, 1, 16, 128), # nkx
(1, 1, 16, 128), # nky
(1, 1, 1, 1), # nkz
'non-uniform', # type_kx
'non-uniform', # type_ky
Expand All @@ -306,7 +306,7 @@ def ismrmrd_cart(ellipse_phantom, tmp_path_factory):
'zero', # type_k2
),
( # (3) 2d Cartesian with irregular sampling
(1, 1, 1, 96, 128), # im_shape
(1, 1, 1, 192, 192), # im_shape
(1, 1, 1, 1, 192), # k_shape
(1, 1, 1, 192), # nkx
(1, 1, 1, 192), # nky
Expand Down Expand Up @@ -344,71 +344,84 @@ def ismrmrd_cart(ellipse_phantom, tmp_path_factory):
'non-uniform', # type_k1
'non-uniform', # type_k2
),
( # (6) 2d non-uniform cine with 8 cardiac phases, 5 coils
(8, 5, 1, 64, 64), # im_shape
(8, 5, 1, 18, 128), # k_shape
(8, 1, 18, 128), # nkx
(8, 1, 18, 128), # nky
(8, 1, 1, 1), # nkz
( # (6) 2d non-uniform cine with 2 cardiac phases, 3 coils
(2, 3, 1, 64, 64), # im_shape
(2, 3, 1, 18, 128), # k_shape
(2, 1, 18, 128), # nkx
(2, 1, 18, 128), # nky
(2, 1, 1, 1), # nkz
'non-uniform', # type_kx
'non-uniform', # type_ky
'zero', # type_kz
'non-uniform', # type_k0
'non-uniform', # type_k1
'zero', # type_k2
),
( # (7) 2d cartesian cine with 9 cardiac phases, 6 coils
(9, 6, 1, 96, 128), # im_shape
(9, 6, 1, 128, 192), # k_shape
(9, 1, 1, 192), # nkx
(9, 1, 128, 1), # nky
(9, 1, 1, 1), # nkz
( # (7) 2d cartesian cine with 2 cardiac phases, 3 coils
(2, 3, 1, 96, 128), # im_shape
(2, 3, 1, 128, 192), # k_shape
(2, 1, 1, 192), # nkx
(2, 1, 128, 1), # nky
(2, 1, 1, 1), # nkz
'uniform', # type_kx
'uniform', # type_ky
'zero', # type_kz
'uniform', # type_k0
'uniform', # type_k1
'zero', # type_k2
),
( # (8) radial phase encoding (RPE), 8 coils, with oversampling in both FFT and non-uniform directions
(2, 8, 64, 32, 48), # im_shape
(2, 8, 8, 64, 96), # k_shape
( # (8) radial phase encoding (RPE), 3 coils, with oversampling in both FFT and non-uniform directions
(2, 3, 64, 32, 48), # im_shape
(2, 3, 8, 64, 96), # k_shape
(2, 1, 1, 96), # nkx
(2, 8, 64, 1), # nky
(2, 8, 64, 1), # nkz
(2, 3, 64, 1), # nky
(2, 3, 64, 1), # nkz
'uniform', # type_kx
'non-uniform', # type_ky
'non-uniform', # type_kz
'uniform', # type_k0
'non-uniform', # type_k1
'non-uniform', # type_k2
),
( # (9) radial phase encoding (RPE), 8 coils with non-Cartesian sampling along readout
(2, 8, 64, 32, 48), # im_shape
(2, 8, 8, 64, 96), # k_shape
( # (9) radial phase encoding (RPE), 3 coils with non-Cartesian sampling along readout
(2, 3, 64, 32, 48), # im_shape
(2, 3, 8, 64, 96), # k_shape
(2, 1, 1, 96), # nkx
(2, 8, 64, 1), # nky
(2, 8, 64, 1), # nkz
(2, 3, 64, 1), # nky
(2, 3, 64, 1), # nkz
'non-uniform', # type_kx
'non-uniform', # type_ky
'non-uniform', # type_kz
'non-uniform', # type_k0
'non-uniform', # type_k1
'non-uniform', # type_k2
),
( # (10) stack of stars, 5 other, 3 coil, oversampling in both FFT and non-uniform directions
(5, 3, 48, 16, 32), # im_shape
(5, 3, 96, 18, 64), # k_shape
(5, 1, 18, 64), # nkx
(5, 1, 18, 64), # nky
(5, 96, 1, 1), # nkz
( # (10) stack of stars, 2 other, 3 coil, oversampling in both FFT and non-uniform directions
(2, 3, 48, 16, 32), # im_shape
(2, 3, 96, 18, 64), # k_shape
(2, 1, 18, 64), # nkx
(2, 1, 18, 64), # nky
(2, 96, 1, 1), # nkz
'non-uniform', # type_kx
'non-uniform', # type_ky
'uniform', # type_kz
'non-uniform', # type_k0
'non-uniform', # type_k1
'uniform', # type_k2
),
( # (11) 2d Cartesian with irregular sampling, oversampling
(1, 1, 1, 96, 128), # im_shape
(1, 1, 1, 1, 192), # k_shape
(1, 1, 1, 192), # nkx
(1, 1, 1, 192), # nky
(1, 1, 1, 1), # nkz
'uniform', # type_kx
'uniform', # type_ky
'zero', # type_kz
'uniform', # type_k0
'zero', # type_k1
'zero', # type_k2
),
],
ids=[
'2d_cartesian_1_coil_no_oversampling',
Expand All @@ -417,10 +430,11 @@ def ismrmrd_cart(ellipse_phantom, tmp_path_factory):
'2d_cartesian_irregular_sampling',
'2d_single_shot_spiral',
'3d_nonuniform_4_coils_2_other',
'2d_nnonuniform_cine_mri_8_cardiac_phases_5_coils',
'2d_cartesian_cine_9_cardiac_phases_6_coils',
'radial_phase_encoding_8_coils_with_oversampling',
'radial_phase_encoding_8_coils_non_cartesian_sampling',
'stack_of_stars_5_other_3_coil_with_oversampling',
'2d_nonuniform_cine_mri_2_cardiac_phases_3_coils',
'2d_cartesian_cine_2_cardiac_phases_3_coils',
'radial_phase_encoding_3_coils_with_oversampling',
'radial_phase_encoding_3_coils_non_cartesian_sampling',
'stack_of_stars_2_other_3_coil_with_oversampling',
'2d_cartesian_irregular_sampling_with_oversampling',
],
)
99 changes: 94 additions & 5 deletions tests/operators/test_fourier_op.py
Original file line number Diff line number Diff line change
Expand Up @@ -3,13 +3,26 @@
import pytest
import torch
from mrpro.data import KData, KTrajectory, SpatialDimension
from mrpro.data.enums import TrajType
from mrpro.data.traj_calculators import KTrajectoryCartesian
from mrpro.operators import FourierOp

from tests import RandomGenerator, dotproduct_adjointness_test
from tests.conftest import COMMON_MR_TRAJECTORIES, create_traj


class NufftTrajektory(KTrajectory):
"""Always returns non-grid trajectory type."""

def _traj_types(
self,
tolerance: float,
) -> tuple[tuple[TrajType, TrajType, TrajType], tuple[TrajType, TrajType, TrajType]]:
true_types = super()._traj_types(tolerance)
modified = tuple([tuple([t & (~TrajType.ONGRID) for t in ts]) for ts in true_types])
return modified


def create_data(im_shape, k_shape, nkx, nky, nkz, type_kx, type_ky, type_kz):
random_generator = RandomGenerator(seed=0)

Expand Down Expand Up @@ -38,16 +51,92 @@
)
fourier_op = FourierOp(recon_matrix=recon_matrix, encoding_matrix=encoding_matrix, traj=trajectory)

# apply forward operator
(kdata,) = fourier_op(img)

# test adjoint property; i.e. <Fu,v> == <u, F^Hv> for all u,v
random_generator = RandomGenerator(seed=0)
u = random_generator.complex64_tensor(size=img.shape)
v = random_generator.complex64_tensor(size=kdata.shape)
u = random_generator.complex64_tensor(size=im_shape)
v = random_generator.complex64_tensor(size=k_shape)
dotproduct_adjointness_test(fourier_op, u, v)

Check failure on line 58 in tests/operators/test_fourier_op.py

View workflow job for this annotation

GitHub Actions / Run Tests and Coverage Report (mrpro_py311)

test_fourier_op_fwd_adj_property[2d_single_shot_spiral] AssertionError

Check failure on line 58 in tests/operators/test_fourier_op.py

View workflow job for this annotation

GitHub Actions / Run Tests and Coverage Report (mrpro_py311)

test_fourier_op_fwd_adj_property[radial_phase_encoding_3_coils_with_oversampling] ValueError: k-space data shape mismatch

Check failure on line 58 in tests/operators/test_fourier_op.py

View workflow job for this annotation

GitHub Actions / Run Tests and Coverage Report (mrpro_py311)

test_fourier_op_fwd_adj_property[radial_phase_encoding_3_coils_non_cartesian_sampling] AssertionError

Check failure on line 58 in tests/operators/test_fourier_op.py

View workflow job for this annotation

GitHub Actions / Run Tests and Coverage Report (mrpro_py312)

test_fourier_op_fwd_adj_property[2d_single_shot_spiral] AssertionError

Check failure on line 58 in tests/operators/test_fourier_op.py

View workflow job for this annotation

GitHub Actions / Run Tests and Coverage Report (mrpro_py312)

test_fourier_op_fwd_adj_property[radial_phase_encoding_3_coils_with_oversampling] ValueError: k-space data shape mismatch

Check failure on line 58 in tests/operators/test_fourier_op.py

View workflow job for this annotation

GitHub Actions / Run Tests and Coverage Report (mrpro_py312)

test_fourier_op_fwd_adj_property[radial_phase_encoding_3_coils_non_cartesian_sampling] AssertionError

Check failure on line 58 in tests/operators/test_fourier_op.py

View workflow job for this annotation

GitHub Actions / Run Tests and Coverage Report (mrpro_py310)

test_fourier_op_fwd_adj_property[2d_single_shot_spiral] AssertionError

Check failure on line 58 in tests/operators/test_fourier_op.py

View workflow job for this annotation

GitHub Actions / Run Tests and Coverage Report (mrpro_py310)

test_fourier_op_fwd_adj_property[radial_phase_encoding_3_coils_with_oversampling] ValueError: k-space data shape mismatch

Check failure on line 58 in tests/operators/test_fourier_op.py

View workflow job for this annotation

GitHub Actions / Run Tests and Coverage Report (mrpro_py310)

test_fourier_op_fwd_adj_property[radial_phase_encoding_3_coils_non_cartesian_sampling] AssertionError


@COMMON_MR_TRAJECTORIES
def test_fourier_op_norm(im_shape, k_shape, nkx, nky, nkz, type_kx, type_ky, type_kz, type_k0, type_k1, type_k2):
"""Test operator norm of Fourier Operator, should be 1."""

# generate random images and k-space trajectories
img, trajectory = create_data(im_shape, k_shape, nkx, nky, nkz, type_kx, type_ky, type_kz)

# create operator
recon_matrix = SpatialDimension(im_shape[-3], im_shape[-2], im_shape[-1])
encoding_matrix = SpatialDimension(
int(trajectory.kz.max() - trajectory.kz.min() + 1),
int(trajectory.ky.max() - trajectory.ky.min() + 1),
int(trajectory.kx.max() - trajectory.kx.min() + 1),
)
fourier_op = FourierOp(recon_matrix=recon_matrix, encoding_matrix=encoding_matrix, traj=trajectory)
(initial_value,) = fourier_op.adjoint(*fourier_op(img))
norm = fourier_op.operator_norm(initial_value, dim=None, max_iterations=4).squeeze()
torch.testing.assert_close(norm, torch.tensor(1.0), atol=0.1, rtol=0.0)

Check failure on line 78 in tests/operators/test_fourier_op.py

View workflow job for this annotation

GitHub Actions / Run Tests and Coverage Report (mrpro_py311)

test_fourier_op_norm[2d_non_cartesian_mri_2_coils] AssertionError: Scalars are not close! Expected 1.0 but got 1.4746354818344116. Absolute difference: 0.4746354818344116 (up to 0.1 allowed) Relative difference: 0.4746354818344116 (up to 0.0 allowed)

Check failure on line 78 in tests/operators/test_fourier_op.py

View workflow job for this annotation

GitHub Actions / Run Tests and Coverage Report (mrpro_py311)

test_fourier_op_norm[2d_single_shot_spiral] AssertionError: Scalars are not close! Expected 1.0 but got 0.7416239380836487. Absolute difference: 0.2583760619163513 (up to 0.1 allowed) Relative difference: 0.2583760619163513 (up to 0.0 allowed)

Check failure on line 78 in tests/operators/test_fourier_op.py

View workflow job for this annotation

GitHub Actions / Run Tests and Coverage Report (mrpro_py311)

test_fourier_op_norm[3d_nonuniform_4_coils_2_other] AssertionError: Scalars are not close! Expected 1.0 but got 0.6982466578483582. Absolute difference: 0.30175334215164185 (up to 0.1 allowed) Relative difference: 0.30175334215164185 (up to 0.0 allowed)

Check failure on line 78 in tests/operators/test_fourier_op.py

View workflow job for this annotation

GitHub Actions / Run Tests and Coverage Report (mrpro_py311)

test_fourier_op_norm[radial_phase_encoding_3_coils_non_cartesian_sampling] AssertionError: Scalars are not close! Expected 1.0 but got 1.54142427444458. Absolute difference: 0.5414242744445801 (up to 0.1 allowed) Relative difference: 0.5414242744445801 (up to 0.0 allowed)

Check failure on line 78 in tests/operators/test_fourier_op.py

View workflow job for this annotation

GitHub Actions / Run Tests and Coverage Report (mrpro_py311)

test_fourier_op_norm[stack_of_stars_2_other_3_coil_with_oversampling] AssertionError: Scalars are not close! Expected 1.0 but got 0.8689930438995361. Absolute difference: 0.13100695610046387 (up to 0.1 allowed) Relative difference: 0.13100695610046387 (up to 0.0 allowed)

Check failure on line 78 in tests/operators/test_fourier_op.py

View workflow job for this annotation

GitHub Actions / Run Tests and Coverage Report (mrpro_py311)

test_fourier_op_norm[2d_cartesian_irregular_sampling_with_oversampling] AssertionError: Scalars are not close! Expected 1.0 but got 0.6897003650665283. Absolute difference: 0.3102996349334717 (up to 0.1 allowed) Relative difference: 0.3102996349334717 (up to 0.0 allowed)

Check failure on line 78 in tests/operators/test_fourier_op.py

View workflow job for this annotation

GitHub Actions / Run Tests and Coverage Report (mrpro_py312)

test_fourier_op_norm[2d_non_cartesian_mri_2_coils] AssertionError: Scalars are not close! Expected 1.0 but got 1.4746354818344116. Absolute difference: 0.4746354818344116 (up to 0.1 allowed) Relative difference: 0.4746354818344116 (up to 0.0 allowed)

Check failure on line 78 in tests/operators/test_fourier_op.py

View workflow job for this annotation

GitHub Actions / Run Tests and Coverage Report (mrpro_py312)

test_fourier_op_norm[2d_single_shot_spiral] AssertionError: Scalars are not close! Expected 1.0 but got 0.7416239380836487. Absolute difference: 0.2583760619163513 (up to 0.1 allowed) Relative difference: 0.2583760619163513 (up to 0.0 allowed)

Check failure on line 78 in tests/operators/test_fourier_op.py

View workflow job for this annotation

GitHub Actions / Run Tests and Coverage Report (mrpro_py312)

test_fourier_op_norm[3d_nonuniform_4_coils_2_other] AssertionError: Scalars are not close! Expected 1.0 but got 0.6982466578483582. Absolute difference: 0.30175334215164185 (up to 0.1 allowed) Relative difference: 0.30175334215164185 (up to 0.0 allowed)

Check failure on line 78 in tests/operators/test_fourier_op.py

View workflow job for this annotation

GitHub Actions / Run Tests and Coverage Report (mrpro_py310)

test_fourier_op_norm[2d_non_cartesian_mri_2_coils] AssertionError: Scalars are not close! Expected 1.0 but got 1.4746354818344116. Absolute difference: 0.4746354818344116 (up to 0.1 allowed) Relative difference: 0.4746354818344116 (up to 0.0 allowed)

Check failure on line 78 in tests/operators/test_fourier_op.py

View workflow job for this annotation

GitHub Actions / Run Tests and Coverage Report (mrpro_py310)

test_fourier_op_norm[2d_single_shot_spiral] AssertionError: Scalars are not close! Expected 1.0 but got 0.7416239380836487. Absolute difference: 0.2583760619163513 (up to 0.1 allowed) Relative difference: 0.2583760619163513 (up to 0.0 allowed)

Check failure on line 78 in tests/operators/test_fourier_op.py

View workflow job for this annotation

GitHub Actions / Run Tests and Coverage Report (mrpro_py310)

test_fourier_op_norm[3d_nonuniform_4_coils_2_other] AssertionError: Scalars are not close! Expected 1.0 but got 0.6982466578483582. Absolute difference: 0.30175334215164185 (up to 0.1 allowed) Relative difference: 0.30175334215164185 (up to 0.0 allowed)


@COMMON_MR_TRAJECTORIES
def test_fourier_op_fft_nufft_forward(
im_shape, k_shape, nkx, nky, nkz, type_kx, type_ky, type_kz, type_k0, type_k1, type_k2
):
"""Test Nufft vs FFT for Fourier operator."""
if not any(t == 'uniform' for t in [type_kx, type_ky, type_kz]):
return # only test for uniform trajectories

img, trajectory = create_data(im_shape, k_shape, nkx, nky, nkz, type_kx, type_ky, type_kz)

recon_matrix = SpatialDimension(im_shape[-3], im_shape[-2], im_shape[-1])
encoding_matrix = SpatialDimension(
int(trajectory.kz.max() - trajectory.kz.min() + 1),
int(trajectory.ky.max() - trajectory.ky.min() + 1),
int(trajectory.kx.max() - trajectory.kx.min() + 1),
)
fourier_op = FourierOp(recon_matrix=recon_matrix, encoding_matrix=encoding_matrix, traj=trajectory)

nufft_fourier_op = FourierOp(
recon_matrix=recon_matrix,
encoding_matrix=encoding_matrix,
traj=NufftTrajektory(trajectory.kz, trajectory.ky, trajectory.kx),
nufft_oversampling=8.0,
)

(result_normal,) = fourier_op(img)
(result_nufft,) = nufft_fourier_op(img)
torch.testing.assert_close(result_normal, result_nufft, atol=1e-4, rtol=5e-3)

Check failure on line 108 in tests/operators/test_fourier_op.py

View workflow job for this annotation

GitHub Actions / Run Tests and Coverage Report (mrpro_py311)

test_fourier_op_fft_nufft_forward[radial_phase_encoding_3_coils_with_oversampling] AssertionError: Tensor-likes are not close! Mismatched elements: 127 / 110592 (0.1%) Greatest absolute difference: 0.000497964269015938 at index (1, 2, 1, 53, 70) (up to 0.0001 allowed) Greatest relative difference: 0.16935928165912628 at index (1, 2, 0, 27, 37) (up to 0.005 allowed)


@COMMON_MR_TRAJECTORIES
def test_fourier_op_fft_nufft_adjoint(
im_shape, k_shape, nkx, nky, nkz, type_kx, type_ky, type_kz, type_k0, type_k1, type_k2
):
"""Test AdjointNufft vs IFFT for Fourier operator."""
if not any(t == 'uniform' for t in [type_kx, type_ky, type_kz]):
return # only test for uniform trajectories
img, trajectory = create_data(im_shape, k_shape, nkx, nky, nkz, type_kx, type_ky, type_kz)
recon_matrix = SpatialDimension(im_shape[-3], im_shape[-2], im_shape[-1])
encoding_matrix = SpatialDimension(
int(trajectory.kz.max() - trajectory.kz.min() + 1),
int(trajectory.ky.max() - trajectory.ky.min() + 1),
int(trajectory.kx.max() - trajectory.kx.min() + 1),
)
fourier_op = FourierOp(recon_matrix=recon_matrix, encoding_matrix=encoding_matrix, traj=trajectory)

nufft_fourier_op = FourierOp(
recon_matrix=recon_matrix,
encoding_matrix=encoding_matrix,
traj=NufftTrajektory(trajectory.kz, trajectory.ky, trajectory.kx),
nufft_oversampling=8.0,
)

k = RandomGenerator(0).complex64_tensor(size=k_shape)
(result_normal,) = fourier_op.H(k)

Check failure on line 135 in tests/operators/test_fourier_op.py

View workflow job for this annotation

GitHub Actions / Run Tests and Coverage Report (mrpro_py310)

test_fourier_op_fft_nufft_adjoint[radial_phase_encoding_3_coils_with_oversampling] ValueError: k-space data shape mismatch
(result_nufft,) = nufft_fourier_op.H(k)
torch.testing.assert_close(result_normal, result_nufft, atol=3e-4, rtol=5e-3)

Check failure on line 137 in tests/operators/test_fourier_op.py

View workflow job for this annotation

GitHub Actions / Run Tests and Coverage Report (mrpro_py310)

test_fourier_op_fft_nufft_adjoint[stack_of_stars_2_other_3_coil_with_oversampling] AssertionError: Tensor-likes are not close! Mismatched elements: 227 / 147456 (0.2%) Greatest absolute difference: 0.0009539038874208927 at index (1, 1, 25, 15, 28) (up to 0.0003 allowed) Greatest relative difference: 0.18000340461730957 at index (1, 1, 19, 6, 29) (up to 0.005 allowed)


@COMMON_MR_TRAJECTORIES
def test_fourier_op_gram(im_shape, k_shape, nkx, nky, nkz, type_kx, type_ky, type_kz, type_k0, type_k1, type_k2):
"""Test gram of Fourier operator."""
Expand All @@ -64,7 +153,7 @@
(expected,) = (fourier_op.H @ fourier_op)(img)
(actual,) = fourier_op.gram(img)

torch.testing.assert_close(actual, expected, rtol=1e-3, atol=1e-3)

Check failure on line 156 in tests/operators/test_fourier_op.py

View workflow job for this annotation

GitHub Actions / Run Tests and Coverage Report (mrpro_py312)

test_fourier_op_gram[2d_non_cartesian_mri_2_coils] AssertionError: Tensor-likes are not close! Mismatched elements: 32768 / 32768 (100.0%) Greatest absolute difference: 0.5523433089256287 at index (0, 0, 0, 3, 35) (up to 0.001 allowed) Greatest relative difference: 0.7502065300941467 at index (0, 1, 0, 28, 116) (up to 0.001 allowed)

Check failure on line 156 in tests/operators/test_fourier_op.py

View workflow job for this annotation

GitHub Actions / Run Tests and Coverage Report (mrpro_py312)

test_fourier_op_gram[2d_single_shot_spiral] AssertionError: Tensor-likes are not close! Mismatched elements: 23945 / 24576 (97.4%) Greatest absolute difference: 0.020737048238515854 at index (0, 0, 0, 21, 80) (up to 0.001 allowed) Greatest relative difference: 0.25013577938079834 at index (0, 1, 0, 71, 127) (up to 0.001 allowed)

Check failure on line 156 in tests/operators/test_fourier_op.py

View workflow job for this annotation

GitHub Actions / Run Tests and Coverage Report (mrpro_py312)

test_fourier_op_gram[radial_phase_encoding_3_coils_with_oversampling] AssertionError: Tensor-likes are not close! Mismatched elements: 589559 / 589824 (100.0%) Greatest absolute difference: 0.17729218304157257 at index (1, 2, 28, 7, 31) (up to 0.001 allowed) Greatest relative difference: 0.5004300475120544 at index (1, 0, 63, 13, 36) (up to 0.001 allowed)

Check failure on line 156 in tests/operators/test_fourier_op.py

View workflow job for this annotation

GitHub Actions / Run Tests and Coverage Report (mrpro_py312)

test_fourier_op_gram[radial_phase_encoding_3_coils_non_cartesian_sampling] AssertionError: Tensor-likes are not close! Mismatched elements: 589632 / 589824 (100.0%) Greatest absolute difference: 0.2131546139717102 at index (1, 2, 28, 7, 31) (up to 0.001 allowed) Greatest relative difference: 0.5005220770835876 at index (0, 0, 62, 2, 29) (up to 0.001 allowed)

Check failure on line 156 in tests/operators/test_fourier_op.py

View workflow job for this annotation

GitHub Actions / Run Tests and Coverage Report (mrpro_py310)

test_fourier_op_gram[2d_non_cartesian_mri_2_coils] AssertionError: Tensor-likes are not close! Mismatched elements: 32768 / 32768 (100.0%) Greatest absolute difference: 0.5523433089256287 at index (0, 0, 0, 3, 35) (up to 0.001 allowed) Greatest relative difference: 0.7502065300941467 at index (0, 1, 0, 28, 116) (up to 0.001 allowed)

Check failure on line 156 in tests/operators/test_fourier_op.py

View workflow job for this annotation

GitHub Actions / Run Tests and Coverage Report (mrpro_py310)

test_fourier_op_gram[2d_single_shot_spiral] AssertionError: Tensor-likes are not close! Mismatched elements: 23945 / 24576 (97.4%) Greatest absolute difference: 0.020737048238515854 at index (0, 0, 0, 21, 80) (up to 0.001 allowed) Greatest relative difference: 0.25013577938079834 at index (0, 1, 0, 71, 127) (up to 0.001 allowed)


@pytest.mark.parametrize(
Expand Down
Loading