diff --git a/docs/source/conf.py b/docs/source/conf.py index 71f22e96..8a8538d9 100644 --- a/docs/source/conf.py +++ b/docs/source/conf.py @@ -1,14 +1,4 @@ # Configuration file for the Sphinx documentation builder. -# -# For the full list of built-in configuration values, see the documentation: -# https://www.sphinx-doc.org/en/master/usage/configuration.html - -# -- Path setup -------------------------------------------------------------- - -# If extensions (or modules to document with autodoc) are in another directory, -# add these directories to sys.path here. If the directory is relative to the -# documentation root, use os.path.abspath to make it absolute, like shown here. - import ast import dataclasses import inspect @@ -16,18 +6,18 @@ import shutil import sys from pathlib import Path +from typing import get_overloads +from sphinx.ext.autodoc import AttributeDocumenter, ClassDocumenter, MethodDocumenter, PropertyDocumenter +from sphinx.util.inspect import isclassmethod, isstaticmethod, signature, stringify_signature from sphinx_pyproject import SphinxConfig -from sphinx_pyproject import SphinxConfig -from sphinx.ext.autodoc import ClassDocumenter, MethodDocumenter, AttributeDocumenter, PropertyDocumenter -from sphinx.util.inspect import isstaticmethod, isclassmethod - from mrpro import __version__ as project_version -config = SphinxConfig('../../pyproject.toml', globalns=globals(), config_overrides={'version': project_version}) sys.path.insert(0, os.path.abspath('../../src')) # Source code dir relative to this file +config = SphinxConfig('../../pyproject.toml', globalns=globals(), config_overrides={'version': project_version}) + project = name copyright = '2023, Physikalisch-Technische Bundesanstalt (PTB) Berlin' @@ -81,6 +71,8 @@ 'dollarmath', ] nb_execution_mode = 'off' +nb_output_stderr ='remove' +nb_output_stdout = 'remove' html_theme = 'sphinx_rtd_theme' html_title = name html_show_sphinx = False @@ -97,8 +89,7 @@ 'github_user': 'PTB-MR', 'github_repo': 'mrpro', 'github_version': 'main', - 'github_url' : 'https://github.com/PTB-MR/mrpro/main' - + 'github_url': 'https://github.com/PTB-MR/mrpro/main', } linkcode_blob = html_context['github_version'] @@ -110,6 +101,7 @@ def get_lambda_source(obj): if isinstance(node, ast.Lambda): return ast.unparse(node.body) + class DefaultValue: """Used to store default values of dataclass fields with default factory.""" @@ -142,36 +134,55 @@ def rewrite_dataclass_init_default_factories(app, obj, bound_method) -> None: defaults[field.name] = DefaultValue(get_lambda_source(field.default_factory)) else: defaults[field.name] = DefaultValue(field.default_factory.__name__ + '()') - new_defaults = tuple(defaults.get(name, param.default) for name, param in parameters.items() if param.default != inspect._empty) + new_defaults = tuple( + defaults.get(name, param.default) for name, param in parameters.items() if param.default != inspect._empty + ) obj.__defaults__ = new_defaults +def autodoc_inherit_overload(app, what, name, obj, options, sig, ret_ann): + """Create overloaded signatures.""" + if what in ('function', 'method') and callable(obj): + try: + overloads = get_overloads(obj) + except: + return (sig, ret_ann) + if overloads: + kwargs = {} + if app.config.autodoc_typehints in ('none', 'description'): + kwargs['show_annotation'] = False + if app.config.autodoc_typehints_format == 'short': + kwargs['unqualified_typehints'] = True + type_aliases = app.config.autodoc_type_aliases + bound_method = what == 'method' + sigs = [] + for overload in overloads: + if hasattr(overload, '__func__'): + overload = overload.__func__ # classmethod or staticmethod + overload_sig = signature(overload, bound_method=bound_method, type_aliases=type_aliases) + sigs.append(stringify_signature(overload_sig, **kwargs)) + return '\n'.join(sigs), None + + class CustomClassDocumenter(ClassDocumenter): - """ - Custom Documenter to reorder class members - """ - - def sort_members( - self, documenters: list[tuple['Documenter', bool]], order: str - ) -> list[tuple['Documenter', bool]]: - """ - Sort the given member list with custom logic for `groupwise` ordering. - """ - if order == "groupwise": + """Custom Documenter to reorder class members.""" + + def sort_members(self, documenters: list[tuple['Documenter', bool]], order: str) -> list[tuple['Documenter', bool]]: + """Sort the given member list with custom logic for `groupwise` ordering.""" + if order == 'groupwise': if not self.parse_name() or not self.import_object(): return documenters # Split members into groups (non-inherited,inherited) - static_methods = [],[] - class_methods = [],[] - special_methods = [],[] - instance_methods = [],[] - attributes = [],[] - properties = [],[] - other=[],[] + static_methods = [], [] + class_methods = [], [] + special_methods = [], [] + instance_methods = [], [] + attributes = [], [] + properties = [], [] + other = [], [] others_methods = [] init_method = [] - for documenter in documenters: doc = documenter[0] parsed = doc.parse_name() and doc.import_object() @@ -180,15 +191,15 @@ def sort_members( attributes[inherited].append(documenter) elif isinstance(doc, PropertyDocumenter): properties[inherited].append(documenter) - elif isinstance(doc,MethodDocumenter): + elif isinstance(doc, MethodDocumenter): if not parsed: others_methods.append(documenter) continue - if doc.object_name == "__init__": + if doc.object_name == '__init__': init_method.append(documenter) - elif dataclasses.is_dataclass(self.object) and doc.object_name=="__new__": + elif dataclasses.is_dataclass(self.object) and doc.object_name == '__new__': ... - elif doc.object_name[:2]=="__": + elif doc.object_name[:2] == '__': special_methods[inherited].append(documenter) elif isclassmethod(doc.object): class_methods[inherited].append(documenter) @@ -201,18 +212,31 @@ def sort_members( continue # Combine groups in the desired order constructors = init_method + class_methods[0] + class_methods[1] - methods = instance_methods[0] + instance_methods[1] + others_methods + static_methods[0] + static_methods[1] + special_methods[0] + special_methods[1] - return constructors+ attributes[0]+attributes[1] + properties[0]+properties[1]+methods + other[0]+other[1] + methods = ( + instance_methods[0] + + instance_methods[1] + + others_methods + + static_methods[0] + + static_methods[1] + + special_methods[0] + + special_methods[1] + ) + return ( + constructors + + attributes[0] + + attributes[1] + + properties[0] + + properties[1] + + methods + + other[0] + + other[1] + ) else: return super().sort_members(documenters, order) - - def sync_notebooks(source_folder, dest_folder): - """ - Synchronize files from the source to the destination folder, copying only new or updated files. - """ + """Synchronize files from the source to the destination folder, copying only new or updated files.""" dest = Path(dest_folder) dest.mkdir(parents=True, exist_ok=True) for src_file in Path(source_folder).iterdir(): @@ -221,9 +245,10 @@ def sync_notebooks(source_folder, dest_folder): if not dest_file.exists() or src_file.stat().st_mtime > dest_file.stat().st_mtime: shutil.copy2(src_file, dest_file) + def setup(app): - app.set_html_assets_policy('always') # forces mathjax on all pages + app.set_html_assets_policy('always') # forces mathjax on all pages app.connect('autodoc-before-process-signature', rewrite_dataclass_init_default_factories) + app.connect('autodoc-process-signature', autodoc_inherit_overload, 0) app.add_autodocumenter(CustomClassDocumenter) - sync_notebooks(app.srcdir.parent.parent/'examples'/'notebooks', app.srcdir/'_notebooks') - + sync_notebooks(app.srcdir.parent.parent / 'examples' / 'notebooks', app.srcdir / '_notebooks') diff --git a/examples/cartesian_reconstruction.ipynb b/examples/cartesian_reconstruction.ipynb deleted file mode 100644 index 7134e1aa..00000000 --- a/examples/cartesian_reconstruction.ipynb +++ /dev/null @@ -1,719 +0,0 @@ -{ - "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": { - "lines_to_next_cell": 2 - }, - "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": "code", - "execution_count": null, - "id": "2f28c8bb", - "metadata": { - "tags": [ - "hide-cell" - ] - }, - "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\n", - "\n", - "print('Downloaded files:')\n", - "for f in data_folder.iterdir():\n", - " print(f.name)" - ] - }, - { - "cell_type": "markdown", - "id": "c8e96446", - "metadata": { - "lines_to_next_cell": 2 - }, - "source": [ - "We have three different scans obtained from the same object with the same FOV and resolution, saved as ISMRMRD\n", - "raw data files (*.mrd):\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": "9eab0f40", - "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 {py:class}`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.\n", - "Simply calling ``print(kdata)`` gives us a basic overview of the `KData` object." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "703daa5a", - "metadata": { - "lines_to_next_cell": 2 - }, - "outputs": [], - "source": [ - "print(kdata)" - ] - }, - { - "cell_type": "markdown", - "id": "016eb8d0", - "metadata": {}, - "source": [ - "We can also have a look at more specific header information like the 1H Lamor frequency" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "960a2c8a", - "metadata": {}, - "outputs": [], - "source": [ - "print('Lamor Frequency:', 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('Shape:', 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('Shape:', 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 = img.abs().square().sum(dim=-4).sqrt().squeeze()\n", - "\n", - "# plot the image" - ] - }, - { - "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": { - "inputHidden": true, - "lines_to_next_cell": 0, - "outputHidden": true, - "tags": [ - "remove-output" - ] - }, - "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)\n", - "plt.show()" - ] - }, - { - "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": {}, - "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": { - "lines_to_next_cell": 0 - }, - "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", - ")\n", - "plt.show()" - ] - }, - { - "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": "inputHidden,outputHidden,tags,-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/scripts/cartesian_reconstruction.py similarity index 89% rename from examples/cartesian_reconstruction.py rename to examples/scripts/cartesian_reconstruction.py index c02b0270..461acf84 100644 --- a/examples/cartesian_reconstruction.py +++ b/examples/scripts/cartesian_reconstruction.py @@ -2,10 +2,8 @@ # # 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 @@ -23,31 +21,27 @@ dataset = '14173489' zenodo_get.zenodo_get([dataset, '-r', 5, '-o', data_folder]) # r: retries -print('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, saved as ISMRMRD -# raw data files (*.mrd): +# raw data files (*.mrd or *.h5): # -# - cart_t1.mrd is a fully sampled Cartesian acquisition +# - ``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_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 +# - ``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 {py:class}`KData` object. +# To read in an ISMRMRD file, we can simply pass on the file name to a {py:class}`~mrpro.data.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. +# calculator (called {py:class}`~mrpro.data.traj_calculator.KTrajectoryCartesian`) without any further parameters. # %% from mrpro.data import KData @@ -57,9 +51,9 @@ # %% [markdown] # Now we can explore this data object. -# Simply calling ``print(kdata)`` gives us a basic overview of the `KData` object. +# Simply printing ``kdata`` gives us a basic overview of the `KData` object. -# %% +# %% tags=["show-output"] print(kdata) @@ -74,9 +68,10 @@ # # 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 +# Let's create an FFT-operator ({py:class}`~mrpro.operator.FastFourierOp` in MRpro) and apply it to our +# {py:class}`~mrpro.data.KData` object. +# Please note thatall MRpro operator work on PyTorch tensors and not on the MRpro objects directly. Therefore, we have +# to call the operator on kdata.data. One other important property 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. # %% @@ -123,7 +118,7 @@ img_fully_sampled = img.abs().square().sum(dim=-4).sqrt().squeeze() # plot the image - +plt.imshow(img_fully_sampled) # %% [markdown] # Great! That was very easy! Let's try to reconstruct the next dataset. @@ -185,8 +180,8 @@ # 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. +# FFT-operator to, we have got the {py:class}`~mrpro.operators.CartesianSamplingOp` in MRpro. This operator performs +# sorting based on the k-space trajectory and the dimensions of the encoding k-space. # # Let's try it out! @@ -216,7 +211,7 @@ # 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. +# {py:class}`~mrpro.operators.SensitivityOp` to combine the data after image reconstruction. # %% [markdown] @@ -267,12 +262,12 @@ # 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. +# we can simply call a `{py:class}`~mrpro.algorithnms.reconstruction.DirectReconstruction`. +# Reconstruction algorithms can be instantiated 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. +# a `{py:class}`~mrpro.data.KData` object and the output is an `{py:class}`~mrpro.data.IData` object containing +# the reconstructed image data. To get its magnitude as tensor, we can call the {py:meth}`~mrpro.data.IData.rss` method. # %% from mrpro.algorithms.reconstruction import DirectReconstruction