Update examples

docs/source/examples.rst

@@ -38,7 +38,6 @@ 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
@@ -137,6 +136,7 @@ Total Variation
    examples/ct_astra_3d_tv_admm
    examples/ct_astra_3d_tv_padmm
    examples/ct_astra_weighted_tv_admm
+   examples/ct_multi_tv_admm
    examples/ct_svmbir_tv_multi
    examples/deconv_circ_tv_admm
    examples/deconv_tv_admm
@@ -208,6 +208,7 @@ ADMM
    examples/ct_tv_admm
    examples/ct_astra_3d_tv_admm
    examples/ct_astra_weighted_tv_admm
+   examples/ct_multi_tv_admm
    examples/ct_svmbir_tv_multi
    examples/ct_svmbir_ppp_bm3d_admm_cg
    examples/ct_svmbir_ppp_bm3d_admm_prox

examples/jnb.py

@@ -33,11 +33,12 @@ def py_file_to_string(src):
             if import_seen:
                 # Once an import statement has been seen, break on encountering a line that
                 # is neither an import statement nor a newline, nor a component of an import
-                # statement extended over multiple lines, nor an os.environ statement, nor
-                # components of a try/except construction (note that handling of these final
-                # two cases is probably not very robust).
+                # statement extended over multiple lines, nor an os.environ statement, nor a
+                # ray.init statement, nor components of a try/except construction (note that
+                # handling of these final two cases is probably not very robust).
                 if not re.match(
-                    r"(^import|^from|^\n$|^\W+[^\W]|^\)$|^os.environ|^try:$|^except)", line
+                    r"(^import|^from|^\n$|^\W+[^\W]|^\)$|^os.environ|^ray.init|^try:$|^except)",
+                    line,
                 ):
                     lines.append(line)
                     break

examples/notebooks_requirements.txt

@@ -1,4 +1,6 @@
 -r examples-requirements.txt
+ipykernel
+ipywidgets
 nbformat
 nbconvert
 nb_conda_kernels

examples/scripts/README.rst

@@ -41,8 +41,6 @@ Computed Tomography
       CT Training and Reconstructions with UNet
    `ct_projector_comparison.py <ct_projector_comparison.py>`_
       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)
 
@@ -176,6 +174,8 @@ Total Variation
       3D TV-Regularized Sparse-View CT Reconstruction (Proximal ADMM Solver)
    `ct_astra_weighted_tv_admm.py <ct_astra_weighted_tv_admm.py>`_
       TV-Regularized Low-Dose CT Reconstruction
+   `ct_multi_tv_admm.py <ct_multi_tv_admm.py>`_
+      TV-Regularized Sparse-View CT Reconstruction (Multiple Projectors)
    `ct_svmbir_tv_multi.py <ct_svmbir_tv_multi.py>`_
       TV-Regularized CT Reconstruction (Multiple Algorithms)
    `deconv_circ_tv_admm.py <deconv_circ_tv_admm.py>`_
@@ -275,6 +275,8 @@ ADMM
       3D TV-Regularized Sparse-View CT Reconstruction (ADMM Solver)
    `ct_astra_weighted_tv_admm.py <ct_astra_weighted_tv_admm.py>`_
       TV-Regularized Low-Dose CT Reconstruction
+   `ct_multi_tv_admm.py <ct_multi_tv_admm.py>`_
+      TV-Regularized Sparse-View CT Reconstruction (Multiple Projectors)
    `ct_svmbir_tv_multi.py <ct_svmbir_tv_multi.py>`_
       TV-Regularized CT Reconstruction (Multiple Algorithms)
    `ct_svmbir_ppp_bm3d_admm_cg.py <ct_svmbir_ppp_bm3d_admm_cg.py>`_

examples/scripts/ct_multi_tv_admm.py

@@ -17,7 +17,8 @@
 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. The solution is computed and compared
