feat: Added medical submodule

docs/source/api/index.rst

@@ -157,6 +157,17 @@ Geophysical subsurface characterization
     prestack.PrestackWaveletModelling
 
 
+Medical imaging
+~~~~~~~~~~~~~~~
+
+.. currentmodule:: pylops.medical
+
+.. autosummary::
+   :toctree: generated/
+
+    CT2D
+
+
 Solvers
 -------
 Template

docs/source/installation.rst

@@ -318,6 +318,20 @@ of GPUs should install it prior to installing PyLops as described in :ref:`Optio
 
 In alphabetic order:
 
+ASTRA
+-----
+`ASTRA <https://www.astra-toolbox.com>`_ is library used to perform computerized
+tomography. It is used in PyLops in the operator :py:class:`pylops.medical.CT2D`
+
+To use this library, install it manually either via ``conda``:
+
+.. code-block:: bash
+   >> conda install --channel astra-toolbox astra-toolbox
+or via pip:
+
+.. code-block:: bash
+   >> pip install astra-toolbox
+
 
 dtcwt
 -----

environment-dev-arm.yml

@@ -26,6 +26,7 @@ dependencies:
   - isort
   - black
   - pip:
+    - astra-toolbox
     - devito
     - dtcwt
     - scikit-fmm

environment-dev.yml

@@ -27,6 +27,7 @@ dependencies:
   - isort
   - black
   - pip:
+    - astra-toolbox
     - devito
     - dtcwt
     - scikit-fmm

pylops/__init__.py

@@ -38,6 +38,8 @@
     Linear Operators for Seismic Reservoir Characterization
 waveeqprocessing
     Linear Operators for Wave Equation oriented processing
+medical
+    Linear Operators for Medical imaging
 optimization
     Solvers
 utils
