Merge branch 'main' into brendt/trace

CHANGES.rst

@@ -3,10 +3,20 @@ SCICO Release Notes
 ===================
 
 
-Version 0.0.6   (unreleased)
+Version 0.0.7   (unreleased)
+----------------------------
+
+• No changes yet.
+
+
+
+Version 0.0.6   (2024-10-25)
 ----------------------------
 
 • Significant changes to ``linop.xray.astra`` API.
+• Rename integrated 2D X-ray transform class to
+  ``linop.xray.XRayTransform2D`` and add filtered back projection method
+  ``fbp``.
 • New integrated 3D X-ray transform via ``linop.xray.XRayTransform3D``.
 • New functional ``functional.IsotropicTVNorm`` and faster implementation
   of ``functional.AnisotropicTVNorm``.
@@ -18,8 +28,8 @@ Version 0.0.6   (unreleased)
 • Rename ``scico.flax.save_weights`` and ``scico.flax.load_weights`` to
   ``scico.flax.save_variables`` and ``scico.flax.load_variables``
   respectively.
-• Support ``jaxlib`` and ``jax`` versions 0.4.3 to 0.4.31.
-• Support ``flax`` versions 0.8.0 to 0.8.3.
+• Support ``jaxlib`` and ``jax`` versions 0.4.13 to 0.4.35.
+• Support ``flax`` versions 0.8.0 to 0.10.0.
 
 
 

MANIFEST.in

@@ -8,10 +8,10 @@ include pyproject.toml
 include pytest.ini
 include requirements.txt
 include dev_requirements.txt
-include examples/scriptcheck.sh
 include docs/docs_requirements.txt
 
 recursive-include scico *.py
-recursive-include scico/data *.png *.npz
-recursive-include docs Makefile *.py *.ipynb *.rst *.bib *.css *.svg *.ico
-recursive-include examples *_requirements.txt *.txt *.rst *.py
+recursive-include scico/data *.png *.mpk *.rst
+recursive-include docs Makefile *.py *.ipynb *.rst *.bib *.css *.svg *.png *.ico
+recursive-include examples *_requirements.txt *.txt *.rst *.py *.sh
+recursive-include misc *.py *.sh *.rst

README.md

@@ -32,11 +32,12 @@ this software for published work, please cite the corresponding [JOSS
 Paper](https://doi.org/10.21105/joss.04722) (see bibtex entry
 `balke-2022-scico` in `docs/source/references.bib`).
 
+
 # Installation
 
