Merge pull request #733 from mrava87/feat-mri

Feat: mri module

Author: mrava87

docs/source/api/index.rst

@@ -171,6 +171,7 @@ Medical imaging
    :toctree: generated/
 
     CT2D
+    MRI2D
 
 
 Solvers

docs/source/faq.rst

@@ -14,7 +14,7 @@ working with linear operators is indeed that you don't really need to access the
 of an operator.
 
 
-**2. Can I have an older version of** ``cupy`` **installed in my system (** ``cupy-cudaXX<10.6.0``)?**
+**2. Can I have an older version of** ``cupy`` **installed in my system (** ``cupy-cudaXX<10.6.0``**)?**
 
 Yes. Nevertheless you need to tell PyLops that you don't want to use its ``cupy``
 backend by setting the environment variable ``CUPY_PYLOPS=0``.

docs/source/gpu.rst

@@ -602,6 +602,10 @@ Medical:
      - |:white_check_mark:|
      - |:white_check_mark:|
      - |:white_check_mark:|
+   * - :class:`pylops.medical.mri.MRI2D`
+     - |:white_check_mark:|
+     - |:red_circle:|
+     - |:white_check_mark:|
 
 .. warning::
 

examples/plot_avo.py

@@ -1,6 +1,6 @@
 r"""
 AVO modelling
-===================
+=============
 This example shows how to create pre-stack angle gathers using
 the :py:class:`pylops.avo.avo.AVOLinearModelling` operator.
 """

examples/plot_mri.py

@@ -0,0 +1,104 @@
+r"""
+MRI modelling
+=============
+This example shows how to use the  :py:class:`pylops.medical.mri.MRI2D` operator
+to create K-space undersampled MRI data.
+"""
+import matplotlib.pyplot as plt
+import numpy as np
+
+import pylops
+
+plt.close("all")
+np.random.seed(0)
+
+###############################################################################
+# Let"s start by loading the Shepp-Logan phantom model.
+x = np.load("../testdata/optimization/shepp_logan_phantom.npy")
+x = x / x.max()
+nx, ny = x.shape
+
+###############################################################################
+# Next, we create a mask to simulate undersampling in K-space and apply it to
+# the phantom model.
+
+# Passing mask as array
+mask = np.zeros((nx, ny))
+mask[:, np.random.randint(0, ny, 2 * ny // 3)] = 1
+mask[:, ny // 2 - 20 : ny // 2 + 10] = 1
+
+Mop = pylops.medical.MRI2D(dims=(nx, ny), mask=mask)
+
+d = Mop @ x
+x_adj = Mop.H @ d
+
+fig, axs = plt.subplots(1, 3, figsize=(12, 5))
+axs[0].imshow(x, cmap="gray", vmin=0, vmax=1)
+axs[0].set_title("Original Image")
+axs[1].imshow(np.abs(d), cmap="jet", vmin=0, vmax=1)
+axs[1].set_title("K-space Data")
+axs[2].imshow(x_adj.real, cmap="gray", vmin=0, vmax=1)
+axs[2].set_title("Adjoint Reconstruction")
+fig.tight_layout()
+
+###############################################################################
+# Alternatively, we can create the same mask by specifying a sampling pattern
+# using the ``mask`` keyword argument. Here, we create a ``vertical-reg`` mask
+# that samples K-space lines in the vertical direction with a regular pattern.
+
+# Vertical uniform with center
+Mop = pylops.medical.MRI2D(
+    dims=(nx, ny), mask="vertical-reg", nlines=ny // 2, perc_center=0.0
+)
+
+d = Mop @ x
+x_adj = (Mop.H @ d).reshape(nx, ny)
+
+fig, axs = plt.subplots(1, 3, figsize=(12, 5))
+axs[0].imshow(x, cmap="gray", vmin=0, vmax=1)
+axs[0].set_title("Original Image")
+axs[1].imshow(np.abs(Mop.ROp.H @ d).reshape(nx, ny), cmap="jet", vmin=0, vmax=1)
+axs[1].set_title("K-space Data")
+axs[2].imshow(x_adj.real, cmap="gray", vmin=0, vmax=1)
+axs[2].set_title("Adjoint Reconstruction")
+fig.tight_layout()
+
+###############################################################################
+# Similarly, we can create a ``vertical-uni`` mask that randomly samples
+# K-space lines in the vertical direction.
+
+# Vertical uniform with center
+Mop = pylops.medical.MRI2D(
+    dims=(nx, ny), mask="vertical-uni", nlines=40, perc_center=0.1
+)
+
+d = Mop @ x
+x_adj = (Mop.H @ d).reshape(nx, ny)
+
+fig, axs = plt.subplots(1, 3, figsize=(12, 5))
+axs[0].imshow(x, cmap="gray", vmin=0, vmax=1)
+axs[0].set_title("Original Image")
+axs[1].imshow(np.abs(Mop.ROp.H @ d).reshape(nx, ny), cmap="jet", vmin=0, vmax=1)
+axs[1].set_title("K-space Data")
+axs[2].imshow(x_adj.real, cmap="gray", vmin=0, vmax=1)
+axs[2].set_title("Adjoint Reconstruction")
+fig.tight_layout()
+
+###############################################################################
+# Finally, we can create a sampling pattern with radial lines using the
+# ``radial-uni`` (or ``radial-reg``) option.
+
+# Radial uniform
+Mop = pylops.medical.MRI2D(dims=(nx, ny), mask="radial-uni", nlines=40)
+
+d = Mop @ x
+x_adj = (Mop.H @ d).reshape(nx, ny)
+
+fig, axs = plt.subplots(1, 3, figsize=(12, 5))
+axs[0].imshow(x, cmap="gray", vmin=0, vmax=1)
+axs[0].set_title("Original Image")
+axs[1].imshow(np.abs(Mop.ROp.H @ d).reshape(nx, ny), cmap="jet", vmin=0, vmax=1)
+axs[1].set_title("K-space Data")
+axs[2].imshow(x_adj.real, cmap="gray", vmin=0, vmax=1)
+axs[2].set_title("Adjoint Reconstruction")
+fig.tight_layout()

pylops/medical/__init__.py

@@ -8,11 +8,12 @@
 A list of operators present in pylops.medical:
 
     CT2D                            2D Computerized Tomography.
+    MRI2D                           2D Magnetic Resonance Imaging.
 
 """
 
 from .ct import *
+from .mri import *
 
-__all__ = [
-    "CT2D",
-]
+
+__all__ = ["CT2D", "MRI2D"]

pylops/medical/ct.py

@@ -63,6 +63,14 @@ class CT2D(LinearOperator):
 
     Attributes
     ----------
+    dims : :obj:`tuple`
+        Shape of the array after the adjoint, but before flattening.
+
+        For example, ``x_reshaped = (Op.H * y.ravel()).reshape(Op.dims)``.
+    dimsd : :obj:`tuple`
+        Shape of the array after the forward, but before flattening.
+
+        For example, ``y_reshaped = (Op * x.ravel()).reshape(Op.dimsd)``.
     shape : :obj:`tuple`
         Operator shape
     explicit : :obj:`bool`