From bc2e05a3048ee9b880d00cdf6a73973d3f344fa8 Mon Sep 17 00:00:00 2001 From: Christoph Kolbitsch Date: Sat, 4 Jan 2025 15:50:59 +0100 Subject: [PATCH 1/5] Make tests faster (#589) Co-authored-by: Felix Zimmermann --- tests/algorithms/dcf/test_dcf_voronoi.py | 16 ++-- tests/algorithms/test_optimizers.py | 9 +- tests/conftest.py | 100 +++++++++++++---------- tests/data/_IsmrmrdRawTestData.py | 12 +-- tests/data/conftest.py | 2 +- tests/data/test_dcf_data.py | 6 +- tests/data/test_kdata.py | 16 ++-- tests/operators/models/conftest.py | 2 +- tests/operators/test_operator_norm.py | 6 +- 9 files changed, 91 insertions(+), 78 deletions(-) diff --git a/tests/algorithms/dcf/test_dcf_voronoi.py b/tests/algorithms/dcf/test_dcf_voronoi.py index d88e7d2c0..208c94278 100644 --- a/tests/algorithms/dcf/test_dcf_voronoi.py +++ b/tests/algorithms/dcf/test_dcf_voronoi.py @@ -25,11 +25,11 @@ def example_traj_rad_2d(n_kr, n_ka, phi0=0.0, broadcast=True): @pytest.mark.parametrize( ('n_kr', 'n_ka', 'phi0', 'broadcast'), [ - (100, 20, 0, True), - (100, 1, 0, True), - (100, 20, torch.pi / 4, True), - (100, 1, torch.pi / 4, True), - (100, 1, 0, False), + (20, 20, 0, True), + (20, 1, 0, True), + (20, 20, torch.pi / 4, True), + (20, 1, torch.pi / 4, True), + (20, 1, 0, False), ], ) def test_dcf_rad_traj_voronoi(n_kr, n_ka, phi0, broadcast): @@ -56,7 +56,7 @@ def test_dcf_rad_traj_voronoi(n_kr, n_ka, phi0, broadcast): assert dcf.shape == traj.broadcasted_shape, 'DCF shape should match broadcasted trajectory shape' -@pytest.mark.parametrize(('n_k2', 'n_k1', 'n_k0'), [(40, 16, 20), (1, 2, 2)]) +@pytest.mark.parametrize(('n_k2', 'n_k1', 'n_k0'), [(4, 6, 8), (1, 2, 2)]) def test_dcf_3d_cart_traj_broadcast_voronoi(n_k2, n_k1, n_k0): """Compare voronoi-based dcf calculation for broadcasted 3D regular Cartesian trajectory to analytical solution which is 1 for each k-space @@ -76,7 +76,7 @@ def test_dcf_3d_cart_traj_broadcast_voronoi(n_k2, n_k1, n_k0): torch.testing.assert_close(dcf[:, 1:-1, 1:-1, 1:-1], dcf_analytical[:, 1:-1, 1:-1, 1:-1]) -@pytest.mark.parametrize(('n_k2', 'n_k1', 'n_k0'), [(40, 16, 20), (1, 2, 2)]) +@pytest.mark.parametrize(('n_k2', 'n_k1', 'n_k0'), [(4, 6, 8), (1, 2, 2)]) def test_dcf_3d_cart_full_traj_voronoi(n_k2, n_k1, n_k0): """Compare voronoi-based dcf calculation for full 3D regular Cartesian trajectory to analytical solution which is 1 for each k-space point.""" @@ -103,7 +103,7 @@ def test_dcf_3d_cart_full_traj_voronoi(n_k2, n_k1, n_k0): @pytest.mark.parametrize( ('n_k2', 'n_k1', 'n_k0', 'k2_steps', 'k1_steps', 'k0_steps'), - [(30, 20, 10, (1.0, 0.5, 0.25), (1.0, 0.5), (1.0,))], + [(4, 6, 8, (1.0, 0.5, 0.25), (1.0, 0.5), (1.0,))], ) def test_dcf_3d_cart_nonuniform_traj_voronoi(n_k2, n_k1, n_k0, k2_steps, k1_steps, k0_steps): """Compare voronoi-based dcf calculation for 3D nonuniform Cartesian diff --git a/tests/algorithms/test_optimizers.py b/tests/algorithms/test_optimizers.py index bda6820e1..ddc355841 100644 --- a/tests/algorithms/test_optimizers.py +++ b/tests/algorithms/test_optimizers.py @@ -9,7 +9,8 @@ @pytest.mark.parametrize('enforce_bounds_on_x1', [True, False]) @pytest.mark.parametrize( - ('optimizer', 'optimizer_kwargs'), [(adam, {'lr': 0.02, 'max_iter': 10000}), (lbfgs, {'lr': 1.0})] + ('optimizer', 'optimizer_kwargs'), + [(adam, {'lr': 0.02, 'max_iter': 2000, 'betas': (0.8, 0.999)}), (lbfgs, {'lr': 1.0, 'max_iter': 20})], ) def test_optimizers_rosenbrock(optimizer, enforce_bounds_on_x1, optimizer_kwargs): # use Rosenbrock function as test case with 2D test data @@ -17,8 +18,8 @@ def test_optimizers_rosenbrock(optimizer, enforce_bounds_on_x1, optimizer_kwargs rosen_brock = Rosenbrock(a, b) # initial point of optimization - x1 = torch.tensor([a / 3.14]) - x2 = torch.tensor([3.14]) + x1 = torch.tensor([a / 1.23]) + x2 = torch.tensor([1.23]) x1.grad = torch.tensor([2.78]) x2.grad = torch.tensor([-1.0]) params_init = [x1, x2] @@ -45,7 +46,7 @@ def test_optimizers_rosenbrock(optimizer, enforce_bounds_on_x1, optimizer_kwargs params_result = constrain_op(*params_result) # obtained solution should match analytical - torch.testing.assert_close(torch.tensor(params_result), analytical_solution) + torch.testing.assert_close(torch.tensor(params_result), analytical_solution, rtol=1e-4, atol=0) for p, before, grad_before in zip(params_init, params_init_before, params_init_grad_before, strict=True): assert p == before, 'the initial parameter should not have changed during optimization' diff --git a/tests/conftest.py b/tests/conftest.py index 30ae9c229..477b93be4 100644 --- a/tests/conftest.py +++ b/tests/conftest.py @@ -192,7 +192,7 @@ def random_acq_info(random_acquisition): return acq_info -@pytest.fixture(params=({'seed': 0, 'n_other': 10, 'n_k2': 40, 'n_k1': 20},)) +@pytest.fixture(params=({'seed': 0, 'n_other': 4, 'n_k2': 12, 'n_k1': 14},)) def random_kheader_shape(request, random_acquisition, random_full_ismrmrd_header): """Random (not necessarily valid) KHeader with defined shape.""" # Get dimensions @@ -263,6 +263,20 @@ def ismrmrd_cart(ellipse_phantom, tmp_path_factory): return ismrmrd_kdata +@pytest.fixture(scope='session') +def ismrmrd_cart_high_res(ellipse_phantom, tmp_path_factory): + """Fully sampled cartesian data set.""" + ismrmrd_filename = tmp_path_factory.mktemp('mrpro') / 'ismrmrd_cart_high_res.h5' + ismrmrd_kdata = IsmrmrdRawTestData( + filename=ismrmrd_filename, + matrix_size=256, + noise_level=0.0, + repetitions=3, + phantom=ellipse_phantom.phantom, + ) + return ismrmrd_kdata + + COMMON_MR_TRAJECTORIES = pytest.mark.parametrize( ('im_shape', 'k_shape', 'nkx', 'nky', 'nkz', 'type_kx', 'type_ky', 'type_kz', 'type_k0', 'type_k1', 'type_k2'), [ @@ -331,12 +345,12 @@ def ismrmrd_cart(ellipse_phantom, tmp_path_factory): 'zero', # type_k1 'zero', # type_k2 ), - ( # (5) 3d non-uniform, 4 coils, 2 other - (2, 4, 16, 32, 64), # im_shape - (2, 4, 16, 32, 64), # k_shape - (2, 16, 32, 64), # nkx - (2, 16, 32, 64), # nky - (2, 16, 32, 64), # nkz + ( # (5) 3d non-uniform, 3 coils, 2 other + (2, 3, 10, 12, 14), # im_shape + (2, 3, 6, 8, 10), # k_shape + (2, 6, 8, 10), # nkx + (2, 6, 8, 10), # nky + (2, 6, 8, 10), # nkz 'non-uniform', # type_kx 'non-uniform', # type_ky 'non-uniform', # type_kz @@ -344,12 +358,12 @@ 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 3 cardiac phases, 2 coils + (3, 2, 1, 64, 64), # im_shape + (3, 2, 1, 18, 128), # k_shape + (3, 1, 18, 128), # nkx + (3, 1, 18, 128), # nky + (3, 1, 1, 1), # nkz 'non-uniform', # type_kx 'non-uniform', # type_ky 'zero', # type_kz @@ -357,12 +371,12 @@ def ismrmrd_cart(ellipse_phantom, tmp_path_factory): '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 @@ -370,12 +384,12 @@ def ismrmrd_cart(ellipse_phantom, tmp_path_factory): '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 - (2, 1, 1, 96), # nkx - (2, 8, 64, 1), # nky - (2, 8, 64, 1), # nkz + ( # (8) radial phase encoding (RPE), 3 coils, with oversampling in both FFT and non-uniform directions + (2, 3, 12, 14, 16), # im_shape + (2, 3, 8, 10, 12), # k_shape + (2, 1, 1, 12), # nkx + (2, 8, 10, 1), # nky + (2, 8, 10, 1), # nkz 'uniform', # type_kx 'non-uniform', # type_ky 'non-uniform', # type_kz @@ -383,12 +397,12 @@ def ismrmrd_cart(ellipse_phantom, tmp_path_factory): '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 - (2, 1, 1, 96), # nkx - (2, 8, 64, 1), # nky - (2, 8, 64, 1), # nkz + ( # (9) radial phase encoding (RPE), 2 coils with non-Cartesian sampling along readout + (2, 2, 12, 14, 16), # im_shape + (2, 2, 8, 10, 12), # k_shape + (2, 2, 1, 12), # nkx + (2, 2, 10, 1), # nky + (2, 2, 10, 1), # nkz 'non-uniform', # type_kx 'non-uniform', # type_ky 'non-uniform', # type_kz @@ -396,12 +410,12 @@ def ismrmrd_cart(ellipse_phantom, tmp_path_factory): '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, 3 other, 2 coil, oversampling in both FFT and non-uniform directions + (3, 2, 24, 16, 32), # im_shape + (3, 2, 48, 12, 14), # k_shape + (3, 1, 12, 14), # nkx + (3, 1, 12, 14), # nky + (3, 48, 1, 1), # nkz 'non-uniform', # type_kx 'non-uniform', # type_ky 'uniform', # type_kz @@ -416,11 +430,11 @@ def ismrmrd_cart(ellipse_phantom, tmp_path_factory): '2d_non_cartesian_mri_2_coils', '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', + '3d_nonuniform_3_coils_2_other', + '2d_nonuniform_cine_mri_3_cardiac_phases_2_coils', + '2d_cartesian_cine_2_cardiac_phases_3_coils', + 'radial_phase_encoding_3_coils_with_oversampling', + 'radial_phase_encoding_2_coils_non_cartesian_sampling', + 'stack_of_stars_3_other_2_coil_with_oversampling', ], ) diff --git a/tests/data/_IsmrmrdRawTestData.py b/tests/data/_IsmrmrdRawTestData.py index efeff6ed1..89e5ac901 100644 --- a/tests/data/_IsmrmrdRawTestData.py +++ b/tests/data/_IsmrmrdRawTestData.py @@ -59,8 +59,8 @@ class IsmrmrdRawTestData: def __init__( self, filename: str | Path, - matrix_size: int = 256, - n_coils: int = 8, + matrix_size: int = 64, + n_coils: int = 4, oversampling: int = 2, repetitions: int = 1, flag_invalid_reps: bool = False, @@ -167,12 +167,12 @@ def __init__( # Encoded and recon spaces encoding_fov = ismrmrd.xsd.fieldOfViewMm() - encoding_fov.x = self.oversampling * 256 - encoding_fov.y = 256 + encoding_fov.x = self.oversampling * matrix_size + encoding_fov.y = matrix_size encoding_fov.z = 5 recon_fov = ismrmrd.xsd.fieldOfViewMm() - recon_fov.x = 256 - recon_fov.y = 256 + recon_fov.x = matrix_size + recon_fov.y = matrix_size recon_fov.z = 5 encoding_matrix = ismrmrd.xsd.matrixSizeType() diff --git a/tests/data/conftest.py b/tests/data/conftest.py index bd37b28c2..1d92b14a5 100644 --- a/tests/data/conftest.py +++ b/tests/data/conftest.py @@ -50,7 +50,7 @@ def random_mandatory_ismrmrd_header(request) -> xsd.ismrmrdschema.ismrmrdHeader: return xsd.ismrmrdschema.ismrmrdHeader(encoding=[encoding], experimentalConditions=experimental_conditions) -@pytest.fixture(params=({'seed': 0, 'n_other': 2, 'n_coils': 16, 'n_z': 32, 'n_y': 128, 'n_x': 256},)) +@pytest.fixture(params=({'seed': 0, 'n_other': 2, 'n_coils': 8, 'n_z': 16, 'n_y': 32, 'n_x': 64},)) def random_test_data(request): seed, n_other, n_coils, n_z, n_y, n_x = ( request.param['seed'], diff --git a/tests/data/test_dcf_data.py b/tests/data/test_dcf_data.py index 314ad2ba2..7df7ea256 100644 --- a/tests/data/test_dcf_data.py +++ b/tests/data/test_dcf_data.py @@ -53,11 +53,9 @@ def test_dcf_spiral_traj_voronoi(n_kr, n_ki, n_ka): def test_dcf_spiral_traj_voronoi_singlespiral(): """For three z-stacked spirals in the x,y plane, the center spiral should be the same as a single 2D spiral. - - Issue #84 """ - n_kr = 100 # points along each spiral ar - n_ki = 5 # turns per spiral arm spirals nka spiral arms + n_kr = 30 # points along each spiral arm + n_ki = 5 # turns per spiral arm trajectory_single = example_traj_spiral_2d(n_kr, n_ki, 1) # A new trajectroy with three spirals stacked in z direction. diff --git a/tests/data/test_kdata.py b/tests/data/test_kdata.py index fa3e4ebd9..42cd061f6 100644 --- a/tests/data/test_kdata.py +++ b/tests/data/test_kdata.py @@ -181,16 +181,16 @@ def test_KData_calibration_lines(ismrmrd_cart_with_calibration_lines): assert kdata.data.shape[-2] == ismrmrd_cart_with_calibration_lines.n_separate_calibration_lines -def test_KData_kspace(ismrmrd_cart): +def test_KData_kspace(ismrmrd_cart_high_res): """Read in data and verify k-space by comparing reconstructed image.""" - kdata = KData.from_file(ismrmrd_cart.filename, DummyTrajectory()) + kdata = KData.from_file(ismrmrd_cart_high_res.filename, DummyTrajectory()) ff_op = FastFourierOp(dim=(-1, -2)) (reconstructed_img,) = ff_op.adjoint(kdata.data) # Due to discretisation artifacts the reconstructed image will be different to the reference image. Using standard # testing functions such as numpy.testing.assert_almost_equal fails because there are few voxels with high # differences along the edges of the elliptic objects. - assert relative_image_difference(reconstructed_img[0, 0, 0, ...], ismrmrd_cart.img_ref) <= 0.05 + assert relative_image_difference(reconstructed_img[0, 0, 0, ...], ismrmrd_cart_high_res.img_ref) <= 0.05 @pytest.mark.parametrize(('field', 'value'), [('lamor_frequency_proton', 42.88 * 1e6), ('tr', torch.tensor([24.3]))]) @@ -423,8 +423,8 @@ def test_KData_split_k2_into_other(consistently_shaped_kdata, monkeypatch, n_oth ('subset_label', 'subset_idx'), [ ('repetition', torch.tensor([1], dtype=torch.int32)), - ('average', torch.tensor([3, 4, 5], dtype=torch.int32)), - ('phase', torch.tensor([2, 2, 8], dtype=torch.int32)), + ('average', torch.tensor([1, 2, 3], dtype=torch.int32)), + ('phase', torch.tensor([2, 2, 3], dtype=torch.int32)), ], ) def test_KData_select_other_subset(consistently_shaped_kdata, monkeypatch, subset_label, subset_idx): @@ -520,9 +520,9 @@ def test_modify_acq_info(random_kheader_shape): assert kheader.acq_info.position.z.shape == (n_other, n_k2, n_k1, 1) -def test_KData_compress_coils(ismrmrd_cart): +def test_KData_compress_coils(ismrmrd_cart_high_res): """Test coil combination does not alter image content (much).""" - kdata = KData.from_file(ismrmrd_cart.filename, DummyTrajectory()) + kdata = KData.from_file(ismrmrd_cart_high_res.filename, DummyTrajectory()) kdata = kdata.compress_coils(n_compressed_coils=4) ff_op = FastFourierOp(dim=(-1, -2)) (reconstructed_img,) = ff_op.adjoint(kdata.data) @@ -530,7 +530,7 @@ def test_KData_compress_coils(ismrmrd_cart): # Image content of each coil is the same. Therefore we only compare one coil image but we need to normalize. reconstructed_img = reconstructed_img[0, 0, 0, ...].abs() reconstructed_img /= reconstructed_img.max() - ref_img = ismrmrd_cart.img_ref[0, 0, 0, ...].abs() + ref_img = ismrmrd_cart_high_res.img_ref[0, 0, 0, ...].abs() ref_img /= ref_img.max() assert relative_image_difference(reconstructed_img, ref_img) <= 0.1 diff --git a/tests/operators/models/conftest.py b/tests/operators/models/conftest.py index 4aab81ae0..7ebc58ad4 100644 --- a/tests/operators/models/conftest.py +++ b/tests/operators/models/conftest.py @@ -47,7 +47,7 @@ def create_parameter_tensor_tuples( - parameter_shape=(10, 5, 100, 100, 100), number_of_tensors=2 + parameter_shape=(10, 5, 12, 14, 16), number_of_tensors=2 ) -> tuple[torch.Tensor, ...]: """Create tuples of tensors as input to operators.""" random_generator = RandomGenerator(seed=0) diff --git a/tests/operators/test_operator_norm.py b/tests/operators/test_operator_norm.py index c6e13dcf9..387254384 100644 --- a/tests/operators/test_operator_norm.py +++ b/tests/operators/test_operator_norm.py @@ -117,17 +117,17 @@ def test_finite_difference_operator_norm(dim): finite_difference_operator = FiniteDifferenceOp(dim=dim, mode='forward') # initialize random image of appropriate shape depending on the dimensionality - image_shape = (1, *tuple([16 for _ in range(len(dim))])) + image_shape = (1, *([8] * len(dim))) random_image = random_generator.complex64_tensor(size=image_shape) # calculate the operator norm - finite_difference_operator_norm = finite_difference_operator.operator_norm(random_image, dim=dim, max_iterations=64) + finite_difference_operator_norm = finite_difference_operator.operator_norm(random_image, dim=dim, max_iterations=32) # closed form solution of the operator norm finite_difference_operator_norm_true = sqrt(len(dim) * 4) torch.testing.assert_close( - finite_difference_operator_norm.item(), finite_difference_operator_norm_true, atol=1e-2, rtol=1e-2 + finite_difference_operator_norm.item(), finite_difference_operator_norm_true, atol=0.1, rtol=0 ) From 4e3caf68018604a0020d92d39c6708b3cb576c5d Mon Sep 17 00:00:00 2001 From: Lunin Leonid Date: Mon, 6 Jan 2025 13:21:57 +0100 Subject: [PATCH 2/5] Remove the binder infrastructure (#538) --- .dockerignore | 11 +++++------ README.md | 2 +- binder/environment.yml | 9 --------- docs/source/contributor_guide.rst | 19 +++++++------------ docs/source/examples.rst | 4 +++- docs/source/user_guide.rst | 8 ++++++-- 6 files changed, 22 insertions(+), 31 deletions(-) delete mode 100755 binder/environment.yml diff --git a/.dockerignore b/.dockerignore index 0981abb00..0e896c7e5 100644 --- a/.dockerignore +++ b/.dockerignore @@ -1,6 +1,5 @@ -.github -.vscode -binder -examples -docs -tests \ No newline at end of file +.github +.vscode +examples +docs +tests diff --git a/README.md b/README.md index 49d7e2c3d..8817da3ae 100644 --- a/README.md +++ b/README.md @@ -11,7 +11,7 @@ MR image reconstruction and processing package specifically developed for PyTorc - **Source code:** - **Documentation:** - **Bug reports:** -- **Try it out:** [![Binder](https://mybinder.org/badge_logo.svg)](https://mybinder.org/v2/gh/PTB-MR/mrpro.git/main?labpath=examples) +- **Try it out:** [![Open In Colab](https://colab.research.google.com/assets/colab-badge.svg)](https://colab.research.google.com/github/PTB-MR/mrpro) ## Awards diff --git a/binder/environment.yml b/binder/environment.yml deleted file mode 100755 index 55c8d56c0..000000000 --- a/binder/environment.yml +++ /dev/null @@ -1,9 +0,0 @@ -name: Binder -channels: - - conda-forge -dependencies: - - python=3.11 - - pytorch::pytorch - - pytorch::cpuonly - - pip: - - mrpro[notebook] @ git+https://github.com/PTB-MR/mrpro diff --git a/docs/source/contributor_guide.rst b/docs/source/contributor_guide.rst index 28bdd4797..cc1458883 100644 --- a/docs/source/contributor_guide.rst +++ b/docs/source/contributor_guide.rst @@ -5,8 +5,6 @@ Contributor Guide Repo structure ============== This repository uses a *pyproject.toml* file to specify all the requirements. -If you need a "normal" *requirements.txt* file, please have a look in *binder*. There you find a *requirements.txt* -automatically created from *pyproject.toml* using GitHub actions. **.github/workflows** Definitions of GitHub action workflows to carry out formatting checks, run tests and automatically create this @@ -18,12 +16,9 @@ automatically created from *pyproject.toml* using GitHub actions. **.docs** Files to create this documentation. -**binder** - *environment.yml* to install MRpro when starting a binder session. - **examples** Python scripts showcasing how MRpro can be used. Any data needed has to be available from - an online repository (e.g. zenodo) such that it can be automatically downloaded. The scripts + an online repository (e.g. zenodo) such that it can be automatically downloaded. The scripts are automatically translated to jupyter notebooks using GitHub actions. Individual cells should be indicated with ``# %%``. For markdown cells use ``# %% [markdown]``. The translation from python script to jupyter notebook is done using @@ -67,14 +62,14 @@ src/mrpro structure Linting ======= -We use Ruff for linting. If you are using VSCode, you can install an -`extension `_, -which we have also added to the list of extensions that VSCode should recommend when you open the code. +We use Ruff for linting. If you are using VSCode, you can install an +`extension `_, +which we have also added to the list of extensions that VSCode should recommend when you open the code. We also run `mypy `_ as a type checker. -In CI, our linting is driven by `pre-commit `_. +In CI, our linting is driven by `pre-commit `_. If you install MRpro via ``pip install -e .[test]``, pre-commit will be installed in your python environment. -You can either add pre-commit to your git pre-commit hooks, requiring it to pass before each commit (``pre-commit install``), +You can either add pre-commit to your git pre-commit hooks, requiring it to pass before each commit (``pre-commit install``), or run it manually using ``pre-commit run --all-files`` after making your changes, before requesting a PR review. Naming convention @@ -115,4 +110,4 @@ We are still in pre-release mode and do not guarantee a stable API / strict sema Compatibility ============= We aim to always be compatible with the latest stable pytorch release and the latest python version supported by pytorch. We are compatible with one previous python version. -Our type hints will usually only be valid with the latest pytorch version. \ No newline at end of file +Our type hints will usually only be valid with the latest pytorch version. diff --git a/docs/source/examples.rst b/docs/source/examples.rst index 95c654781..8d3937d35 100644 --- a/docs/source/examples.rst +++ b/docs/source/examples.rst @@ -1,8 +1,10 @@ Examples ======== +.. |colab-badge| image:: https://colab.research.google.com/assets/colab-badge.svg + :target: https://colab.research.google.com/github/PTB-MR/mrpro Notebooks with examples of how you can use MRpro. -All of the notebooks can directly be run via binder or colab from the repo. +Each notebook can be launched in Colab |colab-badge| .. toctree:: :maxdepth: 1 diff --git a/docs/source/user_guide.rst b/docs/source/user_guide.rst index c3f6bd2f1..1c8b81a55 100644 --- a/docs/source/user_guide.rst +++ b/docs/source/user_guide.rst @@ -6,7 +6,7 @@ MRpro is a MR image reconstruction and processing framework specifically develop The data classes utilize torch tensors for storing data such as MR raw data or reconstructed image data. Where possible batch parallelisation of pytorch is utilized to speed up image reconstruction. -MRpro is designed to work directly from MR raw data using the `MRD `_ data format. +MRpro is designed to work directly from MR raw data using the `MRD `_ data format. A basic pipeline would contain the following steps: @@ -15,9 +15,13 @@ A basic pipeline would contain the following steps: * Data reconstruction * Image processing +.. |colab-badge| image:: https://colab.research.google.com/assets/colab-badge.svg + :target: https://colab.research.google.com/github/PTB-MR/mrpro + The following provides some basic information about these steps. For more detailed information please have a look at the notebooks in the *examples* folder. -You can easily start a binder session via the badge in the *README* and give the notebooks a try without having to + +You can easily launch notebooks via the |colab-badge| badge and give the notebooks a try without having to install anything. Reading in raw data From a400469eb0af1e7fd8295b72281315c5df1e9fb1 Mon Sep 17 00:00:00 2001 From: Patrick Schuenke <37338697+schuenke@users.noreply.github.com> Date: Mon, 6 Jan 2025 14:04:59 +0100 Subject: [PATCH 3/5] Add zenodo reference (#587) --- CITATION.cff | 55 ++++++++++++++++++++++++++++++++++++++++++++++++++++ README.md | 1 + 2 files changed, 56 insertions(+) create mode 100644 CITATION.cff diff --git a/CITATION.cff b/CITATION.cff new file mode 100644 index 000000000..372724793 --- /dev/null +++ b/CITATION.cff @@ -0,0 +1,55 @@ +abstract: +

MRpro is a MR image reconstruction and processing framework specifically + developed to work well with PyTorch. The data classes utilize torch tensors for + storing data such as MR raw data or reconstructed image data. Where possible batch + parallelisation of PyTorch is utilized to speed up image reconstruction.

+authors: + - affiliation: Physikalisch-Technische Bundesanstalt + family-names: Zimmermann + given-names: Felix Frederik + orcid: 0000-0002-0862-8973 + - affiliation: Physikalisch-Technische Bundesanstalt + family-names: Schuenke + given-names: Patrick + orcid: 0000-0002-3179-4830 + - affiliation: Physikalisch-Technische Bundesanstalt + family-names: Brahma + given-names: Sherine + orcid: 0000-0003-4340-6513 + - affiliation: Physikalisch-Technische Bundesanstalt + family-names: Guastini + given-names: Mara + - affiliation: Physikalisch-Technische Bundesanstalt + family-names: Hammacher + given-names: Johannes + - affiliation: Physikalisch-Technische Bundesanstalt + family-names: Kofler + given-names: Andreas + orcid: 0000-0002-7416-4433 + - affiliation: Physikalisch-Technische Bundesanstalt + family-names: Kranich Redshaw + given-names: Catarina + - affiliation: Physikalisch-Technische Bundesanstalt + family-names: Lunin + given-names: Leonid + orcid: 0000-0001-6469-5532 + - affiliation: Physikalisch-Technische Bundesanstalt + family-names: Martin + given-names: Stefan + - affiliation: Physikalisch-Technische Bundesanstalt + family-names: Schote + given-names: David + orcid: 0000-0003-3468-0676 + - affiliation: Physikalisch-Technische Bundesanstalt + family-names: Kolbitsch + given-names: Christoph + orcid: 0000-0002-4355-8368 +cff-version: 1.2.0 +date-released: "2024-12-17" +doi: 10.5281/zenodo.14509598 +license: + - apache-2.0 +repository-code: https://github.com/PTB-MR/mrpro/ +title: MRpro - PyTorch-based MR image reconstruction and processing package +type: software +url: https://doi.org/10.5281/zenodo.14509598 diff --git a/README.md b/README.md index 8817da3ae..4dcce2430 100644 --- a/README.md +++ b/README.md @@ -5,6 +5,7 @@ ![Python](https://img.shields.io/badge/python-3.10%20%7C%203.11%20%7C%203.12-blue) [![License](https://img.shields.io/badge/License-Apache%202.0-blue.svg)](https://opensource.org/licenses/Apache-2.0) ![Coverage Bagde](https://img.shields.io/endpoint?url=https://gist.githubusercontent.com/ckolbPTB/48e334a10caf60e6708d7c712e56d241/raw/coverage.json) +[![DOI](https://zenodo.org/badge/DOI/10.5281/zenodo.14509598.svg)](https://doi.org/10.5281/zenodo.14509598) MR image reconstruction and processing package specifically developed for PyTorch. From 075f1c0c377bbebcf259554f83a0c5c9368bf55e Mon Sep 17 00:00:00 2001 From: Christoph Kolbitsch Date: Mon, 6 Jan 2025 18:03:42 +0100 Subject: [PATCH 4/5] Example for Basics and Cartesian reconstruction (#562) Co-authored-by: Patrick Schuenke --- examples/cartesian_reconstruction.ipynb | 713 ++++++++++++++++++++++++ examples/cartesian_reconstruction.py | 387 +++++++++++++ 2 files changed, 1100 insertions(+) create mode 100644 examples/cartesian_reconstruction.ipynb create mode 100644 examples/cartesian_reconstruction.py diff --git a/examples/cartesian_reconstruction.ipynb b/examples/cartesian_reconstruction.ipynb new file mode 100644 index 000000000..cf5b4c95b --- /dev/null +++ b/examples/cartesian_reconstruction.ipynb @@ -0,0 +1,713 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "1606223d", + "metadata": {}, + "source": [ + "# Basics of MRpro and Cartesian Reconstructions\n", + "Here, we are going to have a look at a few basics of MRpro and reconstruct data acquired with a Cartesian sampling\n", + "pattern." + ] + }, + { + "cell_type": "markdown", + "id": "6442c71d", + "metadata": {}, + "source": [ + "## Overview\n", + "\n", + "In this notebook, we are going to explore the MRpro KData object and the included header parameters. We will then use\n", + "a FFT-operator in order to reconstruct data acquired with a Cartesian sampling scheme. We will also reconstruct data\n", + "acquired on a Cartesian grid but with partial echo and partial Fourier acceleration. Finally, we will reconstruct a\n", + "Cartesian scan with regular undersampling using iterative SENSE." + ] + }, + { + "cell_type": "markdown", + "id": "c8e96446", + "metadata": {}, + "source": [ + "## Import MRpro and download data" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "ac2a1011", + "metadata": {}, + "outputs": [], + "source": [ + "# Get the raw data from zenodo\n", + "import tempfile\n", + "from pathlib import Path\n", + "\n", + "import zenodo_get\n", + "\n", + "data_folder = Path(tempfile.mkdtemp())\n", + "dataset = '14173489'\n", + "zenodo_get.zenodo_get([dataset, '-r', 5, '-o', data_folder]) # r: retries" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "e5089e32", + "metadata": {}, + "outputs": [], + "source": [ + "# List the downloaded files\n", + "for f in data_folder.iterdir():\n", + " print(f.name)" + ] + }, + { + "cell_type": "markdown", + "id": "c217f4f9", + "metadata": { + "lines_to_next_cell": 2 + }, + "source": [ + "We have three different scans obtained from the same object with the same FOV and resolution:\n", + "\n", + "- cart_t1.mrd is a fully sampled Cartesian acquisition\n", + "\n", + "- cart_t1_msense_integrated.mrd is accelerated using regular undersampling and self-calibrated SENSE\n", + "\n", + "- cart_t1_partial_echo_partial_fourier.mrd is accelerated using partial echo and partial Fourier" + ] + }, + { + "cell_type": "markdown", + "id": "c1d2f8e1", + "metadata": {}, + "source": [ + "## Read in raw data and explore header\n", + "\n", + "To read in an ISMRMRD raw data file (*.mrd), we can simply pass on the file name to a `KData` object.\n", + "Additionally, we need to provide information about the trajectory. In MRpro, this is done using trajectory\n", + "calculators. These are functions that calculate the trajectory based on the acquisition information and additional\n", + "parameters provided to the calculators (e.g. the angular step for a radial acquisition).\n", + "\n", + "In this case, we have a Cartesian acquisition. This means that we only need to provide a Cartesian trajectory\n", + "calculator (called `KTrajectoryCartesian` in MRpro) without any further parameters." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "24b9f069", + "metadata": {}, + "outputs": [], + "source": [ + "from mrpro.data import KData\n", + "from mrpro.data.traj_calculators import KTrajectoryCartesian\n", + "\n", + "kdata = KData.from_file(data_folder / 'cart_t1.mrd', KTrajectoryCartesian())" + ] + }, + { + "cell_type": "markdown", + "id": "adda70f1", + "metadata": {}, + "source": [ + "Now we can explore this data object." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "703daa5a", + "metadata": {}, + "outputs": [], + "source": [ + "# Start with simply calling print(kdata), whichs gives us a nice overview of the KData object.\n", + "print(kdata)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "960a2c8a", + "metadata": {}, + "outputs": [], + "source": [ + "# We can also have a look at more specific header information like the 1H Lamor frequency\n", + "print(kdata.header.lamor_frequency_proton)" + ] + }, + { + "cell_type": "markdown", + "id": "b083edfc", + "metadata": {}, + "source": [ + "## Reconstruction of fully sampled acquisition\n", + "\n", + "For the reconstruction of a fully sampled Cartesian acquisition, we can use a simple Fast Fourier Transform (FFT).\n", + "\n", + "Let's create an FFT-operator (called `FastFourierOp` in MRpro) and apply it to our `KData` object. Please note that\n", + "all MRpro operators currently only work on PyTorch tensors and not on the MRpro objects directly. Therefore, we have\n", + "to call the operator on kdata.data. One other important feature of MRpro operators is that they always return a\n", + "tuple of PyTorch tensors, even if the output is only a single tensor. This is why we use the `(img,)` syntax below." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "24e79e3a", + "metadata": {}, + "outputs": [], + "source": [ + "from mrpro.operators import FastFourierOp\n", + "\n", + "fft_op = FastFourierOp(dim=(-2, -1))\n", + "(img,) = fft_op.adjoint(kdata.data)" + ] + }, + { + "cell_type": "markdown", + "id": "9eebf952", + "metadata": {}, + "source": [ + "Let's have a look at the shape of the obtained tensor." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "51158c3a", + "metadata": {}, + "outputs": [], + "source": [ + "print(img.shape)" + ] + }, + { + "cell_type": "markdown", + "id": "2ea1ce99", + "metadata": {}, + "source": [ + "We can see that the second dimension, which is the coil dimension, is 16. This means we still have a coil resolved\n", + "dataset (i.e. one image for each coil element). We can use a simply root-sum-of-squares approach to combine them into\n", + "one. Later, we will do something a bit more sophisticated. We can also see that the x-dimension is 512. This is\n", + "because in MRI we commonly oversample the readout direction by a factor 2 leading to a FOV twice as large as we\n", + "actually need. We can either remove this oversampling along the readout direction or we can simply tell the\n", + "`FastFourierOp` to crop the image by providing the correct output matrix size (recon_matrix)." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "59686881", + "metadata": {}, + "outputs": [], + "source": [ + "# Create FFT-operator with correct output matrix size\n", + "fft_op = FastFourierOp(\n", + " dim=(-2, -1),\n", + " recon_matrix=kdata.header.recon_matrix,\n", + " encoding_matrix=kdata.header.encoding_matrix,\n", + ")\n", + "\n", + "(img,) = fft_op.adjoint(kdata.data)\n", + "print(img.shape)" + ] + }, + { + "cell_type": "markdown", + "id": "9aebc2d1", + "metadata": {}, + "source": [ + "Now, we have an image which is 256 x 256 voxel as we would expect. Let's combine the data from the different receiver\n", + "coils using root-sum-of-squares and then display the image. Note that we usually index from behind in MRpro\n", + "(i.e. -1 for the last, -4 for the fourth last (coil) dimension) to allow for more than one 'other' dimension." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "367e3ea7", + "metadata": {}, + "outputs": [], + "source": [ + "import matplotlib.pyplot as plt\n", + "import torch\n", + "\n", + "# Combine data from different coils\n", + "img_fully_sampled = torch.sqrt(torch.sum(img**2, dim=-4)).abs().squeeze()\n", + "\n", + "# plot the image\n", + "plt.imshow(img_fully_sampled)" + ] + }, + { + "cell_type": "markdown", + "id": "93f1bc4f", + "metadata": {}, + "source": [ + "Great! That was very easy! Let's try to reconstruct the next dataset." + ] + }, + { + "cell_type": "markdown", + "id": "229754fd", + "metadata": {}, + "source": [ + "## Reconstruction of acquisition with partial echo and partial Fourier" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "8d8f961d", + "metadata": { + "lines_to_next_cell": 2 + }, + "outputs": [], + "source": [ + "# Read in the data\n", + "kdata_pe_pf = KData.from_file(data_folder / 'cart_t1_partial_echo_partial_fourier.mrd', KTrajectoryCartesian())\n", + "\n", + "# Create FFT-operator with correct output matrix size\n", + "fft_op = FastFourierOp(\n", + " dim=(-2, -1),\n", + " recon_matrix=kdata.header.recon_matrix,\n", + " encoding_matrix=kdata.header.encoding_matrix,\n", + ")\n", + "\n", + "# Reconstruct coil resolved image(s)\n", + "(img_pe_pf,) = fft_op.adjoint(kdata_pe_pf.data)\n", + "\n", + "# Combine data from different coils using root-sum-of-squares\n", + "img_pe_pf = torch.sqrt(torch.sum(img_pe_pf**2, dim=-4)).abs().squeeze()\n", + "\n", + "# Plot both images\n", + "fig, ax = plt.subplots(1, 2, squeeze=False)\n", + "ax[0, 0].imshow(img_fully_sampled)\n", + "ax[0, 1].imshow(img_pe_pf)" + ] + }, + { + "cell_type": "markdown", + "id": "9f2eaaec", + "metadata": { + "lines_to_next_cell": 2 + }, + "source": [ + "Well, we got an image, but when we compare it to the previous result, it seems like the head has shrunk.\n", + "Since that's extremely unlikely, there's probably a mistake in our reconstruction.\n", + "\n", + "Let's step back and check out the trajectories for both scans." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "90a6a2bb", + "metadata": {}, + "outputs": [], + "source": [ + "print(kdata.traj)" + ] + }, + { + "cell_type": "markdown", + "id": "bfd8a051", + "metadata": { + "lines_to_next_cell": 2 + }, + "source": [ + "We see that the trajectory has kz, ky, and kx components. Kx and ky only vary along one dimension.\n", + "This is because MRpro saves the trajectory in the most efficient way.\n", + "To get the full trajectory as a tensor, we can just call as_tensor()." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "9fd011dd", + "metadata": {}, + "outputs": [], + "source": [ + "# Plot the fully sampled trajectory (in blue)\n", + "plt.plot(kdata.traj.as_tensor()[2, 0, 0, :, :].flatten(), kdata.traj.as_tensor()[1, 0, 0, :, :].flatten(), 'ob')\n", + "\n", + "# Plot the partial echo and partial Fourier trajectory (in red)\n", + "plt.plot(\n", + " kdata_pe_pf.traj.as_tensor()[2, 0, 0, :, :].flatten(), kdata_pe_pf.traj.as_tensor()[1, 0, 0, :, :].flatten(), '+r'\n", + ")" + ] + }, + { + "cell_type": "markdown", + "id": "721db0d8", + "metadata": {}, + "source": [ + "We see that for the fully sampled acquisition, the k-space is covered symmetrically from -256 to 255 along the\n", + "readout direction and from -128 to 127 along the phase encoding direction. For the acquisition with partial Fourier\n", + "and partial echo acceleration, this is of course not the case and the k-space is asymmetrical.\n", + "\n", + "Our FFT-operator does not know about this and simply assumes that the acquisition is symmetric and any difference\n", + "between encoding and recon matrix needs to be zero-padded symmetrically.\n", + "\n", + "To take the asymmetric acquisition into account and sort the data correctly into a matrix where we can apply the\n", + "FFT-operator to, we have got the `CartesianSamplingOp` in MRpro. This operator calculates a sorting index based on the\n", + "k-space trajectory and the dimensions of the encoding k-space.\n", + "\n", + "Let's try it out!" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "2746809f", + "metadata": {}, + "outputs": [], + "source": [ + "from mrpro.operators import CartesianSamplingOp\n", + "\n", + "cart_sampling_op = CartesianSamplingOp(encoding_matrix=kdata_pe_pf.header.encoding_matrix, traj=kdata_pe_pf.traj)" + ] + }, + { + "cell_type": "markdown", + "id": "f0994a1a", + "metadata": {}, + "source": [ + "Now, we first apply the CartesianSamplingOp and then call the FFT-operator." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "7ab5132e", + "metadata": {}, + "outputs": [], + "source": [ + "(img_pe_pf,) = fft_op.adjoint(cart_sampling_op.adjoint(kdata_pe_pf.data)[0])\n", + "img_pe_pf = torch.sqrt(torch.sum(img_pe_pf**2, dim=-4)).abs().squeeze()\n", + "\n", + "fig, ax = plt.subplots(1, 2, squeeze=False)\n", + "ax[0, 0].imshow(img_fully_sampled)\n", + "ax[0, 1].imshow(img_pe_pf)" + ] + }, + { + "cell_type": "markdown", + "id": "1506fd56", + "metadata": { + "lines_to_next_cell": 0 + }, + "source": [] + }, + { + "cell_type": "markdown", + "id": "0effea6e", + "metadata": { + "lines_to_next_cell": 2 + }, + "source": [ + "Voila! We've got the same brains, and they're the same size!\n", + "\n", + "But wait a second—something still looks a bit off. In the bottom left corner, it seems like there's a \"hole\"\n", + "in the brain. That definitely shouldn't be there.\n", + "\n", + "The issue is that we combined the data from the different coils using a root-sum-of-squares approach.\n", + "While it's simple, it's not the ideal method. Typically, coil sensitivity maps are calculated to combine the data\n", + "from different coils. In MRpro, you can do this by calculating coil sensitivity data and then creating a\n", + "`SensitivityOp` to combine the data after image reconstruction." + ] + }, + { + "cell_type": "markdown", + "id": "7dd4d555", + "metadata": {}, + "source": [ + "We have different options for calculating coil sensitivity maps from the image data of the various coils.\n", + "Here, we're going to use the Walsh method." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "860fd313", + "metadata": {}, + "outputs": [], + "source": [ + "from mrpro.algorithms.csm import walsh\n", + "from mrpro.operators import SensitivityOp\n", + "\n", + "# Calculate coil sensitivity maps\n", + "(img_pe_pf,) = fft_op.adjoint(cart_sampling_op.adjoint(kdata_pe_pf.data)[0])\n", + "\n", + "# This algorithms is designed to calculate coil sensitivity maps for each other dimension.\n", + "csm_data = walsh(img_pe_pf[0, ...], smoothing_width=5)[None, ...]\n", + "\n", + "# Create SensitivityOp\n", + "csm_op = SensitivityOp(csm_data)\n", + "\n", + "# Reconstruct coil-combined image\n", + "(img_pe_pf,) = csm_op.adjoint(fft_op.adjoint(cart_sampling_op.adjoint(kdata_pe_pf.data)[0])[0])\n", + "img_pe_pf = img_pe_pf.abs().squeeze()\n", + "\n", + "fig, ax = plt.subplots(1, 2, squeeze=False)\n", + "ax[0, 0].imshow(img_fully_sampled)\n", + "ax[0, 1].imshow(img_pe_pf.squeeze())" + ] + }, + { + "cell_type": "markdown", + "id": "772c3843", + "metadata": {}, + "source": [ + "Tada! The \"hole\" is gone, and the image looks much better.\n", + "\n", + "When we reconstructed the image, we called the adjoint method of several different operators one after the other. That\n", + "was a bit cumbersome. To make our life easier, MRpro allows to combine the operators first and then call the adjoint\n", + "of the composite operator. We have to keep in mind that we have to put them in the order of the forward method of the\n", + "operators. By calling the adjoint, the order will be automatically reversed." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "4653f66a", + "metadata": {}, + "outputs": [], + "source": [ + "# Create composite operator\n", + "acq_op = cart_sampling_op @ fft_op @ csm_op\n", + "(img_pe_pf,) = acq_op.adjoint(kdata_pe_pf.data)\n", + "img_pe_pf = img_pe_pf.abs().squeeze()\n", + "\n", + "fig, ax = plt.subplots(1, 2, squeeze=False)\n", + "ax[0, 0].imshow(img_fully_sampled)\n", + "ax[0, 1].imshow(img_pe_pf)" + ] + }, + { + "cell_type": "markdown", + "id": "76892ff3", + "metadata": {}, + "source": [ + "Although we now have got a nice looking image, it was still a bit cumbersome to create it. We had to define several\n", + "different operators and chain them together. Wouldn't it be nice if this could be done automatically?\n", + "\n", + "That is why we also included some top-level reconstruction algorithms in MRpro. For this whole steps from above,\n", + "we can simply call a `DirectReconstruction`. A `DirectReconstruction` object can be created from only the information\n", + "in the `KData` object.\n", + "\n", + "In contrast to operators, top-level reconstruction algorithms operate on the data objects of MRpro, i.e. the input is\n", + "a `KData` object and the output is an image data (called `IData` in MRpro) object. To get the tensor content of the\n", + "`IData` object, we can call its `rss` method." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "64550e58", + "metadata": {}, + "outputs": [], + "source": [ + "from mrpro.algorithms.reconstruction import DirectReconstruction\n", + "\n", + "# Create DirectReconstruction object from KData object\n", + "direct_recon_pe_pf = DirectReconstruction(kdata_pe_pf)\n", + "\n", + "# Reconstruct image by calling the DirectReconstruction object\n", + "idat_pe_pf = direct_recon_pe_pf(kdata_pe_pf)\n", + "\n", + "fig, ax = plt.subplots(1, 2, squeeze=False)\n", + "ax[0, 0].imshow(img_fully_sampled)\n", + "ax[0, 1].imshow(idat_pe_pf.rss().squeeze())" + ] + }, + { + "cell_type": "markdown", + "id": "95ec919c", + "metadata": { + "lines_to_next_cell": 2 + }, + "source": [ + "This is much simpler — everything happens in the background, so we don't have to worry about it.\n", + "Let's try it on the undersampled dataset now." + ] + }, + { + "cell_type": "markdown", + "id": "4944ac9d", + "metadata": {}, + "source": [ + "## Reconstruction of undersampled data" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "bac530d3", + "metadata": {}, + "outputs": [], + "source": [ + "kdata_us = KData.from_file(data_folder / 'cart_t1_msense_integrated.mrd', KTrajectoryCartesian())\n", + "direct_recon_us = DirectReconstruction(kdata_us)\n", + "idat_us = direct_recon_us(kdata_us)\n", + "\n", + "fig, ax = plt.subplots(1, 2, squeeze=False)\n", + "ax[0, 0].imshow(img_fully_sampled)\n", + "ax[0, 1].imshow(idat_us.rss().squeeze())" + ] + }, + { + "cell_type": "markdown", + "id": "7e4ec563", + "metadata": {}, + "source": [ + "As expected, we can see undersampling artifacts in the image. In order to get rid of them, we can use an iterative\n", + "SENSE algorithm. As you might have guessed, this is also included in MRpro.\n", + "\n", + "Similarly to the `DirectReconstruction`, we can create an `IterativeSENSEReconstruction` and apply it to the\n", + "undersampled data.\n", + "\n", + "One important thing to keep in mind is that this only works if the coil maps that we use do not have any\n", + "undersampling artifacts. Commonly, we would get them from a fully sampled self-calibration reference lines in the\n", + "center of k-space or a separate coil sensitivity scan.\n", + "\n", + "As a first step, we are going to assume that we have got a nice fully sampled reference scan like our partial echo and\n", + "partial Fourier acquisition. We can get the `CsmData`, which is needed for the `IterativeSENSEReconstruction`, from\n", + "the previous reconstruction." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "aad1a1bf", + "metadata": {}, + "outputs": [], + "source": [ + "from mrpro.algorithms.reconstruction import IterativeSENSEReconstruction\n", + "\n", + "it_sense_recon = IterativeSENSEReconstruction(kdata=kdata_us, csm=direct_recon_pe_pf.csm)\n", + "idat_us = it_sense_recon(kdata_us)\n", + "\n", + "fig, ax = plt.subplots(1, 2, squeeze=False)\n", + "ax[0, 0].imshow(img_fully_sampled)\n", + "ax[0, 1].imshow(idat_us.rss().squeeze())" + ] + }, + { + "cell_type": "markdown", + "id": "21237288", + "metadata": {}, + "source": [ + "That worked well, but in practice, we don't want to acquire a fully sampled version of our scan just to\n", + "reconstruct it. A more efficient approach is to get a few self-calibration lines in the center of k-space\n", + "to create a low-resolution, fully sampled image.\n", + "\n", + "In our scan, these lines are part of the dataset, but they aren't used for image reconstruction since\n", + "they're only meant for calibration (i.e., coil sensitivity map calculation). Because they're not labeled\n", + "for imaging, MRpro ignores them by default when reading the data. However, we can set a flag when calling\n", + "`from_file` to read in just those lines for reconstructing the coil sensitivity maps." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "a4d39855", + "metadata": {}, + "outputs": [], + "source": [ + "from mrpro.data.acq_filters import is_coil_calibration_acquisition\n", + "\n", + "kdata_calib_lines = KData.from_file(\n", + " data_folder / 'cart_t1_msense_integrated.mrd',\n", + " KTrajectoryCartesian(),\n", + " acquisition_filter_criterion=lambda acq: is_coil_calibration_acquisition(acq),\n", + ")\n", + "\n", + "direct_recon_calib_lines = DirectReconstruction(kdata_calib_lines)\n", + "im_calib_lines = direct_recon_calib_lines(kdata_calib_lines)\n", + "\n", + "plt.imshow(im_calib_lines.rss().squeeze())" + ] + }, + { + "cell_type": "markdown", + "id": "ea089177", + "metadata": {}, + "source": [ + "Although this only yields a low-resolution image, it is good enough to calculate coil sensitivity maps." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "90da630f", + "metadata": {}, + "outputs": [], + "source": [ + "# Visualize coil sensitivity maps of all 16 coils\n", + "assert direct_recon_calib_lines.csm is not None # needed for type checking\n", + "fig, ax = plt.subplots(4, 4, squeeze=False)\n", + "for idx, cax in enumerate(ax.flatten()):\n", + " cax.imshow(direct_recon_calib_lines.csm.data[0, idx, 0, ...].abs())" + ] + }, + { + "cell_type": "markdown", + "id": "7377ad90", + "metadata": {}, + "source": [ + "Now, we can use these coil sensitivity maps to reconstruct our SENSE scan." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "f6ded207", + "metadata": {}, + "outputs": [], + "source": [ + "it_sense_recon = IterativeSENSEReconstruction(kdata_us, csm=direct_recon_calib_lines.csm)\n", + "idat_us = it_sense_recon(kdata_us)\n", + "\n", + "fig, ax = plt.subplots(1, 2, squeeze=False)\n", + "ax[0, 0].imshow(img_fully_sampled)\n", + "ax[0, 1].imshow(idat_us.rss().squeeze())" + ] + }, + { + "cell_type": "markdown", + "id": "211b9a6f", + "metadata": { + "lines_to_next_cell": 0 + }, + "source": [] + }, + { + "cell_type": "markdown", + "id": "f17eab4e", + "metadata": {}, + "source": [ + "The final image is a little worse (nothing beats fully sampled high-resolution scans for coil map\n", + "calculation), but we've managed to get rid of the undersampling artifacts inside the brain. If you want to\n", + "further improve the coil sensitivity map quality, try:\n", + "- using different methods to calculate them, e.g. `mrpro.algorithms.csm.inati`\n", + "- playing around with the parameters of these methods\n", + "- applying a smoothing filter on the images (or ideally directly in k-space) used to calculate the coil\n", + " sensitivity maps" + ] + } + ], + "metadata": { + "jupytext": { + "cell_metadata_filter": "-all" + }, + "kernelspec": { + "display_name": "Python 3 (ipykernel)", + "language": "python", + "name": "python3" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/examples/cartesian_reconstruction.py b/examples/cartesian_reconstruction.py new file mode 100644 index 000000000..e6069fc9b --- /dev/null +++ b/examples/cartesian_reconstruction.py @@ -0,0 +1,387 @@ +# %% [markdown] +# # Basics of MRpro and Cartesian Reconstructions +# Here, we are going to have a look at a few basics of MRpro and reconstruct data acquired with a Cartesian sampling +# pattern. + +# %% [markdown] +# ## Overview +# +# In this notebook, we are going to explore the MRpro KData object and the included header parameters. We will then use +# a FFT-operator in order to reconstruct data acquired with a Cartesian sampling scheme. We will also reconstruct data +# acquired on a Cartesian grid but with partial echo and partial Fourier acceleration. Finally, we will reconstruct a +# Cartesian scan with regular undersampling using iterative SENSE. + +# %% [markdown] +# ## Import MRpro and download data + +# %% +# Get the raw data from zenodo +import tempfile +from pathlib import Path + +import zenodo_get + +data_folder = Path(tempfile.mkdtemp()) +dataset = '14173489' +zenodo_get.zenodo_get([dataset, '-r', 5, '-o', data_folder]) # r: retries + +# %% +# List the downloaded files +for f in data_folder.iterdir(): + print(f.name) + +# %% [markdown] +# We have three different scans obtained from the same object with the same FOV and resolution: +# +# - cart_t1.mrd is a fully sampled Cartesian acquisition +# +# - cart_t1_msense_integrated.mrd is accelerated using regular undersampling and self-calibrated SENSE +# +# - cart_t1_partial_echo_partial_fourier.mrd is accelerated using partial echo and partial Fourier + + +# %% [markdown] +# ## Read in raw data and explore header +# +# To read in an ISMRMRD raw data file (*.mrd), we can simply pass on the file name to a `KData` object. +# Additionally, we need to provide information about the trajectory. In MRpro, this is done using trajectory +# calculators. These are functions that calculate the trajectory based on the acquisition information and additional +# parameters provided to the calculators (e.g. the angular step for a radial acquisition). +# +# In this case, we have a Cartesian acquisition. This means that we only need to provide a Cartesian trajectory +# calculator (called `KTrajectoryCartesian` in MRpro) without any further parameters. + +# %% +from mrpro.data import KData +from mrpro.data.traj_calculators import KTrajectoryCartesian + +kdata = KData.from_file(data_folder / 'cart_t1.mrd', KTrajectoryCartesian()) + +# %% [markdown] +# Now we can explore this data object. + +# %% +# Start with simply calling print(kdata), whichs gives us a nice overview of the KData object. +print(kdata) + +# %% +# We can also have a look at more specific header information like the 1H Lamor frequency +print(kdata.header.lamor_frequency_proton) + +# %% [markdown] +# ## Reconstruction of fully sampled acquisition +# +# For the reconstruction of a fully sampled Cartesian acquisition, we can use a simple Fast Fourier Transform (FFT). +# +# Let's create an FFT-operator (called `FastFourierOp` in MRpro) and apply it to our `KData` object. Please note that +# all MRpro operators currently only work on PyTorch tensors and not on the MRpro objects directly. Therefore, we have +# to call the operator on kdata.data. One other important feature of MRpro operators is that they always return a +# tuple of PyTorch tensors, even if the output is only a single tensor. This is why we use the `(img,)` syntax below. + +# %% +from mrpro.operators import FastFourierOp + +fft_op = FastFourierOp(dim=(-2, -1)) +(img,) = fft_op.adjoint(kdata.data) + +# %% [markdown] +# Let's have a look at the shape of the obtained tensor. + +# %% +print(img.shape) + +# %% [markdown] +# We can see that the second dimension, which is the coil dimension, is 16. This means we still have a coil resolved +# dataset (i.e. one image for each coil element). We can use a simply root-sum-of-squares approach to combine them into +# one. Later, we will do something a bit more sophisticated. We can also see that the x-dimension is 512. This is +# because in MRI we commonly oversample the readout direction by a factor 2 leading to a FOV twice as large as we +# actually need. We can either remove this oversampling along the readout direction or we can simply tell the +# `FastFourierOp` to crop the image by providing the correct output matrix size (recon_matrix). + +# %% +# Create FFT-operator with correct output matrix size +fft_op = FastFourierOp( + dim=(-2, -1), + recon_matrix=kdata.header.recon_matrix, + encoding_matrix=kdata.header.encoding_matrix, +) + +(img,) = fft_op.adjoint(kdata.data) +print(img.shape) + +# %% [markdown] +# Now, we have an image which is 256 x 256 voxel as we would expect. Let's combine the data from the different receiver +# coils using root-sum-of-squares and then display the image. Note that we usually index from behind in MRpro +# (i.e. -1 for the last, -4 for the fourth last (coil) dimension) to allow for more than one 'other' dimension. + +# %% +import matplotlib.pyplot as plt +import torch + +# Combine data from different coils +img_fully_sampled = torch.sqrt(torch.sum(img**2, dim=-4)).abs().squeeze() + +# plot the image +plt.imshow(img_fully_sampled) + +# %% [markdown] +# Great! That was very easy! Let's try to reconstruct the next dataset. + +# %% [markdown] +# ## Reconstruction of acquisition with partial echo and partial Fourier + +# %% +# Read in the data +kdata_pe_pf = KData.from_file(data_folder / 'cart_t1_partial_echo_partial_fourier.mrd', KTrajectoryCartesian()) + +# Create FFT-operator with correct output matrix size +fft_op = FastFourierOp( + dim=(-2, -1), + recon_matrix=kdata.header.recon_matrix, + encoding_matrix=kdata.header.encoding_matrix, +) + +# Reconstruct coil resolved image(s) +(img_pe_pf,) = fft_op.adjoint(kdata_pe_pf.data) + +# Combine data from different coils using root-sum-of-squares +img_pe_pf = torch.sqrt(torch.sum(img_pe_pf**2, dim=-4)).abs().squeeze() + +# Plot both images +fig, ax = plt.subplots(1, 2, squeeze=False) +ax[0, 0].imshow(img_fully_sampled) +ax[0, 1].imshow(img_pe_pf) + + +# %% [markdown] +# Well, we got an image, but when we compare it to the previous result, it seems like the head has shrunk. +# Since that's extremely unlikely, there's probably a mistake in our reconstruction. +# +# Let's step back and check out the trajectories for both scans. + + +# %% +print(kdata.traj) + +# %% [markdown] +# We see that the trajectory has kz, ky, and kx components. Kx and ky only vary along one dimension. +# This is because MRpro saves the trajectory in the most efficient way. +# To get the full trajectory as a tensor, we can just call as_tensor(). + + +# %% +# Plot the fully sampled trajectory (in blue) +plt.plot(kdata.traj.as_tensor()[2, 0, 0, :, :].flatten(), kdata.traj.as_tensor()[1, 0, 0, :, :].flatten(), 'ob') + +# Plot the partial echo and partial Fourier trajectory (in red) +plt.plot( + kdata_pe_pf.traj.as_tensor()[2, 0, 0, :, :].flatten(), kdata_pe_pf.traj.as_tensor()[1, 0, 0, :, :].flatten(), '+r' +) + +# %% [markdown] +# We see that for the fully sampled acquisition, the k-space is covered symmetrically from -256 to 255 along the +# readout direction and from -128 to 127 along the phase encoding direction. For the acquisition with partial Fourier +# and partial echo acceleration, this is of course not the case and the k-space is asymmetrical. +# +# Our FFT-operator does not know about this and simply assumes that the acquisition is symmetric and any difference +# between encoding and recon matrix needs to be zero-padded symmetrically. +# +# To take the asymmetric acquisition into account and sort the data correctly into a matrix where we can apply the +# FFT-operator to, we have got the `CartesianSamplingOp` in MRpro. This operator calculates a sorting index based on the +# k-space trajectory and the dimensions of the encoding k-space. +# +# Let's try it out! + +# %% +from mrpro.operators import CartesianSamplingOp + +cart_sampling_op = CartesianSamplingOp(encoding_matrix=kdata_pe_pf.header.encoding_matrix, traj=kdata_pe_pf.traj) + +# %% [markdown] +# Now, we first apply the CartesianSamplingOp and then call the FFT-operator. + +# %% +(img_pe_pf,) = fft_op.adjoint(cart_sampling_op.adjoint(kdata_pe_pf.data)[0]) +img_pe_pf = torch.sqrt(torch.sum(img_pe_pf**2, dim=-4)).abs().squeeze() + +fig, ax = plt.subplots(1, 2, squeeze=False) +ax[0, 0].imshow(img_fully_sampled) +ax[0, 1].imshow(img_pe_pf) + +# %% [markdown] +# %% [markdown] +# Voila! We've got the same brains, and they're the same size! +# +# But wait a second—something still looks a bit off. In the bottom left corner, it seems like there's a "hole" +# in the brain. That definitely shouldn't be there. +# +# The issue is that we combined the data from the different coils using a root-sum-of-squares approach. +# While it's simple, it's not the ideal method. Typically, coil sensitivity maps are calculated to combine the data +# from different coils. In MRpro, you can do this by calculating coil sensitivity data and then creating a +# `SensitivityOp` to combine the data after image reconstruction. + + +# %% [markdown] +# We have different options for calculating coil sensitivity maps from the image data of the various coils. +# Here, we're going to use the Walsh method. + +# %% +from mrpro.algorithms.csm import walsh +from mrpro.operators import SensitivityOp + +# Calculate coil sensitivity maps +(img_pe_pf,) = fft_op.adjoint(cart_sampling_op.adjoint(kdata_pe_pf.data)[0]) + +# This algorithms is designed to calculate coil sensitivity maps for each other dimension. +csm_data = walsh(img_pe_pf[0, ...], smoothing_width=5)[None, ...] + +# Create SensitivityOp +csm_op = SensitivityOp(csm_data) + +# Reconstruct coil-combined image +(img_pe_pf,) = csm_op.adjoint(fft_op.adjoint(cart_sampling_op.adjoint(kdata_pe_pf.data)[0])[0]) +img_pe_pf = img_pe_pf.abs().squeeze() + +fig, ax = plt.subplots(1, 2, squeeze=False) +ax[0, 0].imshow(img_fully_sampled) +ax[0, 1].imshow(img_pe_pf.squeeze()) + +# %% [markdown] +# Tada! The "hole" is gone, and the image looks much better. +# +# When we reconstructed the image, we called the adjoint method of several different operators one after the other. That +# was a bit cumbersome. To make our life easier, MRpro allows to combine the operators first and then call the adjoint +# of the composite operator. We have to keep in mind that we have to put them in the order of the forward method of the +# operators. By calling the adjoint, the order will be automatically reversed. + +# %% +# Create composite operator +acq_op = cart_sampling_op @ fft_op @ csm_op +(img_pe_pf,) = acq_op.adjoint(kdata_pe_pf.data) +img_pe_pf = img_pe_pf.abs().squeeze() + +fig, ax = plt.subplots(1, 2, squeeze=False) +ax[0, 0].imshow(img_fully_sampled) +ax[0, 1].imshow(img_pe_pf) + +# %% [markdown] +# Although we now have got a nice looking image, it was still a bit cumbersome to create it. We had to define several +# different operators and chain them together. Wouldn't it be nice if this could be done automatically? +# +# That is why we also included some top-level reconstruction algorithms in MRpro. For this whole steps from above, +# we can simply call a `DirectReconstruction`. A `DirectReconstruction` object can be created from only the information +# in the `KData` object. +# +# In contrast to operators, top-level reconstruction algorithms operate on the data objects of MRpro, i.e. the input is +# a `KData` object and the output is an image data (called `IData` in MRpro) object. To get the tensor content of the +# `IData` object, we can call its `rss` method. + +# %% +from mrpro.algorithms.reconstruction import DirectReconstruction + +# Create DirectReconstruction object from KData object +direct_recon_pe_pf = DirectReconstruction(kdata_pe_pf) + +# Reconstruct image by calling the DirectReconstruction object +idat_pe_pf = direct_recon_pe_pf(kdata_pe_pf) + +fig, ax = plt.subplots(1, 2, squeeze=False) +ax[0, 0].imshow(img_fully_sampled) +ax[0, 1].imshow(idat_pe_pf.rss().squeeze()) + +# %% [markdown] +# This is much simpler — everything happens in the background, so we don't have to worry about it. +# Let's try it on the undersampled dataset now. + + +# %% [markdown] +# ## Reconstruction of undersampled data + +# %% +kdata_us = KData.from_file(data_folder / 'cart_t1_msense_integrated.mrd', KTrajectoryCartesian()) +direct_recon_us = DirectReconstruction(kdata_us) +idat_us = direct_recon_us(kdata_us) + +fig, ax = plt.subplots(1, 2, squeeze=False) +ax[0, 0].imshow(img_fully_sampled) +ax[0, 1].imshow(idat_us.rss().squeeze()) + +# %% [markdown] +# As expected, we can see undersampling artifacts in the image. In order to get rid of them, we can use an iterative +# SENSE algorithm. As you might have guessed, this is also included in MRpro. + +# Similarly to the `DirectReconstruction`, we can create an `IterativeSENSEReconstruction` and apply it to the +# undersampled data. +# +# One important thing to keep in mind is that this only works if the coil maps that we use do not have any +# undersampling artifacts. Commonly, we would get them from a fully sampled self-calibration reference lines in the +# center of k-space or a separate coil sensitivity scan. +# +# As a first step, we are going to assume that we have got a nice fully sampled reference scan like our partial echo and +# partial Fourier acquisition. We can get the `CsmData`, which is needed for the `IterativeSENSEReconstruction`, from +# the previous reconstruction. + +# %% +from mrpro.algorithms.reconstruction import IterativeSENSEReconstruction + +it_sense_recon = IterativeSENSEReconstruction(kdata=kdata_us, csm=direct_recon_pe_pf.csm) +idat_us = it_sense_recon(kdata_us) + +fig, ax = plt.subplots(1, 2, squeeze=False) +ax[0, 0].imshow(img_fully_sampled) +ax[0, 1].imshow(idat_us.rss().squeeze()) + +# %% [markdown] +# That worked well, but in practice, we don't want to acquire a fully sampled version of our scan just to +# reconstruct it. A more efficient approach is to get a few self-calibration lines in the center of k-space +# to create a low-resolution, fully sampled image. +# +# In our scan, these lines are part of the dataset, but they aren't used for image reconstruction since +# they're only meant for calibration (i.e., coil sensitivity map calculation). Because they're not labeled +# for imaging, MRpro ignores them by default when reading the data. However, we can set a flag when calling +# `from_file` to read in just those lines for reconstructing the coil sensitivity maps. + +# %% +from mrpro.data.acq_filters import is_coil_calibration_acquisition + +kdata_calib_lines = KData.from_file( + data_folder / 'cart_t1_msense_integrated.mrd', + KTrajectoryCartesian(), + acquisition_filter_criterion=lambda acq: is_coil_calibration_acquisition(acq), +) + +direct_recon_calib_lines = DirectReconstruction(kdata_calib_lines) +im_calib_lines = direct_recon_calib_lines(kdata_calib_lines) + +plt.imshow(im_calib_lines.rss().squeeze()) + +# %% [markdown] +# Although this only yields a low-resolution image, it is good enough to calculate coil sensitivity maps. + +# %% +# Visualize coil sensitivity maps of all 16 coils +assert direct_recon_calib_lines.csm is not None # needed for type checking +fig, ax = plt.subplots(4, 4, squeeze=False) +for idx, cax in enumerate(ax.flatten()): + cax.imshow(direct_recon_calib_lines.csm.data[0, idx, 0, ...].abs()) + +# %% [markdown] +# Now, we can use these coil sensitivity maps to reconstruct our SENSE scan. + +# %% +it_sense_recon = IterativeSENSEReconstruction(kdata_us, csm=direct_recon_calib_lines.csm) +idat_us = it_sense_recon(kdata_us) + +fig, ax = plt.subplots(1, 2, squeeze=False) +ax[0, 0].imshow(img_fully_sampled) +ax[0, 1].imshow(idat_us.rss().squeeze()) + +# %% [markdown] +# %% [markdown] +# The final image is a little worse (nothing beats fully sampled high-resolution scans for coil map +# calculation), but we've managed to get rid of the undersampling artifacts inside the brain. If you want to +# further improve the coil sensitivity map quality, try: +# - using different methods to calculate them, e.g. `mrpro.algorithms.csm.inati` +# - playing around with the parameters of these methods +# - applying a smoothing filter on the images (or ideally directly in k-space) used to calculate the coil +# sensitivity maps From 3a4ed2a3f94911615ba928b3d6f3c8d06e5c3db3 Mon Sep 17 00:00:00 2001 From: Lunin Leonid Date: Mon, 6 Jan 2025 23:23:39 +0100 Subject: [PATCH 5/5] Add download badge in docs build (#588) Co-authored-by: Patrick Schuenke <37338697+schuenke@users.noreply.github.com> --- .github/workflows/docs.yml | 18 ++++++++++++++++++ 1 file changed, 18 insertions(+) diff --git a/.github/workflows/docs.yml b/.github/workflows/docs.yml index fae4bf3ba..dc683049d 100644 --- a/.github/workflows/docs.yml +++ b/.github/workflows/docs.yml @@ -106,6 +106,24 @@ jobs: run: | echo "current jupyter-notebook: ${{ matrix.notebook }}" + - name: Add nb-myst download badge + run: | + notebook=${{ matrix.notebook }} + notebook_name=$(basename $notebook) + download_badge_md="[![Download notebook](https://img.shields.io/badge/Download-notebook-blue?logo=jupyter)](path:$notebook_name)" + python_command="import nbformat as nbf\n\ + nb = nbf.read(open('$notebook'), as_version=4)\n\ + # if the 1st cell is md and has colab text => add space after\n\ + if nb['cells'][0]['cell_type'] == 'markdown' and 'colab' in nb['cells'][0]['source'].lower():\n\ + nb['cells'][0]['source'] += ' '\n\ + # if there is no md cell with colab => create empty md cell on top\n\ + else:\n\ + nb['cells'].insert(0, nbf.v4.new_markdown_cell())\n\ + nb['cells'][0]['source'] += '$download_badge_md'\n\ + nbf.write(nb, open('$notebook', 'w'))" + + python -c "exec (\"$python_command\")" + - name: Run notebook uses: fzimmermann89/run-notebook@v3 env: