Various improvements related to the new X-ray transform implemementat…

…ion. (#461) * Rename TomographicProjector to XRayTransform * Remove markup from exception messages and other minor edits * Rename modules * Rename files and rename TomographicProjector to XRayTransform * Rename AbelProjector to AbelTransform * Docs edits * Update change summary * Overlooked rename changes * Replace Radon transform label * Update submodule * Renaming of some CT projectors (#453) * Add 2d projector and code to time it * Clean up * Add back projection * Add test * Add timing results to example * Start to add new example * Update data * Address mypy * Address isort * Try to fix tables in notebook * Update submodule * Rename XRayProject to XRayTransform * Minor edits --------- Co-authored-by: Michael McCann <[email protected]> Co-authored-by: Michael-T-McCann <[email protected]> Co-authored-by: Michael McCann <[email protected]> * Restructure X-ray transform modules * Update submodule * Adjust angles to be equivalent between scico and astra projections * Clean up output * New example script * Minor edits * Shorten class name * Clarify Parallel2dProjector angles * Docs edits * Remove problematic jit * New example script * Docs improvement * Rename parameter * Docs improvement * Typo fix * Add noise * Add example script * Change noise level * Update notebooks * Docs fix * Update submodule * Add warning to api docs * Add overlooked change summary entry * Remove unintentionally added file * Remove unintentionally added files * Remove unintentionally added files * Address review comment * Remove unintentionally added file * Fix docs typo * Update submodule * Update submodule --------- Co-authored-by: Michael McCann <[email protected]> Co-authored-by: Michael-T-McCann <[email protected]> Co-authored-by: Michael McCann <[email protected]>
lanl · Nov 2, 2023 · 27e2aec · 27e2aec
1 parent bcb1eab
commit 27e2aec
Show file tree

Hide file tree

Showing 38 changed files with 761 additions and 252 deletions.
diff --git a/CHANGES.rst b/CHANGES.rst
@@ -6,6 +6,11 @@ SCICO Release Notes
 Version 0.0.5   (unreleased)
 ----------------------------
 
+• New integrated Radon/X-ray transform ``linop.XRayTransform``.
+• Rename modules ``radon_astra`` and ``radon_svmbir`` to ``xray.astra`` and
+  ``xray.svmbir`` respectively, and rename ``TomographicProjector`` classes
+  to ``XRayTransform``.
+• Rename ``AbelProjector`` to ``AbelTransform``.
 • Rename ``solver.ATADSolver`` to ``solver.MatrixATADSolver``.
 • Support ``jaxlib`` and ``jax`` versions 0.4.3 to 0.4.19.
 

diff --git a/data b/data
diff --git a/docs/source/examples.rst b/docs/source/examples.rst
@@ -36,7 +36,8 @@ Computed Tomography
    examples/ct_astra_odp_train_foam2
    examples/ct_astra_unet_train_foam2
    examples/ct_projector_comparison
-
+   examples/ct_multi_cs_tv_admm
+   examples/ct_multi_tv_admm
 
 Deconvolution
 ^^^^^^^^^^^^^

diff --git a/docs/source/inverse.rst b/docs/source/inverse.rst
@@ -47,15 +47,15 @@ SCICO provides the :class:`.Operator` and :class:`.LinearOperator`
 classes, which may be subclassed by users, in order to implement the
 forward operator, :math:`A`. It also has several built-in operators,
 most of which are linear, e.g., finite convolutions, discrete Fourier
-transforms, optical propagators, Abel transforms, and Radon
-transforms. For example,
+transforms, optical propagators, Abel transforms, and X-ray transforms
+(the same as Radon transforms in 2D). For example,
 
 .. code:: python
 
        input_shape = (512, 512)
        angles = np.linspace(0, 2 * np.pi, 180, endpoint=False)
        channels = 512
-       A = scico.linop.radon_svmbir.ParallelBeamProjector(input_shape, angles, channels)
+       A = scico.linop.xray.svmbir.XRayTransform(input_shape, angles, channels)
 
 defines a tomographic projection operator.
 

diff --git a/docs/source/notes.rst b/docs/source/notes.rst
@@ -111,13 +111,25 @@ via interfaces to the `bm3d <https://pypi.org/project/bm3d/>`__ and
 when the full benefits of JAX-based code are required.
 
 
-Tomographic Projectors
-----------------------
-
-The :class:`.radon_svmbir.TomographicProjector` class is implemented
+Tomographic Projectors/Radon Transforms
+---------------------------------------
+
+Note that the tomographic projections that are frequently referred
+to as Radon transforms are referred to as X-ray transforms in SCICO.
+While the Radon transform is far more well-known than the X-ray
+transform, which is the same as the Radon transform for projections
+in two dimensions, these two transform differ in higher numbers of
+dimensions, and it is the X-ray transform that is the appropriate
+mathematical model for beam attenuation based imaging in three or
+more dimensions.
+
+SCICO includes three different implementations of X-ray transforms.
+Of these, :class:`.linop.XRayTransform` is an integral component of
+SCICO, while the other two depend on external packages.
+The :class:`.xray.svmbir.XRayTransform` class is implemented
 via an interface to the `svmbir
 <https://svmbir.readthedocs.io/en/latest/>`__ package. The
-:class:`.radon_astra.TomographicProjector` class is implemented via an
+:class:`.xray.astra.XRayTransform` class is implemented via an
 interface to the `ASTRA toolbox
 <https://www.astra-toolbox.com/>`__. This toolbox does provide some
 GPU acceleration support, but efficiency is expected to be lower than

diff --git a/examples/scripts/README.rst b/examples/scripts/README.rst
@@ -36,8 +36,11 @@ Computed Tomography
    `ct_astra_unet_train_foam2.py <ct_astra_unet_train_foam2.py>`_
       CT Training and Reconstructions with UNet
    `ct_projector_comparison.py <ct_projector_comparison.py>`_
-      X-ray Projector Comparison
-
+      X-ray Transform Comparison
+   `ct_multi_cs_tv_admm.py <ct_multi_cs_tv_admm.py>`_
+      TV-Regularized Sparse-View CT Reconstruction (Multiple Projectors, Common Sinogram)
+   `ct_multi_tv_admm.py <ct_multi_tv_admm.py>`_
+      TV-Regularized Sparse-View CT Reconstruction (Multiple Projectors)
 
 Deconvolution
 ^^^^^^^^^^^^^

diff --git a/examples/scripts/ct_abel_tv_admm.py b/examples/scripts/ct_abel_tv_admm.py
@@ -26,7 +26,7 @@
 import scico.numpy as snp
 from scico import functional, linop, loss, metric, plot
 from scico.examples import create_circular_phantom
-from scico.linop.abel import AbelProjector
+from scico.linop.abel import AbelTransform
 from scico.optimize.admm import ADMM, LinearSubproblemSolver
 from scico.util import device_info
 
@@ -40,7 +40,7 @@
 """
 Set up the forward operator and create a test measurement.
 """
-A = AbelProjector(x_gt.shape)
+A = AbelTransform(x_gt.shape)
 y = A @ x_gt
 np.random.seed(12345)
 y = y + np.random.normal(size=y.shape).astype(np.float32)

diff --git a/examples/scripts/ct_abel_tv_admm_tune.py b/examples/scripts/ct_abel_tv_admm_tune.py
@@ -38,7 +38,7 @@
 import scico.numpy as snp
 from scico import functional, linop, loss, metric, plot
 from scico.examples import create_circular_phantom
-from scico.linop.abel import AbelProjector
+from scico.linop.abel import AbelTransform
 from scico.optimize.admm import ADMM, LinearSubproblemSolver
 from scico.ray import tune
 
@@ -52,7 +52,7 @@
 """
 Set up the forward operator and create a test measurement.
 """
-A = AbelProjector(x_gt.shape)
+A = AbelTransform(x_gt.shape)
 y = A @ x_gt
 np.random.seed(12345)
 y = y + np.random.normal(size=y.shape).astype(np.float32)
@@ -84,7 +84,7 @@ def setup(self, config, x_gt, x0, y):
         # Get arrays passed by tune call.
         self.x_gt, self.x0, self.y = snp.array(x_gt), snp.array(x0), snp.array(y)
         # Set up problem to be solved.
-        self.A = AbelProjector(self.x_gt.shape)
+        self.A = AbelTransform(self.x_gt.shape)
         self.f = loss.SquaredL2Loss(y=self.y, A=self.A)
         self.C = linop.FiniteDifference(input_shape=self.x_gt.shape)
         self.reset_config(config)

diff --git a/examples/scripts/ct_astra_3d_tv_admm.py b/examples/scripts/ct_astra_3d_tv_admm.py
@@ -15,9 +15,9 @@
   $$\mathrm{argmin}_{\mathbf{x}} \; (1/2) \| \mathbf{y} - A \mathbf{x}
   \|_2^2 + \lambda \| C \mathbf{x} \|_{2,1} \;,$$
 
-where $A$ is the Radon transform, $\mathbf{y}$ is the sinogram, $C$ is
-a 3D finite difference operator, and $\mathbf{x}$ is the desired
-image.
+where $A$ is the X-ray transform (the CT forward projection operator),
+$\mathbf{y}$ is the sinogram, $C$ is a 3D finite difference operator,
+and $\mathbf{x}$ is the desired image.
 """
 
 
@@ -28,7 +28,7 @@
 import scico.numpy as snp
 from scico import functional, linop, loss, metric, plot
 from scico.examples import create_tangle_phantom
-from scico.linop.radon_astra import TomographicProjector
+from scico.linop.xray.astra import XRayTransform
 from scico.optimize.admm import ADMM, LinearSubproblemSolver
 from scico.util import device_info
 
@@ -43,9 +43,7 @@
 
 n_projection = 10  # number of projections
 angles = np.linspace(0, np.pi, n_projection)  # evenly spaced projection angles
-A = TomographicProjector(
-    tangle.shape, [1.0, 1.0], [Nz, max(Nx, Ny)], angles
-)  # Radon transform operator
+A = XRayTransform(tangle.shape, [1.0, 1.0], [Nz, max(Nx, Ny)], angles)  # CT projection operator
 y = A @ tangle  # sinogram
 
 

diff --git a/examples/scripts/ct_astra_modl_train_foam2.py b/examples/scripts/ct_astra_modl_train_foam2.py
@@ -54,7 +54,7 @@
 from scico import metric, plot
 from scico.flax.examples import load_ct_data
 from scico.flax.train.traversals import clip_positive, construct_traversal
-from scico.linop.radon_astra import TomographicProjector
+from scico.linop.xray.astra import XRayTransform
 
 """
 Prepare parallel processing. Set an arbitrary processor count (only
@@ -81,12 +81,12 @@
 Build CT projection operator.
 """
 angles = np.linspace(0, np.pi, n_projection)  # evenly spaced projection angles
-A = TomographicProjector(
+A = XRayTransform(
     input_shape=(N, N),
     detector_spacing=1,
     det_count=N,
     angles=angles,
-)  # Radon transform operator
+)  # CT projection operator
 A = (1.0 / N) * A  # normalized
 
 

diff --git a/examples/scripts/ct_astra_noreg_pcg.py b/examples/scripts/ct_astra_noreg_pcg.py
@@ -15,8 +15,8 @@
   $$\mathrm{argmin}_{\mathbf{x}} \; (1/2) \| \mathbf{y} - A \mathbf{x}
   \|_2^2 \;,$$
 
-where $A$ is the Radon transform, $\mathbf{y}$ is the sinogram, and
-$\mathbf{x}$ is the reconstructed image.
+where $A$ is the X-ray transform (the CT forward projection operator),
+$\mathbf{y}$ is the sinogram, and $\mathbf{x}$ is the reconstructed image.
 """
 
 from time import time
@@ -29,7 +29,7 @@
 
 from scico import loss, plot
 from scico.linop import CircularConvolve
-from scico.linop.radon_astra import TomographicProjector
+from scico.linop.xray.astra import XRayTransform
 from scico.solver import cg
 
 """
@@ -45,7 +45,7 @@
 """
 n_projection = N  # matches the phantom size so this is not few-view CT
 angles = np.linspace(0, np.pi, n_projection)  # evenly spaced projection angles
-A = 1 / N * TomographicProjector(x_gt.shape, 1, N, angles)  # Radon transform operator
+A = 1 / N * XRayTransform(x_gt.shape, 1, N, angles)  # CT projection operator
 y = A @ x_gt  # sinogram
 
 

diff --git a/examples/scripts/ct_astra_odp_train_foam2.py b/examples/scripts/ct_astra_odp_train_foam2.py
@@ -58,7 +58,7 @@
 from scico import metric, plot
 from scico.flax.examples import load_ct_data
 from scico.flax.train.traversals import clip_positive, construct_traversal
-from scico.linop.radon_astra import TomographicProjector
+from scico.linop.xray.astra import XRayTransform
 
 """
 Prepare parallel processing. Set an arbitrary processor count (only
@@ -85,12 +85,12 @@
 Build CT projection operator.
 """
 angles = np.linspace(0, np.pi, n_projection)  # evenly spaced projection angles
-A = TomographicProjector(
+A = XRayTransform(
     input_shape=(N, N),
     detector_spacing=1,
     det_count=N,
     angles=angles,
-)  # Radon transform operator
+)  # CT projection operator
 A = (1.0 / N) * A  # normalized
 
 

diff --git a/examples/scripts/ct_astra_tv_admm.py b/examples/scripts/ct_astra_tv_admm.py
@@ -14,9 +14,9 @@
   $$\mathrm{argmin}_{\mathbf{x}} \; (1/2) \| \mathbf{y} - A \mathbf{x}
   \|_2^2 + \lambda \| C \mathbf{x} \|_{2,1} \;,$$
 
-where $A$ is the Radon transform, $\mathbf{y}$ is the sinogram, $C$ is
-a 2D finite difference operator, and $\mathbf{x}$ is the desired
-image.
+where $A$ is the X-ray transform (the CT forward projection operator),
+$\mathbf{y}$ is the sinogram, $C$ is a 2D finite difference operator, and
+$\mathbf{x}$ is the desired image.
 """
 
 import numpy as np
@@ -26,7 +26,7 @@
 
 import scico.numpy as snp
 from scico import functional, linop, loss, metric, plot
-from scico.linop.radon_astra import TomographicProjector
+from scico.linop.xray.astra import XRayTransform
 from scico.optimize.admm import ADMM, LinearSubproblemSolver
 from scico.util import device_info
 
@@ -44,7 +44,7 @@
 """
 n_projection = 45  # number of projections
 angles = np.linspace(0, np.pi, n_projection)  # evenly spaced projection angles
-A = TomographicProjector(x_gt.shape, 1, N, angles)  # Radon transform operator
+A = XRayTransform(x_gt.shape, 1, N, angles)  # CT projection operator
 y = A @ x_gt  # sinogram
 
 

diff --git a/examples/scripts/ct_astra_weighted_tv_admm.py b/examples/scripts/ct_astra_weighted_tv_admm.py
@@ -14,11 +14,11 @@
   $$\mathrm{argmin}_{\mathbf{x}} \; (1/2) \| \mathbf{y} - A \mathbf{x}
   \|_W^2 + \lambda \| C \mathbf{x} \|_{2,1} \;,$$
 
-where $A$ is the Radon transform, $\mathbf{y}$ is the sinogram, the norm
-weighting $W$ is chosen so that the weighted norm is an approximation to
-the Poisson negative log likelihood :cite:`sauer-1993-local`, $C$ is
-a 2D finite difference operator, and $\mathbf{x}$ is the desired
-image.
+where $A$ is the X-ray transform (the CT forward projection),
+$\mathbf{y}$ is the sinogram, the norm weighting $W$ is chosen so that
+the weighted norm is an approximation to the Poisson negative log
+likelihood :cite:`sauer-1993-local`, $C$ is a 2D finite difference
+operator, and $\mathbf{x}$ is the desired image.
 """
 
 import numpy as np
@@ -27,7 +27,7 @@
 
 import scico.numpy as snp
 from scico import functional, linop, loss, metric, plot
-from scico.linop.radon_astra import TomographicProjector
+from scico.linop.xray.astra import XRayTransform
 from scico.optimize.admm import ADMM, LinearSubproblemSolver
 from scico.util import device_info
 
@@ -51,7 +51,7 @@
 𝛼 = 1e-2  # attenuation coefficient
 
 angles = np.linspace(0, 2 * np.pi, n_projection)  # evenly spaced projection angles
-A = TomographicProjector(x_gt.shape, 1.0, N, angles)  # Radon transform operator
+A = XRayTransform(x_gt.shape, 1.0, N, angles)  # CT projection operator
 y_c = A @ x_gt  # sinogram
 
 

diff --git a/examples/scripts/ct_fan_svmbir_ppp_bm3d_admm_prox.py b/examples/scripts/ct_fan_svmbir_ppp_bm3d_admm_prox.py
@@ -35,7 +35,7 @@
 from scico import metric, plot
 from scico.functional import BM3D
 from scico.linop import Diagonal, Identity
-from scico.linop.radon_svmbir import SVMBIRExtendedLoss, TomographicProjector
+from scico.linop.xray.svmbir import SVMBIRExtendedLoss, XRayTransform
 from scico.optimize.admm import ADMM, LinearSubproblemSolver
 from scico.util import device_info
 
@@ -65,15 +65,15 @@
 
 dist_source_detector = 1500.0
 magnification = 1.2
-A_fan = TomographicProjector(
+A_fan = XRayTransform(
     x_gt.shape,
     angles,
     num_channels,
     geometry="fan-curved",
     dist_source_detector=dist_source_detector,
     magnification=magnification,
 )
-A_parallel = TomographicProjector(
+A_parallel = XRayTransform(
     x_gt.shape,
     angles,
     num_channels,