@@ -54,6 +56,7 @@
 from . import (
     avo,
     basicoperators,
+    medical,
     optimization,
     signalprocessing,
     utils,

pylops/medical/__init__.py

@@ -0,0 +1,18 @@
+"""
+Medical Operators
+=================
+
+The subpackage medical provides linear operators and applications
+aimed at solving various inverse problems in the area of Medical Imaging.
+
+A list of operators present in pylops.medical:
+
+    CT2D                            2D Computerized Tomography.
+
+"""
+
+from .ct import *
+
+__all__ = [
+    "CT2D",
+]

pylops/medical/ct.py

@@ -0,0 +1,139 @@
+__all__ = [
+    "CT2D",
+]
+
+import logging
+from typing import Optional
+
+import numpy as np
+
+from pylops import LinearOperator
+from pylops.utils import deps
+from pylops.utils.decorators import reshaped
+from pylops.utils.typing import DTypeLike, InputDimsLike, NDArray
+
+logging.basicConfig(format="%(levelname)s: %(message)s", level=logging.WARNING)
+
+
+astra_message = deps.astra_import("the astra module")
+
+if astra_message is None:
+    import astra
+
+
+class CT2D(LinearOperator):
+    r"""2D Computerized Tomography
+
+    Apply 2D computerized tomography operator to model to obtain a
+    2D sinogram.
+
+    Note that the CT2D operator is an overload of the ``astra``
+    implementation of the tomographic operator. Refer to
+    https://www.astra-toolbox.com/ for a detailed description of the
+    input parameters.
+
+    Parameters
+    ----------
+    dims : :obj:`list` or :obj:`int`
+        Number of samples for each dimension. Must be 2-dimensional and of size :math:`n_y \times n_x`
+    det_width : :obj:`float`
+        Detector width
+    det_count : :obj:`int`
+        Number of detectors
+    thetas : :obj:`numpy.ndarray`
+        Vector of angles in degrees
+    proj_geom_type : :obj:`str`, optional
+        Type of projection geometry (``parallel`` or ``fanflat``)
+    source_origin_dist : :obj:`float`, optional
+        Distance between source and origin (only for ``proj_geom_type=fanflat``)
+    origin_detector_dist : :obj:`float`, optional
+        Distance between origin and detector along the source-origin line
+        (only for "proj_geom_type=fanflat")
+    projector_type : :obj:`int`, optional
+        Type of projection geometry (``strip``, or ``line``, or ``linear``)
+    dtype : :obj:`str`, optional
+        Type of elements in input array.
+    name : :obj:`str`, optional
+        Name of operator (to be used by :func:`pylops.utils.describe.describe`)
+
+    Attributes
+    ----------
+    shape : :obj:`tuple`
+        Operator shape
+    explicit : :obj:`bool`
+        Operator contains a matrix that can be solved
+        explicitly (``True``) or not (``False``)
+
+    Notes
+    -----
+    The CT2D operator applies parallel or fan beam computerized tomography operators
+    to 2-dimensional objects and produces their corresponding sinograms.
+
+    Mathematically the forward operator can be described as [1]_:
+
+    .. math::
+        s(r,\theta; i) = \int_l i(l(r,\theta)) dl
+
+    where :math:`l(r,\theta)` is the summation line and :math:`i(x, y)`
+    is the intensity map of the model. Here, :math:`\theta` refers to the angle
+    between the y-axis (:math:`y`) and the summation line and :math:`r` is
+     the distance from the origin of the summation line.
+
+    .. [1] http://people.compute.dtu.dk/pcha/HDtomo/astra-introduction.pdf
+
+    """
+
+    def __init__(
+        self,
+        dims: InputDimsLike,
+        det_width: int,
+        det_count: float,
+        thetas: NDArray,
+        proj_geom_type: Optional[str] = "parallel",
+        source_origin_dist: float = None,
+        origin_detector_dist: float = None,
+        projector_type: Optional[str] = "strip",
+        dtype: DTypeLike = "float64",
+        name: str = "C",
+    ) -> None:
+        if astra_message is not None:
+            raise NotImplementedError(astra_message)
+
+        # create volume and projection geometries
+        self.vol_geom = astra.create_vol_geom(dims)
+        if proj_geom_type == "parallel":
+            self.proj_geom = astra.create_proj_geom(
+                proj_geom_type, det_width, det_count, thetas
+            )
+        else:
+            self.proj_geom = astra.create_proj_geom(
+                proj_geom_type,
+                det_width,
+                det_count,
+                thetas,
+                source_origin_dist,
+                origin_detector_dist,
+            )
+
+        # create projector
+        self.proj_id = astra.create_projector(
+            projector_type, self.proj_geom, self.vol_geom
+        )
+        super().__init__(
+            dtype=np.dtype(dtype), dims=dims, dimsd=(len(thetas), det_count), name=name
+        )
+
+    @reshaped
+    def _matvec(self, x):
+        y_id, y = astra.create_sino(x, self.proj_id)
+        astra.data2d.delete(y_id)
+        return y
+
+    @reshaped
+    def _rmatvec(self, x):
+        y_id, y = astra.create_backprojection(x, self.proj_id)
+        astra.data2d.delete(y_id)
+        return y
+
+    def __del__(self):
+        astra.projector.delete(self.proj_id)

pylops/utils/deps.py

@@ -1,4 +1,5 @@
 __all__ = [
+    "astra_enabled",
     "cupy_enabled",
     "jax_enabled",
     "devito_enabled",
@@ -81,6 +82,23 @@ def jax_import(message: Optional[str] = None) -> str:
     return jax_message
 
 
+def astra_import(message):
+    if astra_enabled:
+        try:
+            import_module("astra")  # noqa: F401
+
+            astra_message = None
+        except Exception as e:
+            astra_message = f"Failed to import astra (error:{e})."
+    else:
+        astra_message = (
+            f"ASTRA not available. "
+            f"In order to be able to use "
+            f'{message} run "pip install astra-toolbox".'
+        )
+    return astra_message
+
+
 def devito_import(message: Optional[str] = None) -> str:
     if devito_enabled:
         try:
@@ -101,7 +119,7 @@ def devito_import(message: Optional[str] = None) -> str:
 def dtcwt_import(message: Optional[str] = None) -> str:
     if dtcwt_enabled:
         try:
-            import dtcwt  # noqa: F401
+            import_module("dtcwt")  # noqa: F401
 
             dtcwt_message = None
         except Exception as e:
@@ -254,6 +272,7 @@ def pytensor_import(message: Optional[str] = None) -> str:
 jax_enabled: bool = (
     True if (jax_import() is None and int(os.getenv("JAX_PYLOPS", 1)) == 1) else False
 )
+astra_enabled = util.find_spec("astra") is not None
 devito_enabled = util.find_spec("devito") is not None
 dtcwt_enabled = util.find_spec("dtcwt") is not None
 numba_enabled = util.find_spec("numba") is not None

pyproject.toml

@@ -44,6 +44,7 @@ advanced = [
     "scikit-fmm",
     "spgl1",
     "dtcwt",
+    "astra-toolbox",
 ]
 
 [tool.setuptools.packages.find]

pytests/test_ct.py

@@ -0,0 +1,52 @@
+import numpy as np
+import pytest
+
+from pylops.medical import CT2D
+from pylops.utils import dottest
+
+par1 = {
+    "ny": 51,
+    "nx": 30,
+    "ntheta": 20,
+    "proj_geom_type": "parallel",
+    "projector_type": "strip",
+    "dtype": "float64",
+}  # parallel, strip
+
+par2 = {
+    "ny": 51,
+    "nx": 30,
+    "ntheta": 20,
+    "proj_geom_type": "parallel",
+    "projector_type": "line",
+    "dtype": "float64",
+}  # parallel, line
+
+par3 = {
+    "ny": 51,
+    "nx": 30,
+    "ntheta": 20,
+    "proj_geom_type": "fanflat",
+    "projector_type": "strip",
+    "dtype": "float64",
+}  # fanflat, strip
+
+
+@pytest.mark.parametrize("par", [(par1), (par2)])
+def test_CT2D(par):
+    """Dot-test for CT2D operator"""
+    theta = np.linspace(0.0, np.pi, par["ntheta"], endpoint=False)
+
+    Cop = CT2D(
+        (par["ny"], par["nx"]),
+        1.0,
+        par["ny"],
+        theta,
+        proj_geom_type=par["proj_geom_type"],
+        projector_type=par["projector_type"],
+    )
+    assert dottest(
+        Cop,
+        par["ny"] * par["ntheta"],
+        par["ny"] * par["nx"],
+    )

requirements-dev.txt

@@ -9,6 +9,7 @@ scikit-fmm
 sympy
 devito
 dtcwt
+astra-toolbox
 matplotlib
 ipython
 pytest

requirements-doc.txt

@@ -15,6 +15,7 @@ scikit-fmm
 sympy
 devito==4.8.6
 dtcwt
+astra-toolbox
 matplotlib
 ipython
 pytest