-See the [online
-documentation](https://scico.rtfd.io/en/latest/install.html) for
-installation instructions.
+The online documentation includes detailed
+[installation instructions](https://scico.rtfd.io/en/latest/install.html).
+
 
 # Usage Examples
 
@@ -47,8 +48,11 @@ Jupyter Notebooks are provided in the
 to `examples/notebooks`. They are also viewable on
 [GitHub](https://github.com/lanl/scico-data/tree/main/notebooks) or
 [nbviewer](https://nbviewer.jupyter.org/github/lanl/scico-data/tree/main/notebooks/index.ipynb),
-or can be run online by
-[binder](https://mybinder.org/v2/gh/lanl/scico-data/binder?labpath=notebooks%2Findex.ipynb).
+and can be run online on
+[binder](https://mybinder.org/v2/gh/lanl/scico-data/binder?labpath=notebooks%2Findex.ipynb)
+or
+[google colab](https://colab.research.google.com/github/lanl/scico-data/blob/colab/notebooks/index.ipynb).
+
 
 # License
 

docs/source/examples.rst

@@ -34,12 +34,11 @@ Computed Tomography
    examples/ct_svmbir_ppp_bm3d_admm_cg
    examples/ct_svmbir_ppp_bm3d_admm_prox
    examples/ct_fan_svmbir_ppp_bm3d_admm_prox
-   examples/ct_astra_modl_train_foam2
-   examples/ct_astra_odp_train_foam2
-   examples/ct_astra_unet_train_foam2
+   examples/ct_modl_train_foam2
+   examples/ct_odp_train_foam2
+   examples/ct_unet_train_foam2
    examples/ct_projector_comparison_2d
    examples/ct_projector_comparison_3d
-   examples/ct_multi_cs_tv_admm
    examples/ct_multi_tv_admm
 
 Deconvolution
@@ -96,7 +95,7 @@ Miscellaneous
    examples/denoise_dncnn_universal
    examples/diffusercam_tv_admm
    examples/video_rpca_admm
-   examples/ct_astra_datagen_foam2
+   examples/ct_datagen_foam2
    examples/deconv_datagen_bsds
    examples/deconv_datagen_foam1
    examples/denoise_datagen_bsds
@@ -181,10 +180,10 @@ Machine Learning
 .. toctree::
    :maxdepth: 1
 
-   examples/ct_astra_datagen_foam2
-   examples/ct_astra_modl_train_foam2
-   examples/ct_astra_odp_train_foam2
-   examples/ct_astra_unet_train_foam2
+   examples/ct_datagen_foam2
+   examples/ct_modl_train_foam2
+   examples/ct_odp_train_foam2
+   examples/ct_unet_train_foam2
    examples/deconv_datagen_bsds
    examples/deconv_datagen_foam1
    examples/deconv_modl_train_foam1

docs/source/references.bib

@@ -396,6 +396,13 @@ @Article {jin-2017-unet
   doi =		 {10.1109/TIP.2017.2713099}
 }
 
+@Book {kak-1988-principles,
+  author = 	 {Avinash C. Kak and Malcolm Slaney},
+  title = 	 {Principles of Computerized Tomographic Imaging},
+  publisher = 	 {IEEE Press},
+  year = 	 1988
+}
+
 @TechReport {kamilov-2016-minimizing,
   author =	 {Ulugbek S. Kamilov},
   title =	 {Minimizing Isotropic Total Variation without
@@ -771,6 +778,7 @@ @Article {zhang-2017-dncnn
   pages =	 {3142--3155}
 }
 
+
 @Article {zhang-2021-plug,
   author =	 {Zhang, Kai and Li, Yawei and Zuo, Wangmeng and
                   Zhang, Lei and Van Gool, Luc and Timofte, Radu},

examples/jnb.py

@@ -62,10 +62,10 @@ def py_file_to_string(src):
 
         # Process remainder of source file
         for line in srcfile:
-            if re.match("^input\(", line):  # end processing when input statement encountered
+            if re.match(r"^input\(", line):  # end processing when input statement encountered
                 break
             line = re.sub('^r"""', '"""', line)  # remove r from r"""
-            line = re.sub(":cite:\`([^`]+)\`", r'<cite data-cite="\1"/>', line)  # fix cite format
+            line = re.sub(r":cite:\`([^`]+)\`", r'<cite data-cite="\1"/>', line)  # fix cite format
             lines.append(line)
 
         # Backtrack through list of lines to remove trailing newlines

examples/scripts/README.rst

@@ -33,11 +33,11 @@ Computed Tomography
       PPP (with BM3D) CT Reconstruction (ADMM with Fast SVMBIR Prox)
    `ct_fan_svmbir_ppp_bm3d_admm_prox.py <ct_fan_svmbir_ppp_bm3d_admm_prox.py>`_
       PPP (with BM3D) Fan-Beam CT Reconstruction
-   `ct_astra_modl_train_foam2.py <ct_astra_modl_train_foam2.py>`_
-      CT Training and Reconstructions with MoDL
-   `ct_astra_odp_train_foam2.py <ct_astra_odp_train_foam2.py>`_
-      CT Training and Reconstructions with ODP
-   `ct_astra_unet_train_foam2.py <ct_astra_unet_train_foam2.py>`_
+   `ct_modl_train_foam2.py <ct_modl_train_foam2.py>`_
+      CT Training and Reconstruction with MoDL
+   `ct_odp_train_foam2.py <ct_odp_train_foam2.py>`_
+      CT Training and Reconstruction with ODP
+   `ct_unet_train_foam2.py <ct_unet_train_foam2.py>`_
       CT Training and Reconstructions with UNet
    `ct_projector_comparison_2d.py <ct_projector_comparison_2d.py>`_
       2D X-ray Transform Comparison
@@ -123,7 +123,7 @@ Miscellaneous
       TV-Regularized 3D DiffuserCam Reconstruction
    `video_rpca_admm.py <video_rpca_admm.py>`_
       Video Decomposition via Robust PCA
-   `ct_astra_datagen_foam2.py <ct_astra_datagen_foam2.py>`_
+   `ct_datagen_foam2.py <ct_datagen_foam2.py>`_
       CT Data Generation for NN Training
    `deconv_datagen_bsds.py <deconv_datagen_bsds.py>`_
       Blurred Data Generation (Natural Images) for NN Training
@@ -239,13 +239,13 @@ Sparsity
 Machine Learning
 ^^^^^^^^^^^^^^^^
 
-   `ct_astra_datagen_foam2.py <ct_astra_datagen_foam2.py>`_
+   `ct_datagen_foam2.py <ct_datagen_foam2.py>`_
       CT Data Generation for NN Training
-   `ct_astra_modl_train_foam2.py <ct_astra_modl_train_foam2.py>`_
-      CT Training and Reconstructions with MoDL
-   `ct_astra_odp_train_foam2.py <ct_astra_odp_train_foam2.py>`_
-      CT Training and Reconstructions with ODP
-   `ct_astra_unet_train_foam2.py <ct_astra_unet_train_foam2.py>`_
+   `ct_modl_train_foam2.py <ct_modl_train_foam2.py>`_
+      CT Training and Reconstruction with MoDL
+   `ct_odp_train_foam2.py <ct_odp_train_foam2.py>`_
+      CT Training and Reconstruction with ODP
+   `ct_unet_train_foam2.py <ct_unet_train_foam2.py>`_
       CT Training and Reconstructions with UNet
    `deconv_datagen_bsds.py <deconv_datagen_bsds.py>`_
       Blurred Data Generation (Natural Images) for NN Training

examples/scripts/ct_astra_3d_tv_admm.py

@@ -44,7 +44,7 @@
 tangle = snp.array(create_tangle_phantom(Nx, Ny, Nz))
 
 n_projection = 10  # number of projections
-angles = np.linspace(0, np.pi, n_projection)  # evenly spaced projection angles
+angles = np.linspace(0, np.pi, n_projection, endpoint=False)  # evenly spaced projection angles
 C = XRayTransform3D(
     tangle.shape, det_count=[Nz, max(Nx, Ny)], det_spacing=[1.0, 1.0], angles=angles
 )  # CT projection operator

examples/scripts/ct_astra_3d_tv_padmm.py

@@ -44,7 +44,7 @@
 tangle = snp.array(create_tangle_phantom(Nx, Ny, Nz))
 
 n_projection = 10  # number of projections
-angles = np.linspace(0, np.pi, n_projection)  # evenly spaced projection angles
+angles = np.linspace(0, np.pi, n_projection, endpoint=False)  # evenly spaced projection angles
 det_spacing = [1.0, 1.0]
 det_count = [Nz, max(Nx, Ny)]
 vectors = angle_to_vector(det_spacing, angles)
@@ -56,7 +56,7 @@
 y = C @ tangle  # sinogram
 
 
-"""
+r"""
 Set up problem and solver. We want to minimize the functional
 
   $$\mathrm{argmin}_{\mathbf{x}} \; (1/2) \| \mathbf{y} - C \mathbf{x}

examples/scripts/ct_astra_noreg_pcg.py

@@ -44,7 +44,7 @@
 Configure a CT projection operator and generate synthetic measurements.
 """
 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
+angles = np.linspace(0, np.pi, n_projection, endpoint=False)  # evenly spaced projection angles
 A = 1 / N * XRayTransform2D(x_gt.shape, N, 1.0, angles)  # CT projection operator
 y = A @ x_gt  # sinogram
 

examples/scripts/ct_astra_tv_admm.py

@@ -38,21 +38,22 @@
 """
 N = 512  # phantom size
 np.random.seed(1234)
-x_gt = discrete_phantom(Foam(size_range=[0.075, 0.0025], gap=1e-3, porosity=1), size=N)
-x_gt = snp.array(x_gt)  # convert to jax type
+x_gt = snp.array(discrete_phantom(Foam(size_range=[0.075, 0.0025], gap=1e-3, porosity=1), size=N))
 
 
 """
 Configure CT projection operator and generate synthetic measurements.
 """
 n_projection = 45  # number of projections
-angles = np.linspace(0, np.pi, n_projection)  # evenly spaced projection angles
-A = XRayTransform2D(x_gt.shape, N, 1.0, angles)  # CT projection operator
+angles = np.linspace(0, np.pi, n_projection, endpoint=False)  # evenly spaced projection angles
+det_count = int(N * 1.05 / np.sqrt(2.0))
+det_spacing = np.sqrt(2)
+A = XRayTransform2D(x_gt.shape, det_count, det_spacing, angles)  # CT projection operator
 y = A @ x_gt  # sinogram
 
 
 """
-Set up ADMM solver object.
+Set up problem functional and ADMM solver object.
 """
 λ = 2e0  # ℓ1 norm regularization parameter
 ρ = 5e0  # ADMM penalty parameter
@@ -65,9 +66,7 @@
 # which is used so that g(Cx) corresponds to isotropic TV.
 C = linop.FiniteDifference(input_shape=x_gt.shape, append=0)
 g = λ * functional.L21Norm()
-
 f = loss.SquaredL2Loss(y=y, A=A)
-
 x0 = snp.clip(A.fbp(y), 0, 1.0)
 
 solver = ADMM(

examples/scripts/ct_astra_weighted_tv_admm.py

@@ -35,7 +35,6 @@
 Create a ground truth image.
 """
 N = 512  # phantom size
-
 np.random.seed(0)
 x_gt = discrete_phantom(Soil(porosity=0.80), size=384)
 x_gt = np.ascontiguousarray(np.pad(x_gt, (64, 64)))
@@ -49,8 +48,7 @@
 n_projection = 360  # number of projections
 Io = 1e3  # source flux
 𝛼 = 1e-2  # attenuation coefficient
-
-angles = np.linspace(0, 2 * np.pi, n_projection)  # evenly spaced projection angles
+angles = np.linspace(0, 2 * np.pi, n_projection, endpoint=False)  # evenly spaced projection angles
 A = XRayTransform2D(x_gt.shape, N, 1.0, angles)  # CT projection operator
 y_c = A @ x_gt  # sinogram
 
@@ -99,13 +97,10 @@ def postprocess(x):
 # shown here).
 ρ = 2.5e3  # ADMM penalty parameter
 lambda_unweighted = 3e2  # regularization strength
-
 maxiter = 100  # number of ADMM iterations
 cg_tol = 1e-5  # CG relative tolerance
 cg_maxiter = 10  # maximum CG iterations per ADMM iteration
-
 f = loss.SquaredL2Loss(y=y, A=A)
-
 admm_unweighted = ADMM(
     f=f,
     g_list=[lambda_unweighted * functional.L21Norm()],
@@ -137,10 +132,8 @@ def postprocess(x):
 $I_0$ changes.
 """
 lambda_weighted = 5e1
-
 weights = snp.array(counts / Io)
 f = loss.SquaredL2Loss(y=y, A=A, W=linop.Diagonal(weights))
-
 admm_weighted = ADMM(
     f=f,
     g_list=[lambda_weighted * functional.L21Norm()],
@@ -151,6 +144,7 @@ def postprocess(x):
     subproblem_solver=LinearSubproblemSolver(cg_kwargs={"tol": cg_tol, "maxiter": cg_maxiter}),
     itstat_options={"display": True, "period": 10},
 )
+print()
 admm_weighted.solve()
 x_weighted = postprocess(admm_weighted.x)
 

examples/scripts/ct_astra_datagen_foam2.py → examples/scripts/ct_datagen_foam2.py

@@ -14,13 +14,17 @@
 """
 
 # isort: off
+import os
 import numpy as np
 
 import logging
 import ray
 
 ray.init(logging_level=logging.ERROR)  # need to call init before jax import: ray-project/ray#44087
 
+# Set an arbitrary processor count (only applies if GPU is not available).
+os.environ["XLA_FLAGS"] = "--xla_force_host_platform_device_count=8"
+
 from scico import plot
 from scico.flax.examples import load_ct_data
 

examples/scripts/ct_astra_modl_train_foam2.py → examples/scripts/ct_modl_train_foam2.py

@@ -5,8 +5,8 @@
 # with the package.
 
 r"""
-CT Training and Reconstructions with MoDL
-=========================================
+CT Training and Reconstruction with MoDL
+========================================
 
 This example demonstrates the training and application of a
 model-based deep learning (MoDL) architecture described in
@@ -65,7 +65,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.xray.astra import XRayTransform2D
+from scico.linop.xray import XRayTransform2D
 
 """
 Prepare parallel processing. Set an arbitrary processor count (only
@@ -89,16 +89,17 @@
 
 
 """
-Build CT projection operator.
+Build CT projection operator. Parameters are chosen so that the operator
+is equivalent to the one used to generate the training data.
 """
-angles = np.linspace(0, np.pi, n_projection)  # evenly spaced projection angles
+angles = np.linspace(0, np.pi, n_projection, endpoint=False)  # evenly spaced projection angles
 A = XRayTransform2D(
     input_shape=(N, N),
-    det_spacing=1,
-    det_count=N,
     angles=angles,
-)  # CT projection operator
-A = (1.0 / N) * A  # normalized
+    det_count=int(N * 1.05 / np.sqrt(2.0)),
+    dx=1.0 / np.sqrt(2),
+)
+A = (1.0 / N) * A  # normalize projection operator
 
 
 """