From b5a3b5b4ecb0c904fdfb79ddce911bf1c1ba8b96 Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Fri, 6 Dec 2024 08:30:28 +0100 Subject: [PATCH 01/22] [pre-commit] pre-commit autoupdate (#573) Co-authored-by: pre-commit-ci[bot] <66853113+pre-commit-ci[bot]@users.noreply.github.com> --- .pre-commit-config.yaml | 4 ++-- src/mrpro/data/AcqInfo.py | 2 +- src/mrpro/data/MoveDataMixin.py | 2 +- src/mrpro/utils/slice_profiles.py | 2 +- src/mrpro/utils/typing.py | 2 +- src/mrpro/utils/unit_conversion.py | 12 ++++++------ 6 files changed, 12 insertions(+), 12 deletions(-) diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml index 2229c194f..3920fb5ae 100644 --- a/.pre-commit-config.yaml +++ b/.pre-commit-config.yaml @@ -14,14 +14,14 @@ repos: - id: mixed-line-ending - repo: https://github.com/astral-sh/ruff-pre-commit - rev: v0.7.2 + rev: v0.8.1 hooks: - id: ruff # linter args: [--fix] - id: ruff-format # formatter - repo: https://github.com/crate-ci/typos - rev: v1.27.0 + rev: typos-dict-v0.11.37 hooks: - id: typos diff --git a/src/mrpro/data/AcqInfo.py b/src/mrpro/data/AcqInfo.py index f5d677f97..b11e07d86 100644 --- a/src/mrpro/data/AcqInfo.py +++ b/src/mrpro/data/AcqInfo.py @@ -188,7 +188,7 @@ def tensor(data: np.ndarray) -> torch.Tensor: data = data.astype(np.int32) case np.uint32 | np.uint64: data = data.astype(np.int64) - # Remove any uncessary dimensions + # Remove any unnecessary dimensions return torch.tensor(np.squeeze(data)) def tensor_2d(data: np.ndarray) -> torch.Tensor: diff --git a/src/mrpro/data/MoveDataMixin.py b/src/mrpro/data/MoveDataMixin.py index 2ac3c1a58..99bcb3df5 100644 --- a/src/mrpro/data/MoveDataMixin.py +++ b/src/mrpro/data/MoveDataMixin.py @@ -121,7 +121,7 @@ def parse1( ) -> parsedType: return device, dtype, non_blocking, copy, memory_format - if args and isinstance(args[0], torch.Tensor) or 'tensor' in kwargs: + if (args and isinstance(args[0], torch.Tensor)) or 'tensor' in kwargs: # overload 3 ("tensor" specifies the dtype and device) device, dtype, non_blocking, copy, memory_format = parse3(*args, **kwargs) elif args and isinstance(args[0], torch.dtype): diff --git a/src/mrpro/utils/slice_profiles.py b/src/mrpro/utils/slice_profiles.py index 0e5fafc1f..8eaf95311 100644 --- a/src/mrpro/utils/slice_profiles.py +++ b/src/mrpro/utils/slice_profiles.py @@ -8,7 +8,7 @@ import torch from torch import Tensor -__all__ = ['SliceProfileBase', 'SliceGaussian', 'SliceSmoothedRectangular', 'SliceInterpolate'] +__all__ = ['SliceGaussian', 'SliceInterpolate', 'SliceProfileBase', 'SliceSmoothedRectangular'] class SliceProfileBase(abc.ABC, torch.nn.Module): diff --git a/src/mrpro/utils/typing.py b/src/mrpro/utils/typing.py index f90e18dab..96ad34d07 100644 --- a/src/mrpro/utils/typing.py +++ b/src/mrpro/utils/typing.py @@ -28,4 +28,4 @@ NestedSequence: TypeAlias = Any NumpyIndexerType: TypeAlias = Any -__all__ = ['TorchIndexerType', 'NumpyIndexerType', 'NestedSequence'] +__all__ = ['NestedSequence', 'NumpyIndexerType', 'TorchIndexerType'] diff --git a/src/mrpro/utils/unit_conversion.py b/src/mrpro/utils/unit_conversion.py index 0115bed47..5a1b5aaae 100644 --- a/src/mrpro/utils/unit_conversion.py +++ b/src/mrpro/utils/unit_conversion.py @@ -6,15 +6,15 @@ import torch __all__ = [ - 'ms_to_s', - 's_to_ms', - 'mm_to_m', - 'm_to_mm', + 'GYROMAGNETIC_RATIO_PROTON', 'deg_to_rad', - 'rad_to_deg', 'lamor_frequency_to_magnetic_field', + 'm_to_mm', 'magnetic_field_to_lamor_frequency', - 'GYROMAGNETIC_RATIO_PROTON', + 'mm_to_m', + 'ms_to_s', + 'rad_to_deg', + 's_to_ms', ] GYROMAGNETIC_RATIO_PROTON = 42.58 * 1e6 From b80141febc17f1d561c2ec0842c577bffdab0459 Mon Sep 17 00:00:00 2001 From: Patrick Schuenke <37338697+schuenke@users.noreply.github.com> Date: Mon, 9 Dec 2024 09:31:09 +0100 Subject: [PATCH 02/22] Fix trajectory scaling in KTrajectoryPulseq (#551) Co-authored-by: Felix F Zimmermann --- .../data/traj_calculators/KTrajectoryPulseq.py | 15 ++++++++++----- 1 file changed, 10 insertions(+), 5 deletions(-) diff --git a/src/mrpro/data/traj_calculators/KTrajectoryPulseq.py b/src/mrpro/data/traj_calculators/KTrajectoryPulseq.py index 7c843572d..598aa2184 100644 --- a/src/mrpro/data/traj_calculators/KTrajectoryPulseq.py +++ b/src/mrpro/data/traj_calculators/KTrajectoryPulseq.py @@ -54,13 +54,18 @@ def __call__(self, kheader: KHeader) -> KTrajectoryRawShape: raise ValueError('We currently only support constant number of samples') n_k0 = int(n_samples.item()) - def reshape_pulseq_traj(k_traj: torch.Tensor, encoding_size: int): - k_traj *= encoding_size / (2 * torch.max(torch.abs(k_traj))) + def rescale_and_reshape_traj(k_traj: torch.Tensor, encoding_size: int): + if encoding_size > 1 and torch.max(torch.abs(k_traj)) > 0: + k_traj = k_traj * encoding_size / (2 * torch.max(torch.abs(k_traj))) + else: + # We force k_traj to be 0 if encoding_size = 1. This is typically the case for kz in 2D sequences. + # However, it happens that seq.calculate_kspace() returns values != 0 (numerical noise) in such cases. + k_traj = torch.zeros_like(k_traj) return rearrange(k_traj, '(other k0) -> other k0', k0=n_k0) # rearrange k-space trajectory to match MRpro convention - kx = reshape_pulseq_traj(k_traj_adc[0], kheader.encoding_matrix.x) - ky = reshape_pulseq_traj(k_traj_adc[1], kheader.encoding_matrix.y) - kz = reshape_pulseq_traj(k_traj_adc[2], kheader.encoding_matrix.z) + kx = rescale_and_reshape_traj(k_traj_adc[0], kheader.encoding_matrix.x) + ky = rescale_and_reshape_traj(k_traj_adc[1], kheader.encoding_matrix.y) + kz = rescale_and_reshape_traj(k_traj_adc[2], kheader.encoding_matrix.z) return KTrajectoryRawShape(kz, ky, kx, self.repeat_detection_tolerance) From 31756e37e2dbd0875b1665bed8d6325246ae2627 Mon Sep 17 00:00:00 2001 From: Lunin Leonid Date: Tue, 10 Dec 2024 09:37:45 +0100 Subject: [PATCH 03/22] Add python 3.10 to badge and fix project keywords (#550) Co-authored-by: Patrick Schuenke <37338697+schuenke@users.noreply.github.com> --- README.md | 2 +- pyproject.toml | 15 ++++++++++++++- 2 files changed, 15 insertions(+), 2 deletions(-) diff --git a/README.md b/README.md index 422c4ef4d..c02b80deb 100644 --- a/README.md +++ b/README.md @@ -1,6 +1,6 @@ # MRpro -![Python](https://img.shields.io/badge/python-3.11%20%7C%203.12-blue) +![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) diff --git a/pyproject.toml b/pyproject.toml index 6f5d09048..70db1cafd 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -22,7 +22,20 @@ description = "MR image reconstruction and processing package specifically devel readme = "README.md" requires-python = ">=3.10,<3.14" dynamic = ["version"] -keywords = ["MRI, reconstruction, processing, PyTorch"] +keywords = ["MRI", + "qMRI", + "medical imaging", + "physics-informed learning", + "model-based reconstruction", + "quantitative", + "signal models", + "machine learning", + "deep learning", + "reconstruction", + "processing", + "Pulseq", + "PyTorch", +] authors = [ { name = "MRpro Team", email = "info@emerpro.de" }, { name = "Christoph Kolbitsch", email = "christoph.kolbitsch@ptb.de" }, From e15faa1673e3eb0789c2ac96e0a8892776145f5b Mon Sep 17 00:00:00 2001 From: Felix F Zimmermann Date: Tue, 10 Dec 2024 09:41:48 +0100 Subject: [PATCH 04/22] Add memory efficient reshape of broadcasted tensors (#557) --- src/mrpro/utils/__init__.py | 3 +- src/mrpro/utils/reshape.py | 101 ++++++++++++++++++++++++++++++++++++ tests/utils/test_reshape.py | 34 +++++++++++- 3 files changed, 136 insertions(+), 2 deletions(-) diff --git a/src/mrpro/utils/__init__.py b/src/mrpro/utils/__init__.py index 944100d54..b6e3c77ef 100644 --- a/src/mrpro/utils/__init__.py +++ b/src/mrpro/utils/__init__.py @@ -6,7 +6,7 @@ from mrpro.utils.remove_repeat import remove_repeat from mrpro.utils.zero_pad_or_crop import zero_pad_or_crop from mrpro.utils.split_idx import split_idx -from mrpro.utils.reshape import broadcast_right, unsqueeze_left, unsqueeze_right, reduce_view +from mrpro.utils.reshape import broadcast_right, unsqueeze_left, unsqueeze_right, reduce_view, reshape_broadcasted import mrpro.utils.unit_conversion __all__ = [ @@ -14,6 +14,7 @@ "fill_range_", "reduce_view", "remove_repeat", + "reshape_broadcasted", "slice_profiles", "smap", "split_idx", diff --git a/src/mrpro/utils/reshape.py b/src/mrpro/utils/reshape.py index 31d495afd..0b208381b 100644 --- a/src/mrpro/utils/reshape.py +++ b/src/mrpro/utils/reshape.py @@ -1,6 +1,8 @@ """Tensor reshaping utilities.""" from collections.abc import Sequence +from functools import lru_cache +from math import prod import torch @@ -99,3 +101,102 @@ def reduce_view(x: torch.Tensor, dim: int | Sequence[int] | None = None) -> torc for d, (oldsize, stride) in enumerate(zip(x.size(), stride, strict=True)) ] return torch.as_strided(x, newsize, stride) + + +@lru_cache +def _reshape_idx(old_shape: tuple[int, ...], new_shape: tuple[int, ...], old_stride: tuple[int, ...]) -> list[slice]: + """Get reshape reduce index (Cached helper function for reshape_broadcasted). + + This function tries to group axes from new_shape and old_shape into the smallest groups that have + the same number of elements, starting from the right. + If all axes of old shape of a group are stride=0 dimensions, we can reduce them. + + Example: + old_shape = (30, 2, 2, 3) + new_shape = (6, 5, 4, 3) + Will results in the groups (starting from the right): + - old: 3 new: 3 + - old: 2, 2 new: 4 + - old: 30 new: 6, 5 + Only the "old" groups are important. + If all axes that are grouped together in an "old" group are stride 0 (=broadcasted) + we can collapse them to singleton dimensions. + This function returns the indexer that either collapses dimensions to singleton or keeps all + elements, i.e. the slices in the returned list are all either slice(1) or slice(None). + """ + idx = [] + pointer_old, pointer_new = len(old_shape) - 1, len(new_shape) - 1 # start from the right + while pointer_old >= 0: + product_new, product_old = 1, 1 # the number of elements in the current "new" and "old" group + group: list[int] = [] + while product_old != product_new or not group: + if product_old <= product_new: + # increase "old" group + product_old *= old_shape[pointer_old] + group.append(pointer_old) + pointer_old -= 1 + else: + # increase "new" group + # we don't need to track the new group, the number of elemeents covered. + product_new *= new_shape[pointer_new] + pointer_new -= 1 + # we found a group. now we need to decide what to do. + if all(old_stride[d] == 0 for d in group): + # all dimensions are broadcasted + # -> reduce to singleton + idx.extend([slice(1)] * len(group)) + else: + # preserve dimension + idx.extend([slice(None)] * len(group)) + idx = idx[::-1] # we worked right to left, but our index should be left to right + return idx + + +def reshape_broadcasted(tensor: torch.Tensor, *shape: int) -> torch.Tensor: + """Reshape a tensor while preserving broadcasted (stride 0) dimensions where possible. + + Parameters + ---------- + tensor + The input tensor to reshape. + shape + The target shape for the tensor. One of the values can be `-1` and its size will be inferred. + + Returns + ------- + A tensor reshaped to the target shape, preserving broadcasted dimensions where feasible. + + """ + try: + # if we can view the tensor directly, it will preserve broadcasting + return tensor.view(shape) + except RuntimeError: + # we cannot do a view, we need to do more work: + + # -1 means infer size, i.e. the remaining elements of the input not already covered by the other axes. + negative_ones = shape.count(-1) + size = tensor.shape.numel() + if not negative_ones: + if prod(shape) != size: + # use same exception as pytorch + raise RuntimeError(f"shape '{list(shape)}' is invalid for input of size {size}") from None + elif negative_ones > 1: + raise RuntimeError('only one dimension can be inferred') from None + elif negative_ones == 1: + # we need to figure out the size of the "-1" dimension + known_size = -prod(shape) # negative, is it includes the -1 + if size % known_size: + # non integer result. no possible size of the -1 axis exists. + raise RuntimeError(f"shape '{list(shape)}' is invalid for input of size {size}") from None + shape = tuple(size // known_size if s == -1 else s for s in shape) + + # most of the broadcasted dimensions can be preserved: only dimensions that are joined with non + # broadcasted dimensions can not be preserved and must be made contiguous. + # all dimensions that can be preserved as broadcasted are first collapsed to singleton, + # such that contiguous does not create copies along these axes. + idx = _reshape_idx(tensor.shape, shape, tensor.stride()) + # make contiguous only in dimensions in which broadcasting cannot be preserved + semicontiguous = tensor[idx].contiguous() + # finally, we can expand the broadcasted dimensions to the requested shape + semicontiguous = semicontiguous.expand(tensor.shape) + return semicontiguous.view(shape) diff --git a/tests/utils/test_reshape.py b/tests/utils/test_reshape.py index dd57b8feb..3a2cc0cf5 100644 --- a/tests/utils/test_reshape.py +++ b/tests/utils/test_reshape.py @@ -1,7 +1,8 @@ """Tests for reshaping utilities.""" +import pytest import torch -from mrpro.utils import broadcast_right, reduce_view, unsqueeze_left, unsqueeze_right +from mrpro.utils import broadcast_right, reduce_view, reshape_broadcasted, unsqueeze_left, unsqueeze_right from tests import RandomGenerator @@ -51,3 +52,34 @@ def test_reduce_view(): reduced_one_pos = reduce_view(tensor, 0) assert reduced_one_pos.shape == (1, 2, 3, 4, 5, 6) assert torch.equal(reduced_one_pos.expand_as(tensor), tensor) + + +@pytest.mark.parametrize( + ('shape', 'expand_shape', 'permute', 'final_shape', 'expected_stride'), + [ + ((1, 2, 3, 1, 1), (1, 2, 3, 4, 5), (0, 2, 1, 3, 4), (1, 6, 2, 2, 5), (6, 1, 0, 0, 0)), + ((1, 2, 1), (100, 2, 2), (0, 1, 2), (100, 4), (0, 1)), + ((1, 1, 1, 1), (2, 3, 4, 5), (2, 3, 0, 1), (1, 2, 6, 10, 1), (0, 0, 0, 0, 0)), + ((1, 2, 3), (1, -1, 3), (0, 1, 2), (6,), (1,)), + ], +) +def test_reshape_broadcasted(shape, expand_shape, permute, final_shape, expected_stride): + """Test reshape_broadcasted""" + rng = RandomGenerator(0) + tensor = rng.float32_tensor(shape).expand(*expand_shape).permute(*permute) + reshaped = reshape_broadcasted(tensor, *final_shape) + expected_values = tensor.reshape(*final_shape) + assert reshaped.shape == expected_values.shape + assert reshaped.stride() == expected_stride + assert torch.equal(reshaped, expected_values) + + +def test_reshape_broadcasted_fail(): + """Test reshape_broadcasted with invalid input""" + a = torch.ones(2) + with pytest.raises(RuntimeError, match='invalid'): + reshape_broadcasted(a, 3) + with pytest.raises(RuntimeError, match='invalid'): + reshape_broadcasted(a, -1, -3) + with pytest.raises(RuntimeError, match='only one dimension'): + reshape_broadcasted(a, -1, -1) From 28038b11d852fb01df76e33f220bac5cb2de7d4c Mon Sep 17 00:00:00 2001 From: Patrick Schuenke <37338697+schuenke@users.noreply.github.com> Date: Tue, 10 Dec 2024 17:31:22 +0100 Subject: [PATCH 05/22] Release v0.241210 (#576) --- src/mrpro/VERSION | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/mrpro/VERSION b/src/mrpro/VERSION index 1e8c1c10e..bcc11aaab 100644 --- a/src/mrpro/VERSION +++ b/src/mrpro/VERSION @@ -1 +1 @@ -0.241126 +0.241210 From e6568dcc56b204cf90fc6cb28c6cb318911e40fe Mon Sep 17 00:00:00 2001 From: Lunin Leonid Date: Wed, 11 Dec 2024 14:08:06 +0100 Subject: [PATCH 06/22] Omit notebook subfolders in docs (#549) --- .github/workflows/docs.yml | 1 + docs/source/examples.rst | 2 +- 2 files changed, 2 insertions(+), 1 deletion(-) diff --git a/.github/workflows/docs.yml b/.github/workflows/docs.yml index cc5ff5458..d88c2ea4b 100644 --- a/.github/workflows/docs.yml +++ b/.github/workflows/docs.yml @@ -155,6 +155,7 @@ jobs: uses: actions/download-artifact@v4 with: path: ./docs/source/_notebooks/ + merge-multiple: true - name: Build docs run: | diff --git a/docs/source/examples.rst b/docs/source/examples.rst index fe24a60a8..95c654781 100644 --- a/docs/source/examples.rst +++ b/docs/source/examples.rst @@ -9,4 +9,4 @@ All of the notebooks can directly be run via binder or colab from the repo. :caption: Contents: :glob: - _notebooks/*/* + _notebooks/* From f6f921d116cef0bcedace7c476bd26c3ab16bb5c Mon Sep 17 00:00:00 2001 From: Lunin Leonid Date: Mon, 16 Dec 2024 15:03:51 +0100 Subject: [PATCH 07/22] Fetch docker names without github token (#555) This will allow the tests to run on forks. Co-authored-by: Christoph Kolbitsch Co-authored-by: Felix F Zimmermann --- .github/workflows/docker.yml | 53 ++++++++++++++---------------------- .github/workflows/docs.yml | 10 +++---- .github/workflows/pytest.yml | 18 ++++++++---- 3 files changed, 37 insertions(+), 44 deletions(-) diff --git a/.github/workflows/docker.yml b/.github/workflows/docker.yml index fe4f2ae03..fffcfe30a 100644 --- a/.github/workflows/docker.yml +++ b/.github/workflows/docker.yml @@ -15,8 +15,7 @@ jobs: runs-on: ubuntu-latest outputs: docker_toml: ${{ steps.filter.outputs.docker_toml }} - dockerfiles: ${{ steps.set-matrix.outputs.dockerfiles }} - imagenames: ${{ steps.set-matrix.outputs.imagenames }} + docker_tasks: ${{ steps.set-matrix.outputs.docker_tasks }} steps: - name: Checkout mrpro repo uses: actions/checkout@v4 @@ -39,20 +38,17 @@ jobs: id: set-matrix if: steps.filter.outputs.docker_toml == 'true' || github.event_name == 'push' run: | - cd ./docker/ - ls - dockerfiles=$(ls Dockerfile_* | jq -R -s -c 'split("\n")[:-1]') - echo "dockerfiles: $dockerfiles" - echo "dockerfiles=$dockerfiles" >> $GITHUB_OUTPUT - imagenames=$(ls Dockerfile_* | sed -e 's/Dockerfile_/ghcr.io\/ptb-mr\/mrpro_/' | jq -R -s -c 'split("\n")[:-1]') - echo "image names: $imagenames" - echo "imagenames=$imagenames" >> $GITHUB_OUTPUT + # docker_tasks is a list of pairs (dictionaries) with keys 'filepath' and 'image_name' like: + # [{"filepath": "docker/Dockerfile_x1", "image_name": "ghcr.io/ptb-mr/mrpro_x1"}, ...] + docker_tasks=$(find docker -type f -name 'Dockerfile*' | jq -R -s -c 'split("\n")[:-1]' | \ + jq -r -c 'map({filepath: ., image_name: . | sub("docker/Dockerfile"; "ghcr.io\/ptb-mr\/mrpro")})') + echo "docker_tasks: $docker_tasks" + echo "docker_tasks=$docker_tasks" >> $GITHUB_OUTPUT - name: Dockerfile overview if: steps.filter.outputs.docker_toml == 'true' || github.event_name == 'push' run: | - echo "final list of dockerfiles: ${{ steps.set-matrix.outputs.dockerfiles }}" - echo "final list of images: ${{ steps.set-matrix.outputs.imagenames }}" + echo "final list of docker_tasks: ${{ steps.set-matrix.outputs.docker_tasks }}" push_test: name: Create test images and push to GCR @@ -64,16 +60,15 @@ jobs: strategy: fail-fast: false matrix: - dockerfile: ${{ fromJson(needs.get_dockerfiles.outputs.dockerfiles) }} + docker_task: ${{ fromJson(needs.get_dockerfiles.outputs.docker_tasks) }} steps: - name: Checkout uses: actions/checkout@v4 - - name: Create image name + - name: Get image basename id: image_name run: | - dockerfile=${{ matrix.dockerfile }} - echo "image_name=${dockerfile/Dockerfile_/ghcr.io/ptb-mr/mrpro_}" >> $GITHUB_OUTPUT + echo "dockerfile_basename=$(basename ${{ matrix.docker_task.filepath }})" >> $GITHUB_OUTPUT - name: Set up Docker Buildx uses: docker/setup-buildx-action@v3 @@ -89,11 +84,11 @@ jobs: uses: docker/build-push-action@v6 with: context: . - cache-from: type=gha,scope=${{ matrix.dockerfile }} - cache-to: type=gha,mode=max,scope=${{ matrix.dockerfile }} - file: ./docker/${{ matrix.dockerfile }} + cache-from: type=gha,scope=$${{ steps.image_name.outputs.dockerfile_basename }}) + cache-to: type=gha,mode=max,scope=${{ steps.image_name.outputs.dockerfile_basename }} + file: ${{ matrix.docker_task.filepath }} push: true - tags: ${{ steps.image_name.outputs.image_name }}:test + tags: ${{ matrix.docker_task.image_name }}:test test: name: Test docker containers @@ -105,10 +100,10 @@ jobs: contents: write strategy: matrix: - imagename: ${{ fromJson(needs.get_dockerfiles.outputs.imagenames) }} + docker_task: ${{ fromJson(needs.get_dockerfiles.outputs.docker_tasks) }} # runs within Docker container container: - image: ${{ matrix.imagename }}:test + image: ${{ matrix.docker_task.image_name }}:test options: --user runner steps: @@ -134,17 +129,11 @@ jobs: packages: write strategy: matrix: - dockerfile: ${{ fromJson(needs.get_dockerfiles.outputs.dockerfiles) }} + docker_task: ${{ fromJson(needs.get_dockerfiles.outputs.docker_tasks) }} steps: - name: Checkout uses: actions/checkout@v4 - - name: Create image name - id: image_name - run: | - dockerfile=${{ matrix.dockerfile }} - echo "image_name=${dockerfile/Dockerfile_/ghcr.io/ptb-mr/mrpro_}" >> $GITHUB_OUTPUT - - name: Login to GitHub Packages uses: docker/login-action@v3 with: @@ -154,6 +143,6 @@ jobs: - name: Pull and push Docker image with new tag run: | - docker pull ${{ steps.image_name.outputs.image_name }}:test - docker tag ${{ steps.image_name.outputs.image_name }}:test ${{ steps.image_name.outputs.image_name }}:latest - docker push ${{ steps.image_name.outputs.image_name }}:latest + docker pull ${{ matrix.docker_task.image_name }}:test + docker tag ${{ matrix.docker_task.image_name }}:test ${{ matrix.docker_task.image_name }}:latest + docker push ${{ matrix.docker_task.image_name}}:latest diff --git a/.github/workflows/docs.yml b/.github/workflows/docs.yml index d88c2ea4b..fae4bf3ba 100644 --- a/.github/workflows/docs.yml +++ b/.github/workflows/docs.yml @@ -71,9 +71,7 @@ jobs: - id: set-matrix run: | - cd ./examples/ - ls - echo "notebooks=$(ls *.ipynb | jq -R -s -c 'split("\n")[:-1]')" >> $GITHUB_OUTPUT + echo "notebooks=$(find examples -type f -name '*.ipynb' | jq -R -s -c 'split("\n")[:-1]')" >> $GITHUB_OUTPUT - name: Notebook overview run: | @@ -113,14 +111,14 @@ jobs: env: RUNNER: ${{ toJson(runner) }} with: - notebook: ./examples/${{ matrix.notebook }} + notebook: ./${{ matrix.notebook }} - name: Get artifact names id: artifact_names run: | notebook=${{ matrix.notebook }} - echo "ARTIFACT_NAME=${notebook/.ipynb/}" >> $GITHUB_OUTPUT - echo "IPYNB_EXECUTED=${notebook}" >> $GITHUB_OUTPUT + echo "ARTIFACT_NAME=$(basename ${notebook/.ipynb})" >> $GITHUB_OUTPUT + echo "IPYNB_EXECUTED=$(basename $notebook)" >> $GITHUB_OUTPUT - name: Upload notebook uses: actions/upload-artifact@v4 diff --git a/.github/workflows/pytest.yml b/.github/workflows/pytest.yml index 53d1343c7..21079c85a 100644 --- a/.github/workflows/pytest.yml +++ b/.github/workflows/pytest.yml @@ -15,18 +15,24 @@ jobs: outputs: imagenames: ${{ steps.set-matrix.outputs.imagenames }} steps: + - name: Checkout + uses: actions/checkout@v4 + - name: Retrieve Docker Image Names id: set-matrix - env: - GH_TOKEN: ${{ secrets.GHCR_TOKEN }} run: | - imagenames=$(curl -s --request GET \ - --url "https://api.github.com/orgs/PTB-MR/packages?package_type=container" \ - --header "Authorization: Bearer $GH_TOKEN" | jq -r '.[].name') + # search for Dockerfile* in the docker directory, replace "Dockerfile" prefix with "mrpro" and to imagenames + imagenames=$(find docker -type f -name 'Dockerfile*' -exec basename {} \; | sed 's/^Dockerfile/mrpro/') echo "image names: $imagenames" + # if imagenames is empty - fail the workflow + if [ -z "$imagenames" ]; then + echo "No Dockerfiles found in the docker directory. Exiting..." + exit 1 + fi + imagenames_latest=() - for image in $(echo $imagenames) + for image in $imagenames do echo "checking $image ..." if docker manifest inspect "ghcr.io/ptb-mr/"$image":latest" >/dev/null; then From da02db8c5bc4b951242c69c0d6e92d0f172e9aeb Mon Sep 17 00:00:00 2001 From: Patrick Schuenke <37338697+schuenke@users.noreply.github.com> Date: Mon, 16 Dec 2024 16:08:15 +0100 Subject: [PATCH 08/22] Add logo to readme and docs (#583) --- README.md | 6 +- docs/source/_static/logo.svg | 218 +++++++++++++++++++++++++++++ docs/source/_static/logo_white.svg | 218 +++++++++++++++++++++++++++++ docs/source/conf.py | 3 +- docs/source/index.rst | 5 + 5 files changed, 446 insertions(+), 4 deletions(-) create mode 100644 docs/source/_static/logo.svg create mode 100644 docs/source/_static/logo_white.svg diff --git a/README.md b/README.md index c02b80deb..49d7e2c3d 100644 --- a/README.md +++ b/README.md @@ -1,4 +1,6 @@ -# MRpro +

+MRpro logo +


![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) @@ -84,5 +86,3 @@ We are looking forward to your contributions via Pull-Requests. 4. Setup pre-commit hook: ``` pre-commit install ``` Please look at our [contributor guide](https://ptb-mr.github.io/mrpro/contributor_guide.html) for more information on the repository structure, naming conventions, and other useful information. - - diff --git a/docs/source/_static/logo.svg b/docs/source/_static/logo.svg new file mode 100644 index 000000000..29c7eeabf --- /dev/null +++ b/docs/source/_static/logo.svg @@ -0,0 +1,218 @@ + + + +MRpr diff --git a/docs/source/_static/logo_white.svg b/docs/source/_static/logo_white.svg new file mode 100644 index 000000000..d0d5b3527 --- /dev/null +++ b/docs/source/_static/logo_white.svg @@ -0,0 +1,218 @@ + + + +MRpr diff --git a/docs/source/conf.py b/docs/source/conf.py index 02f75262b..764b75f13 100644 --- a/docs/source/conf.py +++ b/docs/source/conf.py @@ -60,9 +60,10 @@ html_show_sphinx = False html_static_path = ['_static'] html_css_files = ['custom.css'] +html_logo = '_static/logo_white.svg' html_sidebars = {'**': ['search-field', 'sidebar-nav-bs']} html_theme_options = { - 'logo': {'text': name}, + 'logo_only': True, 'pygment_light_style': 'default', 'pygment_dark_style': 'github-dark', 'show_toc_level': 3, diff --git a/docs/source/index.rst b/docs/source/index.rst index a08619105..294847210 100644 --- a/docs/source/index.rst +++ b/docs/source/index.rst @@ -1,5 +1,10 @@ .. MRpro documentation +.. image:: _static/logo.svg + :align: center + :width: 400 +| + Welcome to MRpro's documentation! ================================= From 4096ead804cd155d19394e6cdbbb40f124619a07 Mon Sep 17 00:00:00 2001 From: Christoph Kolbitsch Date: Tue, 17 Dec 2024 08:20:04 +0100 Subject: [PATCH 09/22] Modification and test bug fix for TrajType (#584) Co-authored-by: Felix F Zimmermann --- src/mrpro/data/enums.py | 8 ++++---- tests/data/test_trajectory.py | 26 ++++++++++++++++++++------ 2 files changed, 24 insertions(+), 10 deletions(-) diff --git a/src/mrpro/data/enums.py b/src/mrpro/data/enums.py index b14de615d..73fc76a90 100644 --- a/src/mrpro/data/enums.py +++ b/src/mrpro/data/enums.py @@ -1,6 +1,6 @@ """All acquisition enums.""" -from enum import Enum, Flag, auto +from enum import Enum, Flag, IntFlag, auto class AcqFlags(Flag): @@ -124,8 +124,8 @@ class CalibrationMode(Enum): OTHER = 'other' -class TrajType(Flag): +class TrajType(IntFlag): """Special Properties of the Trajectory.""" - SINGLEVALUE = 1 - ONGRID = 2 + SINGLEVALUE = auto() + ONGRID = auto() diff --git a/tests/data/test_trajectory.py b/tests/data/test_trajectory.py index 1061a93be..33ae0266b 100644 --- a/tests/data/test_trajectory.py +++ b/tests/data/test_trajectory.py @@ -154,9 +154,16 @@ def test_ktype_along_kzyx(im_shape, k_shape, nkx, nky, nkz, type_kx, type_ky, ty trajectory = create_traj(k_shape, nkx, nky, nkz, type_kx, type_ky, type_kz) # Find out the type of the kz, ky and kz dimensions - single_value_dims = [d for d, s in zip((-3, -2, -1), (type_kz, type_ky, type_kx), strict=True) if s == 'z'] - on_grid_dims = [d for d, s in zip((-3, -2, -1), (type_kz, type_ky, type_kx), strict=True) if s == 'uf'] - not_on_grid_dims = [d for d, s in zip((-3, -2, -1), (type_kz, type_ky, type_kx), strict=True) if s == 'nuf'] + single_value_dims = [d for d, s in zip((-3, -2, -1), (type_kz, type_ky, type_kx), strict=True) if s == 'zero'] + on_grid_dims = [d for d, s in zip((-3, -2, -1), (type_kz, type_ky, type_kx), strict=True) if s == 'uniform'] + not_on_grid_dims = [d for d, s in zip((-3, -2, -1), (type_kz, type_ky, type_kx), strict=True) if s == 'non-uniform'] + + # Make sure all dimensions are covered + if len(single_value_dims) + len(on_grid_dims) + len(not_on_grid_dims) != 3: + raise ValueError( + f'{single_value_dims=}, {on_grid_dims=} and {not_on_grid_dims=} do not cover all dimensions. ', + 'There must be an error in the test itself.', + ) # check dimensions which are of shape 1 and do not need any transform assert all(trajectory.type_along_kzyx[dim] & TrajType.SINGLEVALUE for dim in single_value_dims) @@ -178,9 +185,16 @@ def test_ktype_along_k210(im_shape, k_shape, nkx, nky, nkz, type_kx, type_ky, ty trajectory = create_traj(k_shape, nkx, nky, nkz, type_kx, type_ky, type_kz) # Find out the type of the k2, k1 and k0 dimensions - single_value_dims = [d for d, s in zip((-3, -2, -1), (type_k2, type_k1, type_k0), strict=True) if s == 'z'] - on_grid_dims = [d for d, s in zip((-3, -2, -1), (type_k2, type_k1, type_k0), strict=True) if s == 'uf'] - not_on_grid_dims = [d for d, s in zip((-3, -2, -1), (type_k2, type_k1, type_k0), strict=True) if s == 'nuf'] + single_value_dims = [d for d, s in zip((-3, -2, -1), (type_k2, type_k1, type_k0), strict=True) if s == 'zero'] + on_grid_dims = [d for d, s in zip((-3, -2, -1), (type_k2, type_k1, type_k0), strict=True) if s == 'uniform'] + not_on_grid_dims = [d for d, s in zip((-3, -2, -1), (type_k2, type_k1, type_k0), strict=True) if s == 'non-uniform'] + + # Make sure all dimensions are covered + if len(single_value_dims) + len(on_grid_dims) + len(not_on_grid_dims) != 3: + raise ValueError( + f'{single_value_dims=}, {on_grid_dims=} and {not_on_grid_dims=} do not cover all dimensions. ', + 'There must be an error in the test itself.', + ) # check dimensions which are of shape 1 and do not need any transform assert all(trajectory.type_along_k210[dim] & TrajType.SINGLEVALUE for dim in single_value_dims) From 92e3193624e6ee93fd7baaee08de4dfe42279179 Mon Sep 17 00:00:00 2001 From: Felix F Zimmermann Date: Tue, 17 Dec 2024 17:07:25 +0100 Subject: [PATCH 10/22] Move trajectory scaling into KTrajectory (#582) Co-authored-by: pre-commit-ci[bot] <66853113+pre-commit-ci[bot]@users.noreply.github.com> Co-authored-by: Patrick Schuenke --- src/mrpro/data/KTrajectory.py | 45 +++++++++++++---- src/mrpro/data/KTrajectoryRawShape.py | 49 +++++++++++++++++++ .../traj_calculators/KTrajectoryPulseq.py | 28 ++++------- tests/data/_PulseqRadialTestSeq.py | 4 +- tests/data/test_traj_calculators.py | 8 +-- 5 files changed, 101 insertions(+), 33 deletions(-) diff --git a/src/mrpro/data/KTrajectory.py b/src/mrpro/data/KTrajectory.py index 8cbe207cc..f8bab4bae 100644 --- a/src/mrpro/data/KTrajectory.py +++ b/src/mrpro/data/KTrajectory.py @@ -1,6 +1,7 @@ """KTrajectory dataclass.""" from dataclasses import dataclass +from typing import Literal import numpy as np import torch @@ -8,6 +9,7 @@ from mrpro.data.enums import TrajType from mrpro.data.MoveDataMixin import MoveDataMixin +from mrpro.data.SpatialDimension import SpatialDimension from mrpro.utils import remove_repeat from mrpro.utils.summarize_tensorvalues import summarize_tensorvalues @@ -69,29 +71,52 @@ def from_tensor( cls, tensor: torch.Tensor, stack_dim: int = 0, - repeat_detection_tolerance: float | None = 1e-8, + axes_order: Literal['zxy', 'zyx', 'yxz', 'yzx', 'xyz', 'xzy'] = 'zyx', + repeat_detection_tolerance: float | None = 1e-6, grid_detection_tolerance: float = 1e-3, + scaling_matrix: SpatialDimension | None = None, ) -> Self: """Create a KTrajectory from a tensor representation of the trajectory. - Reduces repeated dimensions to singletons if repeat_detection_tolerance - is not set to None. - + Reduces repeated dimensions to singletons if repeat_detection_tolerance is not set to None. Parameters ---------- tensor The tensor representation of the trajectory. - This should be a 5-dim tensor, with (kz,ky,kx) stacked in this order along stack_dim + This should be a 5-dim tensor, with (kz, ky, kx) stacked in this order along `stack_dim`. stack_dim - The dimension in the tensor the directions have been stacked along. + The dimension in the tensor along which the directions are stacked. + axes_order + The order of the axes in the tensor. The MRpro convention is 'zyx'. repeat_detection_tolerance - detects if broadcasting can be used, i.e. if dimensions are repeated. - Set to None to disable. + Tolerance for detecting repeated dimensions (broadcasting). + If trajectory points differ by less than this value, they are considered identical. + Set to None to disable this feature. grid_detection_tolerance - tolerance to detect if trajectory points are on integer grid positions + Tolerance for detecting whether trajectory points align with integer grid positions. + This tolerance is applied after rescaling if `scaling_matrix` is provided. + scaling_matrix + If a scaling matrix is provided, the trajectory is rescaled to fit within + the dimensions of the matrix. If not provided, the trajectory remains unchanged. + """ - kz, ky, kx = torch.unbind(tensor, dim=stack_dim) + ks = tensor.unbind(dim=stack_dim) + kz, ky, kx = (ks[axes_order.index(axis)] for axis in 'zyx') + + def rescale(k: torch.Tensor, size: float) -> torch.Tensor: + max_abs_range = 2 * k.abs().max() + if size < 2 or max_abs_range < 1e-6: + # a single encoding point should be at zero + # avoid division by zero + return torch.zeros_like(k) + return k * (size / max_abs_range) + + if scaling_matrix is not None: + kz = rescale(kz, scaling_matrix.z) + ky = rescale(ky, scaling_matrix.y) + kx = rescale(kx, scaling_matrix.x) + return cls( kz, ky, diff --git a/src/mrpro/data/KTrajectoryRawShape.py b/src/mrpro/data/KTrajectoryRawShape.py index 3730e6669..e75ad8c16 100644 --- a/src/mrpro/data/KTrajectoryRawShape.py +++ b/src/mrpro/data/KTrajectoryRawShape.py @@ -1,13 +1,16 @@ """KTrajectoryRawShape dataclass.""" from dataclasses import dataclass +from typing import Literal import numpy as np import torch from einops import rearrange +from typing_extensions import Self from mrpro.data.KTrajectory import KTrajectory from mrpro.data.MoveDataMixin import MoveDataMixin +from mrpro.data.SpatialDimension import SpatialDimension @dataclass(slots=True, frozen=True) @@ -32,6 +35,52 @@ class KTrajectoryRawShape(MoveDataMixin): repeat_detection_tolerance: None | float = 1e-3 """tolerance for repeat detection. Set to None to disable.""" + @classmethod + def from_tensor( + cls, + tensor: torch.Tensor, + stack_dim: int = 0, + axes_order: Literal['zxy', 'zyx', 'yxz', 'yzx', 'xyz', 'xzy'] = 'zyx', + repeat_detection_tolerance: float | None = 1e-6, + scaling_matrix: SpatialDimension | None = None, + ) -> Self: + """Create a KTrajectoryRawShape from a tensor representation of the trajectory. + + Parameters + ---------- + tensor + The tensor representation of the trajectory. + This should be a 5-dim tensor, with (kz, ky, kx) stacked in this order along `stack_dim`. + stack_dim + The dimension in the tensor along which the directions are stacked. + axes_order + The order of the axes in the tensor. The MRpro convention is 'zyx'. + repeat_detection_tolerance + Tolerance for detecting repeated dimensions (broadcasting). + If trajectory points differ by less than this value, they are considered identical. + Set to None to disable this feature. + scaling_matrix + If a scaling matrix is provided, the trajectory is rescaled to fit within + the dimensions of the matrix. If not provided, the trajectory remains unchanged. + """ + ks = tensor.unbind(dim=stack_dim) + kz, ky, kx = (ks[axes_order.index(axis)] for axis in 'zyx') + + def rescale(k: torch.Tensor, size: float) -> torch.Tensor: + max_abs_range = 2 * k.abs().max() + if size < 2 or max_abs_range < 1e-6: + # a single encoding point should be at zero + # avoid division by zero + return torch.zeros_like(k) + return k * (size / max_abs_range) + + if scaling_matrix is not None: + kz = rescale(kz, scaling_matrix.z) + ky = rescale(ky, scaling_matrix.y) + kx = rescale(kx, scaling_matrix.x) + + return cls(kz, ky, kx, repeat_detection_tolerance=repeat_detection_tolerance) + def sort_and_reshape( self, sort_idx: np.ndarray, diff --git a/src/mrpro/data/traj_calculators/KTrajectoryPulseq.py b/src/mrpro/data/traj_calculators/KTrajectoryPulseq.py index 598aa2184..a446464a8 100644 --- a/src/mrpro/data/traj_calculators/KTrajectoryPulseq.py +++ b/src/mrpro/data/traj_calculators/KTrajectoryPulseq.py @@ -2,7 +2,6 @@ from pathlib import Path -import pypulseq as pp import torch from einops import rearrange @@ -40,8 +39,10 @@ def __call__(self, kheader: KHeader) -> KTrajectoryRawShape: ------- trajectory of type KTrajectoryRawShape """ + from pypulseq import Sequence + # create PyPulseq Sequence object and read .seq file - seq = pp.Sequence() + seq = Sequence() seq.read(file_path=str(self.seq_path)) # calculate k-space trajectory using PyPulseq @@ -52,20 +53,11 @@ def __call__(self, kheader: KHeader) -> KTrajectoryRawShape: n_samples = torch.unique(n_samples) if len(n_samples) > 1: raise ValueError('We currently only support constant number of samples') - n_k0 = int(n_samples.item()) - - def rescale_and_reshape_traj(k_traj: torch.Tensor, encoding_size: int): - if encoding_size > 1 and torch.max(torch.abs(k_traj)) > 0: - k_traj = k_traj * encoding_size / (2 * torch.max(torch.abs(k_traj))) - else: - # We force k_traj to be 0 if encoding_size = 1. This is typically the case for kz in 2D sequences. - # However, it happens that seq.calculate_kspace() returns values != 0 (numerical noise) in such cases. - k_traj = torch.zeros_like(k_traj) - return rearrange(k_traj, '(other k0) -> other k0', k0=n_k0) - - # rearrange k-space trajectory to match MRpro convention - kx = rescale_and_reshape_traj(k_traj_adc[0], kheader.encoding_matrix.x) - ky = rescale_and_reshape_traj(k_traj_adc[1], kheader.encoding_matrix.y) - kz = rescale_and_reshape_traj(k_traj_adc[2], kheader.encoding_matrix.z) - return KTrajectoryRawShape(kz, ky, kx, self.repeat_detection_tolerance) + k_traj_reshaped = rearrange(k_traj_adc, 'xyz (other k0) -> xyz other k0', k0=int(n_samples.item())) + return KTrajectoryRawShape.from_tensor( + k_traj_reshaped, + axes_order='xyz', + scaling_matrix=kheader.encoding_matrix, + repeat_detection_tolerance=self.repeat_detection_tolerance, + ) diff --git a/tests/data/_PulseqRadialTestSeq.py b/tests/data/_PulseqRadialTestSeq.py index 82cab5577..0cd9b6032 100644 --- a/tests/data/_PulseqRadialTestSeq.py +++ b/tests/data/_PulseqRadialTestSeq.py @@ -29,7 +29,9 @@ def __init__(self, seq_filename: str, n_x=256, n_spokes=10): system = pypulseq.Opts() rf, gz, _ = pypulseq.make_sinc_pulse(flip_angle=0.1, slice_thickness=1e-3, system=system, return_gz=True) - gx = pypulseq.make_trapezoid(channel='x', flat_area=n_x * delta_k, flat_time=2e-3, system=system) + gx = pypulseq.make_trapezoid( + channel='x', flat_area=n_x * delta_k, flat_time=n_x * system.grad_raster_time, system=system + ) adc = pypulseq.make_adc(num_samples=n_x, duration=gx.flat_time, delay=gx.rise_time, system=system) gx_pre = pypulseq.make_trapezoid(channel='x', area=-gx.area / 2 - delta_k / 2, duration=2e-3, system=system) gz_reph = pypulseq.make_trapezoid(channel='z', area=-gz.area / 2, duration=2e-3, system=system) diff --git a/tests/data/test_traj_calculators.py b/tests/data/test_traj_calculators.py index 7ddcb30aa..8ac3d0b71 100644 --- a/tests/data/test_traj_calculators.py +++ b/tests/data/test_traj_calculators.py @@ -260,11 +260,11 @@ def test_KTrajectoryPulseq_validseq_random_header(pulseq_example_rad_seq, valid_ trajectory_calculator = KTrajectoryPulseq(seq_path=pulseq_example_rad_seq.seq_filename) trajectory = trajectory_calculator(kheader=valid_rad2d_kheader) - kx_test = pulseq_example_rad_seq.traj_analytical.kx.squeeze(0).squeeze(0) - kx_test *= valid_rad2d_kheader.encoding_matrix.x / (2 * torch.max(torch.abs(kx_test))) + kx_test = pulseq_example_rad_seq.traj_analytical.kx.squeeze() + kx_test = kx_test * valid_rad2d_kheader.encoding_matrix.x / (2 * kx_test.abs().max()) - ky_test = pulseq_example_rad_seq.traj_analytical.ky.squeeze(0).squeeze(0) - ky_test *= valid_rad2d_kheader.encoding_matrix.y / (2 * torch.max(torch.abs(ky_test))) + ky_test = pulseq_example_rad_seq.traj_analytical.ky.squeeze() + ky_test = ky_test * valid_rad2d_kheader.encoding_matrix.y / (2 * ky_test.abs().max()) torch.testing.assert_close(trajectory.kx.to(torch.float32), kx_test.to(torch.float32), atol=1e-2, rtol=1e-3) torch.testing.assert_close(trajectory.ky.to(torch.float32), ky_test.to(torch.float32), atol=1e-2, rtol=1e-3) From 9dc116731154cade0029c8c58ee97013918d1a98 Mon Sep 17 00:00:00 2001 From: Patrick Schuenke <37338697+schuenke@users.noreply.github.com> Date: Tue, 17 Dec 2024 17:18:01 +0100 Subject: [PATCH 11/22] Release v0.241217 (#586) --- src/mrpro/VERSION | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/mrpro/VERSION b/src/mrpro/VERSION index bcc11aaab..8c5b741c6 100644 --- a/src/mrpro/VERSION +++ b/src/mrpro/VERSION @@ -1 +1 @@ -0.241210 +0.241217 From bc2e05a3048ee9b880d00cdf6a73973d3f344fa8 Mon Sep 17 00:00:00 2001 From: Christoph Kolbitsch Date: Sat, 4 Jan 2025 15:50:59 +0100 Subject: [PATCH 12/22] 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 13/22] 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 14/22] 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 15/22] 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 16/22] 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: From d50b5aa6bad16de70a34c0c71c75c3079357eda0 Mon Sep 17 00:00:00 2001 From: Felix F Zimmermann Date: Tue, 7 Jan 2025 23:48:43 +0100 Subject: [PATCH 17/22] Move functions from mixins into KData (#559) also fixes the issue with the KData.__init__ documentation --- src/mrpro/algorithms/prewhiten_kspace.py | 2 +- .../reconstruction/DirectReconstruction.py | 2 +- .../IterativeSENSEReconstruction.py | 2 +- .../reconstruction/Reconstruction.py | 2 +- ...RegularizedIterativeSENSEReconstruction.py | 2 +- src/mrpro/data/{_kdata => }/KData.py | 300 +++++++++++++++++- src/mrpro/data/__init__.py | 2 +- src/mrpro/data/_kdata/KDataProtocol.py | 41 --- src/mrpro/data/_kdata/KDataRearrangeMixin.py | 41 --- src/mrpro/data/_kdata/KDataRemoveOsMixin.py | 75 ----- src/mrpro/data/_kdata/KDataSelectMixin.py | 65 ---- src/mrpro/data/_kdata/KDataSplitMixin.py | 161 ---------- src/mrpro/operators/FourierOp.py | 2 +- 13 files changed, 300 insertions(+), 397 deletions(-) rename src/mrpro/data/{_kdata => }/KData.py (57%) delete mode 100644 src/mrpro/data/_kdata/KDataProtocol.py delete mode 100644 src/mrpro/data/_kdata/KDataRearrangeMixin.py delete mode 100644 src/mrpro/data/_kdata/KDataRemoveOsMixin.py delete mode 100644 src/mrpro/data/_kdata/KDataSelectMixin.py delete mode 100644 src/mrpro/data/_kdata/KDataSplitMixin.py diff --git a/src/mrpro/algorithms/prewhiten_kspace.py b/src/mrpro/algorithms/prewhiten_kspace.py index ab50f325e..3e1dde12d 100644 --- a/src/mrpro/algorithms/prewhiten_kspace.py +++ b/src/mrpro/algorithms/prewhiten_kspace.py @@ -5,7 +5,7 @@ import torch from einops import einsum, parse_shape, rearrange -from mrpro.data._kdata.KData import KData +from mrpro.data.KData import KData from mrpro.data.KNoise import KNoise diff --git a/src/mrpro/algorithms/reconstruction/DirectReconstruction.py b/src/mrpro/algorithms/reconstruction/DirectReconstruction.py index 3feab5cdb..a201b86c2 100644 --- a/src/mrpro/algorithms/reconstruction/DirectReconstruction.py +++ b/src/mrpro/algorithms/reconstruction/DirectReconstruction.py @@ -3,10 +3,10 @@ from collections.abc import Callable from mrpro.algorithms.reconstruction.Reconstruction import Reconstruction -from mrpro.data._kdata.KData import KData from mrpro.data.CsmData import CsmData from mrpro.data.DcfData import DcfData from mrpro.data.IData import IData +from mrpro.data.KData import KData from mrpro.data.KNoise import KNoise from mrpro.operators.FourierOp import FourierOp from mrpro.operators.LinearOperator import LinearOperator diff --git a/src/mrpro/algorithms/reconstruction/IterativeSENSEReconstruction.py b/src/mrpro/algorithms/reconstruction/IterativeSENSEReconstruction.py index 32785a91a..444d85712 100644 --- a/src/mrpro/algorithms/reconstruction/IterativeSENSEReconstruction.py +++ b/src/mrpro/algorithms/reconstruction/IterativeSENSEReconstruction.py @@ -7,9 +7,9 @@ from mrpro.algorithms.reconstruction.RegularizedIterativeSENSEReconstruction import ( RegularizedIterativeSENSEReconstruction, ) -from mrpro.data._kdata.KData import KData from mrpro.data.CsmData import CsmData from mrpro.data.DcfData import DcfData +from mrpro.data.KData import KData from mrpro.data.KNoise import KNoise from mrpro.operators.LinearOperator import LinearOperator diff --git a/src/mrpro/algorithms/reconstruction/Reconstruction.py b/src/mrpro/algorithms/reconstruction/Reconstruction.py index c4208157e..54a1f6af2 100644 --- a/src/mrpro/algorithms/reconstruction/Reconstruction.py +++ b/src/mrpro/algorithms/reconstruction/Reconstruction.py @@ -8,10 +8,10 @@ from typing_extensions import Self from mrpro.algorithms.prewhiten_kspace import prewhiten_kspace -from mrpro.data._kdata.KData import KData from mrpro.data.CsmData import CsmData from mrpro.data.DcfData import DcfData from mrpro.data.IData import IData +from mrpro.data.KData import KData from mrpro.data.KNoise import KNoise from mrpro.operators.FourierOp import FourierOp from mrpro.operators.LinearOperator import LinearOperator diff --git a/src/mrpro/algorithms/reconstruction/RegularizedIterativeSENSEReconstruction.py b/src/mrpro/algorithms/reconstruction/RegularizedIterativeSENSEReconstruction.py index c9a307ebe..e3c1c49ce 100644 --- a/src/mrpro/algorithms/reconstruction/RegularizedIterativeSENSEReconstruction.py +++ b/src/mrpro/algorithms/reconstruction/RegularizedIterativeSENSEReconstruction.py @@ -9,10 +9,10 @@ from mrpro.algorithms.optimizers.cg import cg from mrpro.algorithms.prewhiten_kspace import prewhiten_kspace from mrpro.algorithms.reconstruction.DirectReconstruction import DirectReconstruction -from mrpro.data._kdata.KData import KData from mrpro.data.CsmData import CsmData from mrpro.data.DcfData import DcfData from mrpro.data.IData import IData +from mrpro.data.KData import KData from mrpro.data.KNoise import KNoise from mrpro.operators.IdentityOp import IdentityOp from mrpro.operators.LinearOperator import LinearOperator diff --git a/src/mrpro/data/_kdata/KData.py b/src/mrpro/data/KData.py similarity index 57% rename from src/mrpro/data/_kdata/KData.py rename to src/mrpro/data/KData.py index 4b5df6250..47ddc27d2 100644 --- a/src/mrpro/data/_kdata/KData.py +++ b/src/mrpro/data/KData.py @@ -1,23 +1,21 @@ """MR raw data / k-space data class.""" +import copy import dataclasses import datetime import warnings from collections.abc import Callable, Sequence from pathlib import Path from types import EllipsisType +from typing import Literal, cast import h5py import ismrmrd import numpy as np import torch -from einops import rearrange -from typing_extensions import Self +from einops import rearrange, repeat +from typing_extensions import Self, TypeVar -from mrpro.data._kdata.KDataRearrangeMixin import KDataRearrangeMixin -from mrpro.data._kdata.KDataRemoveOsMixin import KDataRemoveOsMixin -from mrpro.data._kdata.KDataSelectMixin import KDataSelectMixin -from mrpro.data._kdata.KDataSplitMixin import KDataSplitMixin from mrpro.data.acq_filters import has_n_coils, is_image_acquisition from mrpro.data.AcqInfo import AcqInfo, rearrange_acq_info_fields from mrpro.data.EncodingLimits import Limits @@ -29,6 +27,8 @@ from mrpro.data.traj_calculators.KTrajectoryCalculator import KTrajectoryCalculator from mrpro.data.traj_calculators.KTrajectoryIsmrmrd import KTrajectoryIsmrmrd +RotationOrTensor = TypeVar('RotationOrTensor', bound=torch.Tensor | Rotation) + KDIM_SORT_LABELS = ( 'k1', 'k2', @@ -63,7 +63,9 @@ @dataclasses.dataclass(slots=True, frozen=True) -class KData(KDataSplitMixin, KDataRearrangeMixin, KDataSelectMixin, KDataRemoveOsMixin, MoveDataMixin): +class KData( + MoveDataMixin, +): """MR raw data / k-space data class.""" header: KHeader @@ -366,3 +368,287 @@ def compress_coils( ).permute(*np.argsort(permute_order)) return type(self)(self.header.clone(), kdata_coil_compressed, self.traj.clone()) + + def rearrange_k2_k1_into_k1(self: Self) -> Self: + """Rearrange kdata from (... k2 k1 ...) to (... 1 (k2 k1) ...). + + Parameters + ---------- + kdata + K-space data (other coils k2 k1 k0) + + Returns + ------- + K-space data (other coils 1 (k2 k1) k0) + """ + # Rearrange data + kdat = rearrange(self.data, '... coils k2 k1 k0->... coils 1 (k2 k1) k0') + + # Rearrange trajectory + ktraj = rearrange(self.traj.as_tensor(), 'dim ... k2 k1 k0-> dim ... 1 (k2 k1) k0') + + # Create new header with correct shape + kheader = copy.deepcopy(self.header) + + # Update shape of acquisition info index + kheader.acq_info.apply_( + lambda field: rearrange_acq_info_fields(field, 'other k2 k1 ... -> other 1 (k2 k1) ...') + ) + + return type(self)(kheader, kdat, type(self.traj).from_tensor(ktraj)) + + def remove_readout_os(self: Self) -> Self: + """Remove any oversampling along the readout (k0) direction [GAD]_. + + Returns a copy of the data. + + Parameters + ---------- + kdata + K-space data + + Returns + ------- + Copy of K-space data with oversampling removed. + + Raises + ------ + ValueError + If the recon matrix along x is larger than the encoding matrix along x. + + References + ---------- + .. [GAD] Gadgetron https://github.com/gadgetron/gadgetron-python + """ + from mrpro.operators.FastFourierOp import FastFourierOp + + # Ratio of k0/x between encoded and recon space + x_ratio = self.header.recon_matrix.x / self.header.encoding_matrix.x + if x_ratio == 1: + # If the encoded and recon space is the same we don't have to do anything + return self + elif x_ratio > 1: + raise ValueError('Recon matrix along x should be equal or larger than encoding matrix along x.') + + # Starting and end point of image after removing oversampling + start_cropped_readout = (self.header.encoding_matrix.x - self.header.recon_matrix.x) // 2 + end_cropped_readout = start_cropped_readout + self.header.recon_matrix.x + + def crop_readout(data_to_crop: torch.Tensor) -> torch.Tensor: + # returns a cropped copy + return data_to_crop[..., start_cropped_readout:end_cropped_readout].clone() + + # Transform to image space along readout, crop to reconstruction matrix size and transform back + fourier_k0_op = FastFourierOp(dim=(-1,)) + (cropped_data,) = fourier_k0_op(crop_readout(*fourier_k0_op.H(self.data))) + + # Adapt trajectory + ks = [self.traj.kz, self.traj.ky, self.traj.kx] + # only cropped ks that are not broadcasted/singleton along k0 + cropped_ks = [crop_readout(k) if k.shape[-1] > 1 else k.clone() for k in ks] + cropped_traj = KTrajectory(cropped_ks[0], cropped_ks[1], cropped_ks[2]) + + # Adapt header parameters + header = copy.deepcopy(self.header) + header.acq_info.center_sample -= start_cropped_readout + header.acq_info.number_of_samples[:] = cropped_data.shape[-1] + header.encoding_matrix.x = cropped_data.shape[-1] + + header.acq_info.discard_post = (header.acq_info.discard_post * x_ratio).to(torch.int32) + header.acq_info.discard_pre = (header.acq_info.discard_pre * x_ratio).to(torch.int32) + + return type(self)(header, cropped_data, cropped_traj) + + def select_other_subset( + self: Self, + subset_idx: torch.Tensor, + subset_label: Literal['average', 'slice', 'contrast', 'phase', 'repetition', 'set'], + ) -> Self: + """Select a subset from the other dimension of KData. + + Parameters + ---------- + kdata + K-space data (other coils k2 k1 k0) + subset_idx + Index which elements of the other subset to use, e.g. phase 0,1,2 and 5 + subset_label + Name of the other label, e.g. phase + + Returns + ------- + K-space data (other_subset coils k2 k1 k0) + + Raises + ------ + ValueError + If the subset indices are not available in the data + """ + # Make a copy such that the original kdata.header remains the same + kheader = copy.deepcopy(self.header) + ktraj = self.traj.as_tensor() + + # Verify that the subset_idx is available + label_idx = getattr(kheader.acq_info.idx, subset_label) + if not all(el in torch.unique(label_idx) for el in subset_idx): + raise ValueError('Subset indices are outside of the available index range') + + # Find subset index in acq_info index + other_idx = torch.cat([torch.where(idx == label_idx[:, 0, 0])[0] for idx in subset_idx], dim=0) + + # Adapt header + kheader.acq_info.apply_( + lambda field: field[other_idx, ...] if isinstance(field, torch.Tensor | Rotation) else field + ) + + # Select data + kdat = self.data[other_idx, ...] + + # Select ktraj + if ktraj.shape[1] > 1: + ktraj = ktraj[:, other_idx, ...] + + return type(self)(kheader, kdat, type(self.traj).from_tensor(ktraj)) + + def _split_k2_or_k1_into_other( + self, + split_idx: torch.Tensor, + other_label: Literal['average', 'slice', 'contrast', 'phase', 'repetition', 'set'], + split_dir: Literal['k2', 'k1'], + ) -> Self: + """Based on an index tensor, split the data in e.g. phases. + + Parameters + ---------- + split_idx + 2D index describing the k2 or k1 points in each block to be moved to the other dimension + (other_split, k1_per_split) or (other_split, k2_per_split) + other_label + Label of other dimension, e.g. repetition, phase + split_dir + Dimension to split, either 'k1' or 'k2' + + Returns + ------- + K-space data with new shape + ((other other_split) coils k2 k1_per_split k0) or ((other other_split) coils k2_per_split k1 k0) + + Raises + ------ + ValueError + Already existing "other_label" can only be of length 1 + """ + # Number of other + n_other = split_idx.shape[0] + + # Verify that the specified label of the other dimension is unused + if getattr(self.header.encoding_limits, other_label).length > 1: + raise ValueError(f'{other_label} is already used to encode different parts of the scan.') + + # Set-up splitting + if split_dir == 'k1': + # Split along k1 dimensions + def split_data_traj(dat_traj: torch.Tensor) -> torch.Tensor: + return dat_traj[:, :, :, split_idx, :] + + def split_acq_info(acq_info: RotationOrTensor) -> RotationOrTensor: + # cast due to https://github.com/python/mypy/issues/10817 + return cast(RotationOrTensor, acq_info[:, :, split_idx, ...]) + + # Rearrange other_split and k1 dimension + rearrange_pattern_data = 'other coils k2 other_split k1 k0->(other other_split) coils k2 k1 k0' + rearrange_pattern_traj = 'dim other k2 other_split k1 k0->dim (other other_split) k2 k1 k0' + rearrange_pattern_acq_info = 'other k2 other_split k1 ... -> (other other_split) k2 k1 ...' + + elif split_dir == 'k2': + # Split along k2 dimensions + def split_data_traj(dat_traj: torch.Tensor) -> torch.Tensor: + return dat_traj[:, :, split_idx, :, :] + + def split_acq_info(acq_info: RotationOrTensor) -> RotationOrTensor: + return cast(RotationOrTensor, acq_info[:, split_idx, ...]) + + # Rearrange other_split and k1 dimension + rearrange_pattern_data = 'other coils other_split k2 k1 k0->(other other_split) coils k2 k1 k0' + rearrange_pattern_traj = 'dim other other_split k2 k1 k0->dim (other other_split) k2 k1 k0' + rearrange_pattern_acq_info = 'other other_split k2 k1 ... -> (other other_split) k2 k1 ...' + + else: + raise ValueError('split_dir has to be "k1" or "k2"') + + # Split data + kdat = rearrange(split_data_traj(self.data), rearrange_pattern_data) + + # First we need to make sure the other dimension is the same as data then we can split the trajectory + ktraj = self.traj.as_tensor() + # Verify that other dimension of trajectory is 1 or matches data + if ktraj.shape[1] > 1 and ktraj.shape[1] != self.data.shape[0]: + raise ValueError(f'other dimension of trajectory has to be 1 or match data ({self.data.shape[0]})') + elif ktraj.shape[1] == 1 and self.data.shape[0] > 1: + ktraj = repeat(ktraj, 'dim other k2 k1 k0->dim (other_data other) k2 k1 k0', other_data=self.data.shape[0]) + ktraj = rearrange(split_data_traj(ktraj), rearrange_pattern_traj) + + # Create new header with correct shape + kheader = self.header.clone() + + # Update shape of acquisition info index + kheader.acq_info.apply_( + lambda field: rearrange_acq_info_fields(split_acq_info(field), rearrange_pattern_acq_info) + if isinstance(field, Rotation | torch.Tensor) + else field + ) + + # Update other label limits and acquisition info + setattr(kheader.encoding_limits, other_label, Limits(min=0, max=n_other - 1, center=0)) + + # acq_info for new other dimensions + acq_info_other_split = repeat( + torch.linspace(0, n_other - 1, n_other), 'other-> other k2 k1', k2=kdat.shape[-3], k1=kdat.shape[-2] + ) + setattr(kheader.acq_info.idx, other_label, acq_info_other_split) + + return type(self)(kheader, kdat, type(self.traj).from_tensor(ktraj)) + + def split_k1_into_other( + self: Self, + split_idx: torch.Tensor, + other_label: Literal['average', 'slice', 'contrast', 'phase', 'repetition', 'set'], + ) -> Self: + """Based on an index tensor, split the data in e.g. phases. + + Parameters + ---------- + kdata + K-space data (other coils k2 k1 k0) + split_idx + 2D index describing the k1 points in each block to be moved to other dimension (other_split, k1_per_split) + other_label + Label of other dimension, e.g. repetition, phase + + Returns + ------- + K-space data with new shape ((other other_split) coils k2 k1_per_split k0) + """ + return self._split_k2_or_k1_into_other(split_idx, other_label, split_dir='k1') + + def split_k2_into_other( + self: Self, + split_idx: torch.Tensor, + other_label: Literal['average', 'slice', 'contrast', 'phase', 'repetition', 'set'], + ) -> Self: + """Based on an index tensor, split the data in e.g. phases. + + Parameters + ---------- + kdata + K-space data (other coils k2 k1 k0) + split_idx + 2D index describing the k2 points in each block to be moved to other dimension (other_split, k2_per_split) + other_label + Label of other dimension, e.g. repetition, phase + + Returns + ------- + K-space data with new shape ((other other_split) coils k2_per_split k1 k0) + """ + return self._split_k2_or_k1_into_other(split_idx, other_label, split_dir='k2') diff --git a/src/mrpro/data/__init__.py b/src/mrpro/data/__init__.py index d5667a5bc..89cc6695b 100644 --- a/src/mrpro/data/__init__.py +++ b/src/mrpro/data/__init__.py @@ -6,7 +6,7 @@ from mrpro.data.EncodingLimits import EncodingLimits, Limits from mrpro.data.IData import IData from mrpro.data.IHeader import IHeader -from mrpro.data._kdata.KData import KData +from mrpro.data.KData import KData from mrpro.data.KHeader import KHeader from mrpro.data.KNoise import KNoise from mrpro.data.KTrajectory import KTrajectory diff --git a/src/mrpro/data/_kdata/KDataProtocol.py b/src/mrpro/data/_kdata/KDataProtocol.py deleted file mode 100644 index 485a8fc4d..000000000 --- a/src/mrpro/data/_kdata/KDataProtocol.py +++ /dev/null @@ -1,41 +0,0 @@ -"""Protocol for KData.""" - -from typing import Literal - -import torch -from typing_extensions import Protocol, Self - -from mrpro.data.KHeader import KHeader -from mrpro.data.KTrajectory import KTrajectory - - -class _KDataProtocol(Protocol): - """Protocol for KData used for type hinting in KData mixins. - - Note that the actual KData class can have more properties and methods than those defined here. - - If you want to use a property or method of KData in a new KDataMixin class, - you must add it to this Protocol to make sure that the type hinting works [PRO]_. - - References - ---------- - .. [PRO] Protocols https://typing.readthedocs.io/en/latest/spec/protocol.html#protocols - """ - - @property - def header(self) -> KHeader: ... - - @property - def data(self) -> torch.Tensor: ... - - @property - def traj(self) -> KTrajectory: ... - - def __init__(self, header: KHeader, data: torch.Tensor, traj: KTrajectory): ... - - def _split_k2_or_k1_into_other( - self, - split_idx: torch.Tensor, - other_label: Literal['average', 'slice', 'contrast', 'phase', 'repetition', 'set'], - split_dir: Literal['k1', 'k2'], - ) -> Self: ... diff --git a/src/mrpro/data/_kdata/KDataRearrangeMixin.py b/src/mrpro/data/_kdata/KDataRearrangeMixin.py deleted file mode 100644 index 23a58dea6..000000000 --- a/src/mrpro/data/_kdata/KDataRearrangeMixin.py +++ /dev/null @@ -1,41 +0,0 @@ -"""Rearrange KData.""" - -import copy - -from einops import rearrange -from typing_extensions import Self - -from mrpro.data._kdata.KDataProtocol import _KDataProtocol -from mrpro.data.AcqInfo import rearrange_acq_info_fields - - -class KDataRearrangeMixin(_KDataProtocol): - """Rearrange KData.""" - - def rearrange_k2_k1_into_k1(self: Self) -> Self: - """Rearrange kdata from (... k2 k1 ...) to (... 1 (k2 k1) ...). - - Parameters - ---------- - kdata - K-space data (other coils k2 k1 k0) - - Returns - ------- - K-space data (other coils 1 (k2 k1) k0) - """ - # Rearrange data - kdat = rearrange(self.data, '... coils k2 k1 k0->... coils 1 (k2 k1) k0') - - # Rearrange trajectory - ktraj = rearrange(self.traj.as_tensor(), 'dim ... k2 k1 k0-> dim ... 1 (k2 k1) k0') - - # Create new header with correct shape - kheader = copy.deepcopy(self.header) - - # Update shape of acquisition info index - kheader.acq_info.apply_( - lambda field: rearrange_acq_info_fields(field, 'other k2 k1 ... -> other 1 (k2 k1) ...') - ) - - return type(self)(kheader, kdat, type(self.traj).from_tensor(ktraj)) diff --git a/src/mrpro/data/_kdata/KDataRemoveOsMixin.py b/src/mrpro/data/_kdata/KDataRemoveOsMixin.py deleted file mode 100644 index 555f56a39..000000000 --- a/src/mrpro/data/_kdata/KDataRemoveOsMixin.py +++ /dev/null @@ -1,75 +0,0 @@ -"""Remove oversampling along readout dimension.""" - -from copy import deepcopy - -import torch -from typing_extensions import Self - -from mrpro.data._kdata.KDataProtocol import _KDataProtocol -from mrpro.data.KTrajectory import KTrajectory - - -class KDataRemoveOsMixin(_KDataProtocol): - """Remove oversampling along readout dimension.""" - - def remove_readout_os(self: Self) -> Self: - """Remove any oversampling along the readout (k0) direction [GAD]_. - - Returns a copy of the data. - - Parameters - ---------- - kdata - K-space data - - Returns - ------- - Copy of K-space data with oversampling removed. - - Raises - ------ - ValueError - If the recon matrix along x is larger than the encoding matrix along x. - - References - ---------- - .. [GAD] Gadgetron https://github.com/gadgetron/gadgetron-python - """ - from mrpro.operators.FastFourierOp import FastFourierOp - - # Ratio of k0/x between encoded and recon space - x_ratio = self.header.recon_matrix.x / self.header.encoding_matrix.x - if x_ratio == 1: - # If the encoded and recon space is the same we don't have to do anything - return self - elif x_ratio > 1: - raise ValueError('Recon matrix along x should be equal or larger than encoding matrix along x.') - - # Starting and end point of image after removing oversampling - start_cropped_readout = (self.header.encoding_matrix.x - self.header.recon_matrix.x) // 2 - end_cropped_readout = start_cropped_readout + self.header.recon_matrix.x - - def crop_readout(data_to_crop: torch.Tensor) -> torch.Tensor: - # returns a cropped copy - return data_to_crop[..., start_cropped_readout:end_cropped_readout].clone() - - # Transform to image space along readout, crop to reconstruction matrix size and transform back - fourier_k0_op = FastFourierOp(dim=(-1,)) - (cropped_data,) = fourier_k0_op(crop_readout(*fourier_k0_op.H(self.data))) - - # Adapt trajectory - ks = [self.traj.kz, self.traj.ky, self.traj.kx] - # only cropped ks that are not broadcasted/singleton along k0 - cropped_ks = [crop_readout(k) if k.shape[-1] > 1 else k.clone() for k in ks] - cropped_traj = KTrajectory(cropped_ks[0], cropped_ks[1], cropped_ks[2]) - - # Adapt header parameters - header = deepcopy(self.header) - header.acq_info.center_sample -= start_cropped_readout - header.acq_info.number_of_samples[:] = cropped_data.shape[-1] - header.encoding_matrix.x = cropped_data.shape[-1] - - header.acq_info.discard_post = (header.acq_info.discard_post * x_ratio).to(torch.int32) - header.acq_info.discard_pre = (header.acq_info.discard_pre * x_ratio).to(torch.int32) - - return type(self)(header, cropped_data, cropped_traj) diff --git a/src/mrpro/data/_kdata/KDataSelectMixin.py b/src/mrpro/data/_kdata/KDataSelectMixin.py deleted file mode 100644 index 8f8a452cf..000000000 --- a/src/mrpro/data/_kdata/KDataSelectMixin.py +++ /dev/null @@ -1,65 +0,0 @@ -"""Select subset along other dimensions of KData.""" - -import copy -from typing import Literal - -import torch -from typing_extensions import Self - -from mrpro.data._kdata.KDataProtocol import _KDataProtocol -from mrpro.data.Rotation import Rotation - - -class KDataSelectMixin(_KDataProtocol): - """Select subset of KData.""" - - def select_other_subset( - self: Self, - subset_idx: torch.Tensor, - subset_label: Literal['average', 'slice', 'contrast', 'phase', 'repetition', 'set'], - ) -> Self: - """Select a subset from the other dimension of KData. - - Parameters - ---------- - kdata - K-space data (other coils k2 k1 k0) - subset_idx - Index which elements of the other subset to use, e.g. phase 0,1,2 and 5 - subset_label - Name of the other label, e.g. phase - - Returns - ------- - K-space data (other_subset coils k2 k1 k0) - - Raises - ------ - ValueError - If the subset indices are not available in the data - """ - # Make a copy such that the original kdata.header remains the same - kheader = copy.deepcopy(self.header) - ktraj = self.traj.as_tensor() - - # Verify that the subset_idx is available - label_idx = getattr(kheader.acq_info.idx, subset_label) - if not all(el in torch.unique(label_idx) for el in subset_idx): - raise ValueError('Subset indices are outside of the available index range') - - # Find subset index in acq_info index - other_idx = torch.cat([torch.where(idx == label_idx[:, 0, 0])[0] for idx in subset_idx], dim=0) - - # Adapt header - kheader.acq_info.apply_( - lambda field: field[other_idx, ...] if isinstance(field, torch.Tensor | Rotation) else field - ) - - # Select data - kdat = self.data[other_idx, ...] - - # Select ktraj - if ktraj.shape[1] > 1: - ktraj = ktraj[:, other_idx, ...] - - return type(self)(kheader, kdat, type(self.traj).from_tensor(ktraj)) diff --git a/src/mrpro/data/_kdata/KDataSplitMixin.py b/src/mrpro/data/_kdata/KDataSplitMixin.py deleted file mode 100644 index c28004af4..000000000 --- a/src/mrpro/data/_kdata/KDataSplitMixin.py +++ /dev/null @@ -1,161 +0,0 @@ -"""Mixin class to split KData into other subsets.""" - -from typing import Literal, TypeVar, cast - -import torch -from einops import rearrange, repeat -from typing_extensions import Self - -from mrpro.data._kdata.KDataProtocol import _KDataProtocol -from mrpro.data.AcqInfo import rearrange_acq_info_fields -from mrpro.data.EncodingLimits import Limits -from mrpro.data.Rotation import Rotation - -RotationOrTensor = TypeVar('RotationOrTensor', bound=torch.Tensor | Rotation) - - -class KDataSplitMixin(_KDataProtocol): - """Split KData into other subsets.""" - - def _split_k2_or_k1_into_other( - self, - split_idx: torch.Tensor, - other_label: Literal['average', 'slice', 'contrast', 'phase', 'repetition', 'set'], - split_dir: Literal['k2', 'k1'], - ) -> Self: - """Based on an index tensor, split the data in e.g. phases. - - Parameters - ---------- - split_idx - 2D index describing the k2 or k1 points in each block to be moved to the other dimension - (other_split, k1_per_split) or (other_split, k2_per_split) - other_label - Label of other dimension, e.g. repetition, phase - split_dir - Dimension to split, either 'k1' or 'k2' - - Returns - ------- - K-space data with new shape - ((other other_split) coils k2 k1_per_split k0) or ((other other_split) coils k2_per_split k1 k0) - - Raises - ------ - ValueError - Already existing "other_label" can only be of length 1 - """ - # Number of other - n_other = split_idx.shape[0] - - # Verify that the specified label of the other dimension is unused - if getattr(self.header.encoding_limits, other_label).length > 1: - raise ValueError(f'{other_label} is already used to encode different parts of the scan.') - - # Set-up splitting - if split_dir == 'k1': - # Split along k1 dimensions - def split_data_traj(dat_traj: torch.Tensor) -> torch.Tensor: - return dat_traj[:, :, :, split_idx, :] - - def split_acq_info(acq_info: RotationOrTensor) -> RotationOrTensor: - # cast due to https://github.com/python/mypy/issues/10817 - return cast(RotationOrTensor, acq_info[:, :, split_idx, ...]) - - # Rearrange other_split and k1 dimension - rearrange_pattern_data = 'other coils k2 other_split k1 k0->(other other_split) coils k2 k1 k0' - rearrange_pattern_traj = 'dim other k2 other_split k1 k0->dim (other other_split) k2 k1 k0' - rearrange_pattern_acq_info = 'other k2 other_split k1 ... -> (other other_split) k2 k1 ...' - - elif split_dir == 'k2': - # Split along k2 dimensions - def split_data_traj(dat_traj: torch.Tensor) -> torch.Tensor: - return dat_traj[:, :, split_idx, :, :] - - def split_acq_info(acq_info: RotationOrTensor) -> RotationOrTensor: - return cast(RotationOrTensor, acq_info[:, split_idx, ...]) - - # Rearrange other_split and k1 dimension - rearrange_pattern_data = 'other coils other_split k2 k1 k0->(other other_split) coils k2 k1 k0' - rearrange_pattern_traj = 'dim other other_split k2 k1 k0->dim (other other_split) k2 k1 k0' - rearrange_pattern_acq_info = 'other other_split k2 k1 ... -> (other other_split) k2 k1 ...' - - else: - raise ValueError('split_dir has to be "k1" or "k2"') - - # Split data - kdat = rearrange(split_data_traj(self.data), rearrange_pattern_data) - - # First we need to make sure the other dimension is the same as data then we can split the trajectory - ktraj = self.traj.as_tensor() - # Verify that other dimension of trajectory is 1 or matches data - if ktraj.shape[1] > 1 and ktraj.shape[1] != self.data.shape[0]: - raise ValueError(f'other dimension of trajectory has to be 1 or match data ({self.data.shape[0]})') - elif ktraj.shape[1] == 1 and self.data.shape[0] > 1: - ktraj = repeat(ktraj, 'dim other k2 k1 k0->dim (other_data other) k2 k1 k0', other_data=self.data.shape[0]) - ktraj = rearrange(split_data_traj(ktraj), rearrange_pattern_traj) - - # Create new header with correct shape - kheader = self.header.clone() - - # Update shape of acquisition info index - kheader.acq_info.apply_( - lambda field: rearrange_acq_info_fields(split_acq_info(field), rearrange_pattern_acq_info) - if isinstance(field, Rotation | torch.Tensor) - else field - ) - - # Update other label limits and acquisition info - setattr(kheader.encoding_limits, other_label, Limits(min=0, max=n_other - 1, center=0)) - - # acq_info for new other dimensions - acq_info_other_split = repeat( - torch.linspace(0, n_other - 1, n_other), 'other-> other k2 k1', k2=kdat.shape[-3], k1=kdat.shape[-2] - ) - setattr(kheader.acq_info.idx, other_label, acq_info_other_split) - - return type(self)(kheader, kdat, type(self.traj).from_tensor(ktraj)) - - def split_k1_into_other( - self: Self, - split_idx: torch.Tensor, - other_label: Literal['average', 'slice', 'contrast', 'phase', 'repetition', 'set'], - ) -> Self: - """Based on an index tensor, split the data in e.g. phases. - - Parameters - ---------- - kdata - K-space data (other coils k2 k1 k0) - split_idx - 2D index describing the k1 points in each block to be moved to other dimension (other_split, k1_per_split) - other_label - Label of other dimension, e.g. repetition, phase - - Returns - ------- - K-space data with new shape ((other other_split) coils k2 k1_per_split k0) - """ - return self._split_k2_or_k1_into_other(split_idx, other_label, split_dir='k1') - - def split_k2_into_other( - self: Self, - split_idx: torch.Tensor, - other_label: Literal['average', 'slice', 'contrast', 'phase', 'repetition', 'set'], - ) -> Self: - """Based on an index tensor, split the data in e.g. phases. - - Parameters - ---------- - kdata - K-space data (other coils k2 k1 k0) - split_idx - 2D index describing the k2 points in each block to be moved to other dimension (other_split, k2_per_split) - other_label - Label of other dimension, e.g. repetition, phase - - Returns - ------- - K-space data with new shape ((other other_split) coils k2_per_split k1 k0) - """ - return self._split_k2_or_k1_into_other(split_idx, other_label, split_dir='k2') diff --git a/src/mrpro/operators/FourierOp.py b/src/mrpro/operators/FourierOp.py index 2e2087f78..c32adb71d 100644 --- a/src/mrpro/operators/FourierOp.py +++ b/src/mrpro/operators/FourierOp.py @@ -8,8 +8,8 @@ from torchkbnufft import KbNufft, KbNufftAdjoint from typing_extensions import Self -from mrpro.data._kdata.KData import KData from mrpro.data.enums import TrajType +from mrpro.data.KData import KData from mrpro.data.KTrajectory import KTrajectory from mrpro.data.SpatialDimension import SpatialDimension from mrpro.operators.CartesianSamplingOp import CartesianSamplingOp From cb9629cc7620b173c8263a101da20c6aa4f84ccc Mon Sep 17 00:00:00 2001 From: Felix F Zimmermann Date: Tue, 7 Jan 2025 23:57:46 +0100 Subject: [PATCH 18/22] Ensure compressed number of coils is not greater than existing number (#567) --- src/mrpro/data/KData.py | 26 ++++++++++++++++---------- tests/data/test_kdata.py | 7 +++++++ 2 files changed, 23 insertions(+), 10 deletions(-) diff --git a/src/mrpro/data/KData.py b/src/mrpro/data/KData.py index 47ddc27d2..719141921 100644 --- a/src/mrpro/data/KData.py +++ b/src/mrpro/data/KData.py @@ -332,6 +332,13 @@ def compress_coils( from mrpro.operators import PCACompressionOp coil_dim = -4 % self.data.ndim + + if n_compressed_coils > (n_current_coils := self.data.shape[coil_dim]): + raise ValueError( + f'Number of compressed coils ({n_compressed_coils}) cannot be greater ' + f'than the number of current coils ({n_current_coils}).' + ) + if batch_dims is not None and joint_dims is not Ellipsis: raise ValueError('Either batch_dims or joint_dims can be defined not both.') @@ -349,22 +356,21 @@ def compress_coils( # reshape to (*batch dimension, -1, coils) permute_order = ( - batch_dims_normalized - + [i for i in range(self.data.ndim) if i != coil_dim and i not in batch_dims_normalized] - + [coil_dim] + *batch_dims_normalized, + *[i for i in range(self.data.ndim) if i != coil_dim and i not in batch_dims_normalized], + coil_dim, ) - kdata_coil_compressed = self.data.permute(permute_order) - permuted_kdata_shape = kdata_coil_compressed.shape - kdata_coil_compressed = kdata_coil_compressed.flatten( + kdata_permuted = self.data.permute(permute_order) + kdata_flattened = kdata_permuted.flatten( start_dim=len(batch_dims_normalized), end_dim=-2 ) # keep separate dimensions and coil - pca_compression_op = PCACompressionOp(data=kdata_coil_compressed, n_components=n_compressed_coils) - (kdata_coil_compressed,) = pca_compression_op(kdata_coil_compressed) - + pca_compression_op = PCACompressionOp(data=kdata_flattened, n_components=n_compressed_coils) + (kdata_coil_compressed_flattened,) = pca_compression_op(kdata_flattened) + del kdata_flattened # reshape to original dimensions and undo permutation kdata_coil_compressed = torch.reshape( - kdata_coil_compressed, [*permuted_kdata_shape[:-1], n_compressed_coils] + kdata_coil_compressed_flattened, [*kdata_permuted.shape[:-1], n_compressed_coils] ).permute(*np.argsort(permute_order)) return type(self)(self.header.clone(), kdata_coil_compressed, self.traj.clone()) diff --git a/tests/data/test_kdata.py b/tests/data/test_kdata.py index 42cd061f6..158b7ede6 100644 --- a/tests/data/test_kdata.py +++ b/tests/data/test_kdata.py @@ -574,3 +574,10 @@ def test_KData_compress_coils_error_coil_dim(consistently_shaped_kdata): with pytest.raises(ValueError, match='Coil dimension must not'): consistently_shaped_kdata.compress_coils(n_compressed_coils=3, joint_dims=(-4,)) + + +def test_KData_compress_coils_error_n_coils(consistently_shaped_kdata): + """Test if error is raised if new coils would be larger than existing coils""" + existing_coils = consistently_shaped_kdata.data.shape[-4] + with pytest.raises(ValueError, match='greater'): + consistently_shaped_kdata.compress_coils(existing_coils + 1) From f040c001d68985c25674724794743e0f53764903 Mon Sep 17 00:00:00 2001 From: Felix F Zimmermann Date: Wed, 8 Jan 2025 00:07:27 +0100 Subject: [PATCH 19/22] Release v0.250107 (#599) --- src/mrpro/VERSION | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/mrpro/VERSION b/src/mrpro/VERSION index 8c5b741c6..5c4cfd72c 100644 --- a/src/mrpro/VERSION +++ b/src/mrpro/VERSION @@ -1 +1 @@ -0.241217 +0.250107 From a1873a0ae0a1db04d8549d426ffb448e515245b3 Mon Sep 17 00:00:00 2001 From: Lunin Leonid Date: Wed, 8 Jan 2025 09:54:12 +0100 Subject: [PATCH 20/22] Add dark and light logo to README (#600) --- README.md | 8 +++++++- 1 file changed, 7 insertions(+), 1 deletion(-) diff --git a/README.md b/README.md index 4dcce2430..2e53443b2 100644 --- a/README.md +++ b/README.md @@ -1,5 +1,11 @@

-MRpro logo + + + + + MRpro logo + +


![Python](https://img.shields.io/badge/python-3.10%20%7C%203.11%20%7C%203.12-blue) From a9724e9a144d1ce53e51fdfc6b76963f6c7d8c19 Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Thu, 9 Jan 2025 21:25:25 +0000 Subject: [PATCH 21/22] [pre-commit] pre-commit autoupdate (#597) Co-authored-by: pre-commit-ci[bot] <66853113+pre-commit-ci[bot]@users.noreply.github.com> Co-authored-by: Felix F Zimmermann --- .pre-commit-config.yaml | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml index 3920fb5ae..1dc87ab60 100644 --- a/.pre-commit-config.yaml +++ b/.pre-commit-config.yaml @@ -14,14 +14,14 @@ repos: - id: mixed-line-ending - repo: https://github.com/astral-sh/ruff-pre-commit - rev: v0.8.1 + rev: v0.8.6 hooks: - id: ruff # linter args: [--fix] - id: ruff-format # formatter - repo: https://github.com/crate-ci/typos - rev: typos-dict-v0.11.37 + rev: v1.29.4 hooks: - id: typos @@ -33,7 +33,7 @@ repos: exclude: ^tests/ - repo: https://github.com/pre-commit/mirrors-mypy - rev: v1.13.0 + rev: v1.14.1 hooks: - id: mypy pass_filenames: false From e029d5ede0b0a7c92d485f6cd7670939c9ab5f47 Mon Sep 17 00:00:00 2001 From: Felix F Zimmermann Date: Fri, 10 Jan 2025 07:57:20 +0100 Subject: [PATCH 22/22] Fix license in citation for zenodo (#601) --- CITATION.cff | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/CITATION.cff b/CITATION.cff index 372724793..825faa7c7 100644 --- a/CITATION.cff +++ b/CITATION.cff @@ -47,8 +47,7 @@ authors: cff-version: 1.2.0 date-released: "2024-12-17" doi: 10.5281/zenodo.14509598 -license: - - apache-2.0 +license: Apache-2.0 repository-code: https://github.com/PTB-MR/mrpro/ title: MRpro - PyTorch-based MR image reconstruction and processing package type: software