-for all three 2D CT projectors available in scico.
+for all three 2D CT projectors available in scico, using a sinogram
+computed with the astra projector.
 """
 
 import numpy as np
@@ -37,44 +38,64 @@
 np.random.seed(1234)
 x_gt = snp.array(discrete_phantom(Foam(size_range=[0.075, 0.0025], gap=1e-3, porosity=1), size=N))
 
+det_count = N
+det_spacing = np.sqrt(2)
+
 
 """
 Define CT geometry and construct array of (approximately) equivalent projectors.
 """
 n_projection = 45  # number of projections
 angles = np.linspace(0, np.pi, n_projection)  # evenly spaced projection angles
 projectors = {
-    "astra": astra.XRayTransform2D(x_gt.shape, N, 1.0, angles - np.pi / 2.0),  # astra
-    "svmbir": svmbir.XRayTransform(x_gt.shape, 2 * np.pi - angles, N),  # svmbir
-    "scico": XRayTransform(Parallel2dProjector((N, N), angles, det_count=N)),  # scico
+    "astra": astra.XRayTransform2D(
+        x_gt.shape, det_count, det_spacing, angles - np.pi / 2.0
+    ),  # astra
+    "svmbir": svmbir.XRayTransform(
+        x_gt.shape, 2 * np.pi - angles, det_count, delta_pixel=1.0, delta_channel=det_spacing
+    ),  # svmbir
+    "scico": XRayTransform(
+        Parallel2dProjector((N, N), angles, det_count=det_count, dx=1 / det_spacing)
+    ),  # scico
 }
 
 
+"""
+Compute common sinogram using astra projector.
+"""
+A = projectors["astra"]
+noise = np.random.normal(size=(n_projection, det_count)).astype(np.float32)
+y = A @ x_gt + 2.0 * noise
+
+
+"""
+Construct initial solution for regularized problem.
+"""
+x0 = A.fbp(y)
+
+
 """
 Solve the same problem using the different projectors.
 """
 print(f"Solving on {device_info()}")
-y, x_rec, hist = {}, {}, {}
-noise = np.random.normal(size=(n_projection, N)).astype(np.float32)
-for p in ("astra", "svmbir", "scico"):
+x_rec, hist = {}, {}
+for p in projectors.keys():
     print(f"\nSolving with {p} projector")
-    A = projectors[p]
-    y[p] = A @ x_gt + 2.0 * noise  # sinogram
 
     # Set up ADMM solver object.
-    λ = 2e0  # L1 norm regularization parameter
-    ρ = 5e0  # ADMM penalty parameter
-    maxiter = 25  # number of ADMM iterations
+    λ = 2e1  # L1 norm regularization parameter
+    ρ = 1e3  # ADMM penalty parameter
+    maxiter = 100  # number of ADMM iterations
     cg_tol = 1e-4  # CG relative tolerance
-    cg_maxiter = 25  # maximum CG iterations per ADMM iteration
+    cg_maxiter = 50  # maximum CG iterations per ADMM iteration
 
     # The append=0 option makes the results of horizontal and vertical
     # finite differences the same shape, which is required for the L21Norm,
     # 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[p], A=A)
-    x0 = snp.clip(A.T(y[p]), 0, 1.0)
+    A = projectors[p]
+    f = loss.SquaredL2Loss(y=y, A=A)
 
     # Set up the solver.
     solver = ADMM(
@@ -91,15 +112,25 @@
     # Run the solver.
     solver.solve()
     hist[p] = solver.itstat_object.history(transpose=True)
-    x_rec[p] = snp.clip(solver.x, 0, 1.0)
+    x_rec[p] = solver.x
+
+    if p == "scico":
+        x_rec[p] = x_rec[p] * det_spacing  # to match ASTRA's scaling
+
+
+"""
+Compare reconstruction results.
+"""
+print("Reconstruction SNR:")
+for p in projectors.keys():
+    print(f"  {(p + ':'):7s}  {metric.snr(x_gt, x_rec[p]):5.2f} dB")
 
 
 """
-Compare sinograms.
+Display sinogram.
 """
-fig, ax = plot.subplots(nrows=3, ncols=1, figsize=(15, 10))
-for idx, name in enumerate(projectors.keys()):
-    plot.imview(y[name], title=f"{name} sinogram", cbar=None, fig=fig, ax=ax[idx])
+fig, ax = plot.subplots(nrows=1, ncols=1, figsize=(15, 3))
+plot.imview(y, title="sinogram", fig=fig, ax=ax)
 fig.show()
 
 
@@ -147,6 +178,8 @@
         fig=fig,
         ax=ax[n + 1],
     )
+for ax in ax:
+    ax.get_images()[0].set_clim(-0.1, 1.1)
 fig.show()