From f3ad2b99ef8a8d4259012703a9491d8c72fa4cc7 Mon Sep 17 00:00:00 2001 From: Brendt Wohlberg Date: Tue, 4 Jun 2024 14:19:15 -0600 Subject: [PATCH 01/33] Minor comment edits --- examples/scripts/ct_multi_cs_tv_admm.py | 2 +- examples/scripts/ct_multi_tv_admm.py | 2 +- examples/scripts/deconv_circ_tv_admm.py | 2 +- examples/scripts/deconv_tv_admm.py | 2 +- examples/scripts/sparsecode_conv_admm.py | 2 +- examples/scripts/sparsecode_conv_md_admm.py | 2 +- examples/scripts/video_rpca_admm.py | 2 +- 7 files changed, 7 insertions(+), 7 deletions(-) diff --git a/examples/scripts/ct_multi_cs_tv_admm.py b/examples/scripts/ct_multi_cs_tv_admm.py index 214098d4a..04d531569 100644 --- a/examples/scripts/ct_multi_cs_tv_admm.py +++ b/examples/scripts/ct_multi_cs_tv_admm.py @@ -68,7 +68,7 @@ print(f"\nSolving with {p} projector") # Set up ADMM solver object. - λ = 2e0 # L1 norm regularization parameter + λ = 2e0 # ℓ2,1 norm regularization parameter ρ = 5e0 # ADMM penalty parameter maxiter = 25 # number of ADMM iterations cg_tol = 1e-4 # CG relative tolerance diff --git a/examples/scripts/ct_multi_tv_admm.py b/examples/scripts/ct_multi_tv_admm.py index 084ad99ac..caa4270e3 100644 --- a/examples/scripts/ct_multi_tv_admm.py +++ b/examples/scripts/ct_multi_tv_admm.py @@ -62,7 +62,7 @@ y[p] = A @ x_gt + 2.0 * noise # sinogram # Set up ADMM solver object. - λ = 2e0 # L1 norm regularization parameter + λ = 2e0 # ℓ2,1 norm regularization parameter ρ = 5e0 # ADMM penalty parameter maxiter = 25 # number of ADMM iterations cg_tol = 1e-4 # CG relative tolerance diff --git a/examples/scripts/deconv_circ_tv_admm.py b/examples/scripts/deconv_circ_tv_admm.py index b2ba83202..90e898dc7 100644 --- a/examples/scripts/deconv_circ_tv_admm.py +++ b/examples/scripts/deconv_circ_tv_admm.py @@ -54,7 +54,7 @@ """ Set up an ADMM solver object. """ -λ = 2e-2 # L21 norm regularization parameter +λ = 2e-2 # ℓ2,1 norm regularization parameter ρ = 5e-1 # ADMM penalty parameter maxiter = 50 # number of ADMM iterations diff --git a/examples/scripts/deconv_tv_admm.py b/examples/scripts/deconv_tv_admm.py index 874b70e3c..8acb54d36 100644 --- a/examples/scripts/deconv_tv_admm.py +++ b/examples/scripts/deconv_tv_admm.py @@ -77,7 +77,7 @@ f = loss.SquaredL2Loss(y=y, A=C) # Penalty parameters must be accounted for in the gi functions, not as # additional inputs. -λ = 2.1e-2 # L21 norm regularization parameter +λ = 2.1e-2 # ℓ2,1 norm regularization parameter g = λ * functional.L21Norm() # The append=0 option makes the results of horizontal and vertical # finite differences the same shape, which is required for the L21Norm, diff --git a/examples/scripts/sparsecode_conv_admm.py b/examples/scripts/sparsecode_conv_admm.py index 611e7238b..0c46356a3 100644 --- a/examples/scripts/sparsecode_conv_admm.py +++ b/examples/scripts/sparsecode_conv_admm.py @@ -75,7 +75,7 @@ """ Set functional and solver parameters. """ -λ = 1e0 # l1-l2 norm regularization parameter +λ = 1e0 # ℓ1-ℓ2 norm regularization parameter ρ = 2e0 # ADMM penalty parameter maxiter = 200 # number of ADMM iterations diff --git a/examples/scripts/sparsecode_conv_md_admm.py b/examples/scripts/sparsecode_conv_md_admm.py index 9f70e9ce3..77d24f19b 100644 --- a/examples/scripts/sparsecode_conv_md_admm.py +++ b/examples/scripts/sparsecode_conv_md_admm.py @@ -96,7 +96,7 @@ """ Set functional and solver parameters. """ -λ = 1e0 # l1-l2 norm regularization parameter +λ = 1e0 # ℓ1-ℓ2 norm regularization parameter ρ0 = 1e0 # ADMM penalty parameters ρ1 = 3e0 maxiter = 200 # number of ADMM iterations diff --git a/examples/scripts/video_rpca_admm.py b/examples/scripts/video_rpca_admm.py index e6c71cb8e..ee8fa5c14 100644 --- a/examples/scripts/video_rpca_admm.py +++ b/examples/scripts/video_rpca_admm.py @@ -62,7 +62,7 @@ Set up an ADMM solver object. """ λ0 = 1e1 # nuclear norm regularization parameter -λ1 = 3e1 # l1 norm regularization parameter +λ1 = 3e1 # ℓ1 norm regularization parameter ρ0 = 2e1 # ADMM penalty parameter ρ1 = 2e1 # ADMM penalty parameter maxiter = 50 # number of ADMM iterations From c2efef529e6f4ca854d03d680229a2b98e48d1e0 Mon Sep 17 00:00:00 2001 From: Brendt Wohlberg Date: Wed, 5 Jun 2024 17:26:05 -0600 Subject: [PATCH 02/33] Fix missing typing import --- scico/numpy/util.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/scico/numpy/util.py b/scico/numpy/util.py index af1071b61..6917cd620 100644 --- a/scico/numpy/util.py +++ b/scico/numpy/util.py @@ -11,7 +11,7 @@ from __future__ import annotations from math import prod -from typing import Any, List, Optional, Tuple, Union +from typing import Any, List, Optional, Sequence, Tuple, Union import numpy as np From fd012fa3d20e22b7f077f8bccd0e3dee189fe7aa Mon Sep 17 00:00:00 2001 From: Brendt Wohlberg Date: Wed, 5 Jun 2024 17:29:16 -0600 Subject: [PATCH 03/33] Docs fixes --- scico/optimize/_pgm.py | 22 ++++++++++++---------- 1 file changed, 12 insertions(+), 10 deletions(-) diff --git a/scico/optimize/_pgm.py b/scico/optimize/_pgm.py index 8375c0b2d..0ae51623d 100644 --- a/scico/optimize/_pgm.py +++ b/scico/optimize/_pgm.py @@ -1,5 +1,5 @@ # -*- coding: utf-8 -*- -# Copyright (C) 2020-2023 by SCICO Developers +# Copyright (C) 2020-2024 by SCICO Developers # All rights reserved. BSD 3-clause License. # This file is part of the SCICO package. Details of the copyright and # user license can be found in the 'LICENSE' file distributed with the @@ -35,9 +35,9 @@ class PGM(Optimizer): Minimize a function of the form :math:`f(\mb{x}) + g(\mb{x})`, where :math:`f` and the :math:`g` are instances of :class:`.Functional`. - Uses helper :class:`StepSize` to provide an estimate of the Lipschitz + Uses :class:`.PGMStepSize` helper classes to estimate the Lipschitz constant :math:`L` of :math:`f`. The step size :math:`\alpha` is the - reciprocal of :math:`L`, i.e.: :math:`\alpha = 1 / L`. + reciprocal of :math:`L`, i.e. :math:`\alpha = 1 / L`. """ def __init__( @@ -52,12 +52,13 @@ def __init__( r""" Args: - f: Loss or Functional object with `grad` defined. - g: Instance of Functional with defined prox method. + f: Instance of :class:`.Loss` or :class:`.Functional` with + defined `grad` method. + g: Instance of :class:`.Functional` with defined prox method. L0: Initial estimate of Lipschitz constant of f. x0: Starting point for :math:`\mb{x}`. - step_size: helper :class:`StepSize` to estimate the Lipschitz - constant of f. + step_size: helper :class:`.PGMStepSize` to estimate the + Lipschitz constant of f. **kwargs: Additional optional parameters handled by initializer of base class :class:`.Optimizer`. """ @@ -173,11 +174,12 @@ def __init__( ): r""" Args: - f: Loss or Functional object with `grad` defined. - g: Instance of Functional with defined prox method. + f: Instance of :class:`.Loss` or :class:`.Functional` with + defined `grad` method. + g: Instance of :class:`.Functional` with defined prox method. L0: Initial estimate of Lipschitz constant of f. x0: Starting point for :math:`\mb{x}`. - step_size: helper :class:`StepSize` to estimate the Lipschitz + step_size: helper :class:`.PGMStepSize` to estimate the Lipschitz constant of f. **kwargs: Additional optional parameters handled by initializer of base class :class:`.Optimizer`. From d5820e3edac9b89e5a5eb8551b8eeba9e927db20 Mon Sep 17 00:00:00 2001 From: Brendt Wohlberg Date: Wed, 5 Jun 2024 18:44:21 -0600 Subject: [PATCH 04/33] Improve docs --- scico/optimize/_pgm.py | 39 +++++++++++++++++++++++++++++---------- 1 file changed, 29 insertions(+), 10 deletions(-) diff --git a/scico/optimize/_pgm.py b/scico/optimize/_pgm.py index 0ae51623d..53cf5ac4d 100644 --- a/scico/optimize/_pgm.py +++ b/scico/optimize/_pgm.py @@ -30,14 +30,35 @@ class PGM(Optimizer): - r"""Proximal Gradient Method (PGM) base class. + r"""Proximal gradient method (PGM) algorithm. - Minimize a function of the form :math:`f(\mb{x}) + g(\mb{x})`, where - :math:`f` and the :math:`g` are instances of :class:`.Functional`. + Minimize a functional of the form :math:`f(\mb{x}) + g(\mb{x})`, + where :math:`f` and the :math:`g` are instances of + :class:`.Functional`. Functional :math:`f` should be differentiable + and have a Lipschitz continuous derivative, and functional :math:`g` + should have a proximal operator defined. + + The step size :math:`\alpha` of the algorithm is defined in terms of + its reciprocal :math:`L`, i.e. :math:`\alpha = 1 / L`. The initial + value for this parameter, `L0`, is required to satisfy + + .. math:: + L_0 \geq K(\nabla f) \;, - Uses :class:`.PGMStepSize` helper classes to estimate the Lipschitz - constant :math:`L` of :math:`f`. The step size :math:`\alpha` is the - reciprocal of :math:`L`, i.e. :math:`\alpha = 1 / L`. + where :math:`K(\nabla f)` denotes the Lipschitz constant of the + gradient of :math:`f`. When `f` is an instance of + :class:`.SquaredL2Loss` with a :class:`.LinearOperator` `A`, + + .. math:: + K(\nabla f) = \lambda_{ \mathrm{max} }( A^H A ) = \| A \|_2^2 \;, + + where :math:`\lambda_{\mathrm{max}}(B)` denotes the largest + eigenvalue of :math:`B`. + + The evolution of the step size is controlled by auxiliary class + :class:`.PGMStepSize` and derived classes. The default + :class:`.PGMStepSize` simply sets :math:`L = L_0`, while the derived + classes implement a variety of adaptive strategies. """ def __init__( @@ -151,16 +172,14 @@ def step(self): class AcceleratedPGM(PGM): - r"""Accelerated Proximal Gradient Method (AcceleratedPGM) base class. - - Minimize a function of the form :math:`f(\mb{x}) + g(\mb{x})`. + r"""Accelerated proximal gradient method (APGM) algorithm. Minimize a function of the form :math:`f(\mb{x}) + g(\mb{x})`, where :math:`f` and the :math:`g` are instances of :class:`.Functional`. The accelerated form of PGM is also known as FISTA :cite:`beck-2009-fast`. - For documentation on inherited attributes, see :class:`.PGM`. + See :class:`.PGM` for more detailed documentation. """ def __init__( From d518af1a730e1da81be6ceda5f2e3f50e5f14a71 Mon Sep 17 00:00:00 2001 From: Brendt Wohlberg Date: Wed, 5 Jun 2024 18:49:16 -0600 Subject: [PATCH 05/33] Docs fixes --- scico/optimize/_pgm.py | 16 +++++++++------- 1 file changed, 9 insertions(+), 7 deletions(-) diff --git a/scico/optimize/_pgm.py b/scico/optimize/_pgm.py index 53cf5ac4d..6188f7480 100644 --- a/scico/optimize/_pgm.py +++ b/scico/optimize/_pgm.py @@ -76,10 +76,11 @@ def __init__( f: Instance of :class:`.Loss` or :class:`.Functional` with defined `grad` method. g: Instance of :class:`.Functional` with defined prox method. - L0: Initial estimate of Lipschitz constant of f. + L0: Initial estimate of Lipschitz constant of gradient of `f`. x0: Starting point for :math:`\mb{x}`. - step_size: helper :class:`.PGMStepSize` to estimate the - Lipschitz constant of f. + step_size: Instance of an auxiliary class of type + :class:`.PGMStepSize` determining the evolution of the + algorithm step size. **kwargs: Additional optional parameters handled by initializer of base class :class:`.Optimizer`. """ @@ -97,7 +98,7 @@ def __init__( step_size = PGMStepSize() self.step_size: PGMStepSize = step_size self.step_size.internal_init(self) - self.L: float = L0 # reciprocal of step size (estimate of Lipschitz constant of f) + self.L: float = L0 # reciprocal of step size (estimate of Lipschitz constant of ∇f) self.fixed_point_residual = snp.inf def x_step(v: Union[Array, BlockArray], L: float) -> Union[Array, BlockArray]: @@ -196,10 +197,11 @@ def __init__( f: Instance of :class:`.Loss` or :class:`.Functional` with defined `grad` method. g: Instance of :class:`.Functional` with defined prox method. - L0: Initial estimate of Lipschitz constant of f. + L0: Initial estimate of Lipschitz constant of gradient of `f`. x0: Starting point for :math:`\mb{x}`. - step_size: helper :class:`.PGMStepSize` to estimate the Lipschitz - constant of f. + step_size: Instance of an auxiliary class of type + :class:`.PGMStepSize` determining the evolution of the + algorithm step size. **kwargs: Additional optional parameters handled by initializer of base class :class:`.Optimizer`. """ From e01f1b6ac31876205fa2902200a73eeace5749b6 Mon Sep 17 00:00:00 2001 From: Brendt Wohlberg Date: Wed, 5 Jun 2024 19:13:06 -0600 Subject: [PATCH 06/33] Add missing dependency (due to recent changes in other packages) for building notebooks --- examples/notebooks_requirements.txt | 1 + 1 file changed, 1 insertion(+) diff --git a/examples/notebooks_requirements.txt b/examples/notebooks_requirements.txt index bcb9e03c9..dab8d40df 100644 --- a/examples/notebooks_requirements.txt +++ b/examples/notebooks_requirements.txt @@ -1,4 +1,5 @@ -r examples-requirements.txt +ipykernel nbformat nbconvert nb_conda_kernels From 9c4131bec1f5aaaf69b55a4aac9fafe062bfef97 Mon Sep 17 00:00:00 2001 From: Brendt Wohlberg Date: Wed, 5 Jun 2024 19:54:54 -0600 Subject: [PATCH 07/33] Update submodule --- data | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/data b/data index b085cef32..5473ae209 160000 --- a/data +++ b/data @@ -1 +1 @@ -Subproject commit b085cef32b1866cb87f5d10734355c665491e400 +Subproject commit 5473ae209af24fc3b6fdaedc24e63903737856b5 From e40615a05df676e3001dbed77b54df3442098a0b Mon Sep 17 00:00:00 2001 From: Brendt Wohlberg Date: Wed, 5 Jun 2024 20:12:49 -0600 Subject: [PATCH 08/33] Fix script to notebook conversion bug --- examples/jnb.py | 9 +++++---- 1 file changed, 5 insertions(+), 4 deletions(-) diff --git a/examples/jnb.py b/examples/jnb.py index 5b8ead6ee..9135add1d 100644 --- a/examples/jnb.py +++ b/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 From cbec0cf991ebcf4f313aab27f609dac1e19245e9 Mon Sep 17 00:00:00 2001 From: Brendt Wohlberg Date: Wed, 5 Jun 2024 20:20:27 -0600 Subject: [PATCH 09/33] Update submodule --- data | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/data b/data index 5473ae209..881c536f9 160000 --- a/data +++ b/data @@ -1 +1 @@ -Subproject commit 5473ae209af24fc3b6fdaedc24e63903737856b5 +Subproject commit 881c536f96ba7a4750dd2e7ac5232b791c0b64a5 From 5124e627f8bd4dc9dfa8df7ac28fa81c78170bb5 Mon Sep 17 00:00:00 2001 From: Brendt Wohlberg Date: Thu, 6 Jun 2024 08:19:40 -0600 Subject: [PATCH 10/33] Update submodule --- data | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/data b/data index 881c536f9..04d2cbdd3 160000 --- a/data +++ b/data @@ -1 +1 @@ -Subproject commit 881c536f96ba7a4750dd2e7ac5232b791c0b64a5 +Subproject commit 04d2cbdd3c6ea3944422bd89a8514fd7f629c20e From b64f97745e47ef19de4cab2557a2280d3183bd84 Mon Sep 17 00:00:00 2001 From: Brendt Wohlberg Date: Thu, 6 Jun 2024 08:39:53 -0600 Subject: [PATCH 11/33] Update submodule --- data | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/data b/data index 04d2cbdd3..b3bbcf91e 160000 --- a/data +++ b/data @@ -1 +1 @@ -Subproject commit 04d2cbdd3c6ea3944422bd89a8514fd7f629c20e +Subproject commit b3bbcf91e35591d1567f578e0742e9b71f64ac39 From a46200cfcb9b940d494acc62164195712c71e0c0 Mon Sep 17 00:00:00 2001 From: Brendt Wohlberg Date: Thu, 6 Jun 2024 11:59:10 -0600 Subject: [PATCH 12/33] Update submodule --- data | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/data b/data index b3bbcf91e..f31b7ef7b 160000 --- a/data +++ b/data @@ -1 +1 @@ -Subproject commit b3bbcf91e35591d1567f578e0742e9b71f64ac39 +Subproject commit f31b7ef7b81cde59ddb7404cf129d7b4c63dcae8 From 9c5d6a7fb7ef2a9367e3dbcff1979009785a7d54 Mon Sep 17 00:00:00 2001 From: Brendt Wohlberg Date: Fri, 7 Jun 2024 08:18:42 -0600 Subject: [PATCH 13/33] Add missing dependency for building notebooks --- examples/notebooks_requirements.txt | 1 + 1 file changed, 1 insertion(+) diff --git a/examples/notebooks_requirements.txt b/examples/notebooks_requirements.txt index dab8d40df..644a4db21 100644 --- a/examples/notebooks_requirements.txt +++ b/examples/notebooks_requirements.txt @@ -1,5 +1,6 @@ -r examples-requirements.txt ipykernel +ipywidgets nbformat nbconvert nb_conda_kernels From 4c3a195fa8d3a1ad230f9752f0130d6e9e617861 Mon Sep 17 00:00:00 2001 From: Brendt Wohlberg Date: Wed, 12 Jun 2024 19:10:34 -0600 Subject: [PATCH 14/33] Update submodule --- data | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/data b/data index f31b7ef7b..ff536809a 160000 --- a/data +++ b/data @@ -1 +1 @@ -Subproject commit f31b7ef7b81cde59ddb7404cf129d7b4c63dcae8 +Subproject commit ff536809a7cc8b0086644e27a8e9f955e8553ba8 From 5e38ca8b68c8eb675736eb914dd3446173e8428e Mon Sep 17 00:00:00 2001 From: Brendt Wohlberg Date: Thu, 13 Jun 2024 08:48:44 -0600 Subject: [PATCH 15/33] Improve example script docstring --- examples/scripts/deconv_microscopy_allchn_tv_admm.py | 8 ++++---- examples/scripts/deconv_microscopy_tv_admm.py | 8 ++++---- 2 files changed, 8 insertions(+), 8 deletions(-) diff --git a/examples/scripts/deconv_microscopy_allchn_tv_admm.py b/examples/scripts/deconv_microscopy_allchn_tv_admm.py index 13ec27b87..2606a7be1 100644 --- a/examples/scripts/deconv_microscopy_allchn_tv_admm.py +++ b/examples/scripts/deconv_microscopy_allchn_tv_admm.py @@ -37,10 +37,10 @@ from scico.optimize.admm import ADMM, CircularConvolveSolver """ -Get and preprocess data. We downsample the data for the for purposes of -the example. Reducing the downsampling rate will make the example slower -and more memory-intensive. To run this example on a GPU it may be -necessary to set environment variables +Get and preprocess data. The data is downsampled to limit the memory +requirements and run time of the example. Reducing the downsampling rate +will make the example slower and more memory-intensive. To run this +example on a GPU it may be necessary to set environment variables `XLA_PYTHON_CLIENT_ALLOCATOR=platform` and `XLA_PYTHON_CLIENT_PREALLOCATE=false`. If your GPU does not have enough memory, you can try setting the environment variable diff --git a/examples/scripts/deconv_microscopy_tv_admm.py b/examples/scripts/deconv_microscopy_tv_admm.py index eee4d8f42..5c9180b60 100644 --- a/examples/scripts/deconv_microscopy_tv_admm.py +++ b/examples/scripts/deconv_microscopy_tv_admm.py @@ -34,10 +34,10 @@ from scico.optimize.admm import ADMM, CircularConvolveSolver """ -Get and preprocess data. We downsample the data for the for purposes of -the example. Reducing the downsampling rate will make the example slower -and more memory-intensive. To run this example on a GPU it may be -necessary to set environment variables +Get and preprocess data. The data is downsampled to limit the memory +requirements and run time of the example. Reducing the downsampling rate +will make the example slower and more memory-intensive. To run this +example on a GPU it may be necessary to set environment variables `XLA_PYTHON_CLIENT_ALLOCATOR=platform` and `XLA_PYTHON_CLIENT_PREALLOCATE=false`. If your GPU does not have enough memory, you can try setting the environment variable From abd683da3edb53f748ff4e729e3a5a5b662a1313 Mon Sep 17 00:00:00 2001 From: Brendt Wohlberg Date: Thu, 13 Jun 2024 08:59:33 -0600 Subject: [PATCH 16/33] Update submodule --- data | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/data b/data index ff536809a..390a272a4 160000 --- a/data +++ b/data @@ -1 +1 @@ -Subproject commit ff536809a7cc8b0086644e27a8e9f955e8553ba8 +Subproject commit 390a272a4b96afb29ae3b0b9f8abd7b4d36caab9 From 80ac6289d6d70a11b399685ade251407b738634c Mon Sep 17 00:00:00 2001 From: Brendt Wohlberg Date: Thu, 13 Jun 2024 09:05:00 -0600 Subject: [PATCH 17/33] Avoid ray complaints --- examples/scripts/deconv_microscopy_allchn_tv_admm.py | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/examples/scripts/deconv_microscopy_allchn_tv_admm.py b/examples/scripts/deconv_microscopy_allchn_tv_admm.py index 2606a7be1..ac8ee583d 100644 --- a/examples/scripts/deconv_microscopy_allchn_tv_admm.py +++ b/examples/scripts/deconv_microscopy_allchn_tv_admm.py @@ -28,9 +28,14 @@ non-negativity constraint, and $\mathbf{x}$ is the desired image. """ +# isort: off 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 + import scico.numpy as snp from scico import functional, linop, loss, plot from scico.examples import downsample_volume, epfl_deconv_data, tile_volume_slices From 822cef883fd4c6cc1cebf5ca9bdac1faa8bc3475 Mon Sep 17 00:00:00 2001 From: Brendt Wohlberg Date: Thu, 13 Jun 2024 09:30:11 -0600 Subject: [PATCH 18/33] Remove second ray.init call --- examples/scripts/deconv_microscopy_allchn_tv_admm.py | 6 ++---- 1 file changed, 2 insertions(+), 4 deletions(-) diff --git a/examples/scripts/deconv_microscopy_allchn_tv_admm.py b/examples/scripts/deconv_microscopy_allchn_tv_admm.py index ac8ee583d..ccd8c9471 100644 --- a/examples/scripts/deconv_microscopy_allchn_tv_admm.py +++ b/examples/scripts/deconv_microscopy_allchn_tv_admm.py @@ -86,11 +86,9 @@ """ -Initialize ray, determine available computing resources, and put large arrays -in object store. +Determine available computing resources, and put large arrays in ray +object store. """ -ray.init() - ngpu = 0 ar = ray.available_resources() ncpu = max(int(ar["CPU"]) // 3, 1) From f72adbebe4a05b86a85dcccf8ab252cad038bb3a Mon Sep 17 00:00:00 2001 From: Brendt Wohlberg Date: Mon, 17 Jun 2024 16:47:34 -0600 Subject: [PATCH 19/33] Update submodule --- data | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/data b/data index 390a272a4..196ee56ca 160000 --- a/data +++ b/data @@ -1 +1 @@ -Subproject commit 390a272a4b96afb29ae3b0b9f8abd7b4d36caab9 +Subproject commit 196ee56ca70eda3a96962d8e472951082c5bae6b From 6bee347dab16333daf374430ca44c7a666312249 Mon Sep 17 00:00:00 2001 From: Brendt Wohlberg Date: Tue, 18 Jun 2024 15:04:20 -0600 Subject: [PATCH 20/33] Update submodule --- data | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/data b/data index 196ee56ca..25ed52ba8 160000 --- a/data +++ b/data @@ -1 +1 @@ -Subproject commit 196ee56ca70eda3a96962d8e472951082c5bae6b +Subproject commit 25ed52ba8fe47ba12e3d62ef9d73982f7a266a62 From 020214a3cf0e616e87c62098f2ee13d6977d7fa4 Mon Sep 17 00:00:00 2001 From: Brendt Wohlberg Date: Tue, 18 Jun 2024 15:10:56 -0600 Subject: [PATCH 21/33] Update submodule --- data | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/data b/data index 196ee56ca..c55c26df3 160000 --- a/data +++ b/data @@ -1 +1 @@ -Subproject commit 196ee56ca70eda3a96962d8e472951082c5bae6b +Subproject commit c55c26df3b69ff4d331384a6a07460b07a2f5932 From 503b6d540a564076638803206dafe4b87c2f6e5d Mon Sep 17 00:00:00 2001 From: Brendt Wohlberg Date: Wed, 19 Jun 2024 09:18:48 -0600 Subject: [PATCH 22/33] Replace deprecated matplotlib method --- scico/plot.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/scico/plot.py b/scico/plot.py index 1cdd1adc7..6c0375be1 100644 --- a/scico/plot.py +++ b/scico/plot.py @@ -1,5 +1,5 @@ # -*- coding: utf-8 -*- -# Copyright (C) 2020-2022 by SCICO Developers +# Copyright (C) 2020-2024 by SCICO Developers # All rights reserved. BSD 3-clause License. # This file is part of the SCICO package. Details of the copyright and # user license can be found in the 'LICENSE' file distributed with the @@ -360,7 +360,7 @@ def surf( # https://stackoverflow.com/a/35221116 if ax.name != "3d": ax.remove() - ax = fig.add_subplot(*ax.get_geometry(), projection="3d") + ax = fig.add_subplot(ax.get_subplotspec(), projection="3d") if elev is not None or azim is not None: ax.view_init(elev=elev, azim=azim) From e7a27a960591009fc345d0494be0e209f89c3246 Mon Sep 17 00:00:00 2001 From: Brendt Wohlberg Date: Tue, 9 Jul 2024 14:25:27 -0600 Subject: [PATCH 23/33] Update submodule --- data | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/data b/data index c55c26df3..07fb707e3 160000 --- a/data +++ b/data @@ -1 +1 @@ -Subproject commit c55c26df3b69ff4d331384a6a07460b07a2f5932 +Subproject commit 07fb707e3c12f8a885778a4a06efafd93535c4b9 From 2829fcbd8c3183f970cb3fdc2606766f9b3b18a1 Mon Sep 17 00:00:00 2001 From: Brendt Wohlberg Date: Tue, 9 Jul 2024 14:48:04 -0600 Subject: [PATCH 24/33] Update submodule --- data | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/data b/data index 07fb707e3..ee18f0b77 160000 --- a/data +++ b/data @@ -1 +1 @@ -Subproject commit 07fb707e3c12f8a885778a4a06efafd93535c4b9 +Subproject commit ee18f0b77da681694ee941503913c2da8aefc1ed From 0eb31fdfb6638f2793e7dea10c87163690ea9c77 Mon Sep 17 00:00:00 2001 From: Brendt Wohlberg Date: Tue, 9 Jul 2024 16:27:24 -0600 Subject: [PATCH 25/33] Remove largely redundant example --- examples/scripts/ct_multi_tv_admm.py | 153 --------------------------- examples/scripts/index.rst | 3 +- 2 files changed, 2 insertions(+), 154 deletions(-) delete mode 100644 examples/scripts/ct_multi_tv_admm.py diff --git a/examples/scripts/ct_multi_tv_admm.py b/examples/scripts/ct_multi_tv_admm.py deleted file mode 100644 index caa4270e3..000000000 --- a/examples/scripts/ct_multi_tv_admm.py +++ /dev/null @@ -1,153 +0,0 @@ -#!/usr/bin/env python -# -*- coding: utf-8 -*- -# This file is part of the SCICO package. Details of the copyright -# and user license can be found in the 'LICENSE.txt' file distributed -# with the package. - -r""" -TV-Regularized Sparse-View CT Reconstruction (Multiple Projectors) -================================================================== - -This example demonstrates solution of a sparse-view CT reconstruction -problem with isotropic total variation (TV) regularization - - $$\mathrm{argmin}_{\mathbf{x}} \; (1/2) \| \mathbf{y} - A \mathbf{x} - \|_2^2 + \lambda \| C \mathbf{x} \|_{2,1} \;,$$ - -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. -""" - -import numpy as np - -from xdesign import Foam, discrete_phantom - -import scico.numpy as snp -from scico import functional, linop, loss, metric, plot -from scico.linop.xray import Parallel2dProjector, XRayTransform, astra, svmbir -from scico.optimize.admm import ADMM, LinearSubproblemSolver -from scico.util import device_info - -""" -Create a ground truth image. -""" -N = 512 # phantom size -np.random.seed(1234) -x_gt = snp.array(discrete_phantom(Foam(size_range=[0.075, 0.0025], gap=1e-3, porosity=1), size=N)) - - -""" -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 -} - - -""" -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"): - print(f"\nSolving with {p} projector") - A = projectors[p] - y[p] = A @ x_gt + 2.0 * noise # sinogram - - # Set up ADMM solver object. - λ = 2e0 # ℓ2,1 norm regularization parameter - ρ = 5e0 # ADMM penalty parameter - maxiter = 25 # number of ADMM iterations - cg_tol = 1e-4 # CG relative tolerance - cg_maxiter = 25 # 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) - - # Set up the solver. - solver = ADMM( - f=f, - g_list=[g], - C_list=[C], - rho_list=[ρ], - x0=x0, - maxiter=maxiter, - subproblem_solver=LinearSubproblemSolver(cg_kwargs={"tol": cg_tol, "maxiter": cg_maxiter}), - itstat_options={"display": True, "period": 5}, - ) - - # Run the solver. - solver.solve() - hist[p] = solver.itstat_object.history(transpose=True) - x_rec[p] = snp.clip(solver.x, 0, 1.0) - - -""" -Compare sinograms. -""" -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.show() - - -""" -Plot convergence statistics. -""" -fig, ax = plot.subplots(nrows=1, ncols=3, figsize=(12, 5)) -plot.plot( - np.vstack([hist[p].Objective for p in projectors.keys()]).T, - title="Objective function", - xlbl="Iteration", - ylbl="Functional value", - lgnd=projectors.keys(), - fig=fig, - ax=ax[0], -) -plot.plot( - np.vstack([hist[p].Prml_Rsdl for p in projectors.keys()]).T, - ptyp="semilogy", - title="Primal Residual", - xlbl="Iteration", - fig=fig, - ax=ax[1], -) -plot.plot( - np.vstack([hist[p].Dual_Rsdl for p in projectors.keys()]).T, - ptyp="semilogy", - title="Dual Residual", - xlbl="Iteration", - fig=fig, - ax=ax[2], -) -fig.show() - - -""" -Show the recovered images. -""" -fig, ax = plot.subplots(nrows=1, ncols=4, figsize=(15, 5)) -plot.imview(x_gt, title="Ground truth", fig=fig, ax=ax[0]) -for n, p in enumerate(projectors.keys()): - plot.imview( - x_rec[p], - title="%s SNR: %.2f (dB)" % (p, metric.snr(x_gt, x_rec[p])), - fig=fig, - ax=ax[n + 1], - ) -fig.show() - - -input("\nWaiting for input to close figures and exit") diff --git a/examples/scripts/index.rst b/examples/scripts/index.rst index bbc697a95..7a4f83abc 100644 --- a/examples/scripts/index.rst +++ b/examples/scripts/index.rst @@ -26,7 +26,6 @@ Computed Tomography - ct_astra_unet_train_foam2.py - ct_projector_comparison.py - ct_multi_cs_tv_admm.py - - ct_multi_tv_admm.py Deconvolution ^^^^^^^^^^^^^ @@ -106,6 +105,7 @@ Total Variation - ct_astra_3d_tv_admm.py - ct_astra_3d_tv_padmm.py - ct_astra_weighted_tv_admm.py + - ct_multi_cs_tv_admm.py - ct_svmbir_tv_multi.py - deconv_circ_tv_admm.py - deconv_tv_admm.py @@ -165,6 +165,7 @@ ADMM - ct_tv_admm.py - ct_astra_3d_tv_admm.py - ct_astra_weighted_tv_admm.py + - ct_multi_cs_tv_admm.py - ct_svmbir_tv_multi.py - ct_svmbir_ppp_bm3d_admm_cg.py - ct_svmbir_ppp_bm3d_admm_prox.py From ab90d40a739f12806a48f0ac19c429737d700b39 Mon Sep 17 00:00:00 2001 From: Brendt Wohlberg Date: Tue, 9 Jul 2024 16:28:54 -0600 Subject: [PATCH 26/33] Update example index --- docs/source/examples.rst | 3 ++- examples/scripts/README.rst | 6 ++++-- 2 files changed, 6 insertions(+), 3 deletions(-) diff --git a/docs/source/examples.rst b/docs/source/examples.rst index da89c8289..a5aa278e3 100644 --- a/docs/source/examples.rst +++ b/docs/source/examples.rst @@ -39,7 +39,6 @@ Computed Tomography 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_cs_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_cs_tv_admm examples/ct_svmbir_tv_multi examples/ct_svmbir_ppp_bm3d_admm_cg examples/ct_svmbir_ppp_bm3d_admm_prox diff --git a/examples/scripts/README.rst b/examples/scripts/README.rst index 95d6aed8c..26a8f0a73 100644 --- a/examples/scripts/README.rst +++ b/examples/scripts/README.rst @@ -43,8 +43,6 @@ Computed Tomography X-ray Transform Comparison `ct_multi_cs_tv_admm.py `_ TV-Regularized Sparse-View CT Reconstruction (Multiple Projectors, Common Sinogram) - `ct_multi_tv_admm.py `_ - TV-Regularized Sparse-View CT Reconstruction (Multiple Projectors) Deconvolution ^^^^^^^^^^^^^ @@ -176,6 +174,8 @@ Total Variation 3D TV-Regularized Sparse-View CT Reconstruction (Proximal ADMM Solver) `ct_astra_weighted_tv_admm.py `_ TV-Regularized Low-Dose CT Reconstruction + `ct_multi_cs_tv_admm.py `_ + TV-Regularized Sparse-View CT Reconstruction (Multiple Projectors, Common Sinogram) `ct_svmbir_tv_multi.py `_ TV-Regularized CT Reconstruction (Multiple Algorithms) `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 `_ TV-Regularized Low-Dose CT Reconstruction + `ct_multi_cs_tv_admm.py `_ + TV-Regularized Sparse-View CT Reconstruction (Multiple Projectors, Common Sinogram) `ct_svmbir_tv_multi.py `_ TV-Regularized CT Reconstruction (Multiple Algorithms) `ct_svmbir_ppp_bm3d_admm_cg.py `_ From be7fb76ed994c6482d8720c477cda110950df9f6 Mon Sep 17 00:00:00 2001 From: Brendt Wohlberg Date: Tue, 9 Jul 2024 16:31:11 -0600 Subject: [PATCH 27/33] Rename example script --- examples/scripts/ct_multi_cs_tv_admm.py | 186 ------------------------ examples/scripts/index.rst | 6 +- 2 files changed, 3 insertions(+), 189 deletions(-) delete mode 100644 examples/scripts/ct_multi_cs_tv_admm.py diff --git a/examples/scripts/ct_multi_cs_tv_admm.py b/examples/scripts/ct_multi_cs_tv_admm.py deleted file mode 100644 index cb8ddc6b6..000000000 --- a/examples/scripts/ct_multi_cs_tv_admm.py +++ /dev/null @@ -1,186 +0,0 @@ -#!/usr/bin/env python -# -*- coding: utf-8 -*- -# This file is part of the SCICO package. Details of the copyright -# and user license can be found in the 'LICENSE.txt' file distributed -# with the package. - -r""" -TV-Regularized Sparse-View CT Reconstruction (Multiple Projectors, Common Sinogram) -=================================================================================== - -This example demonstrates solution of a sparse-view CT reconstruction -problem with isotropic total variation (TV) regularization - - $$\mathrm{argmin}_{\mathbf{x}} \; (1/2) \| \mathbf{y} - A \mathbf{x} - \|_2^2 + \lambda \| C \mathbf{x} \|_{2,1} \;,$$ - -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, using a sinogram -computed with the svmbir projector. -""" - -import numpy as np - -from xdesign import Foam, discrete_phantom - -import scico.numpy as snp -from scico import functional, linop, loss, metric, plot -from scico.linop.xray import Parallel2dProjector, XRayTransform, astra, svmbir -from scico.optimize.admm import ADMM, LinearSubproblemSolver -from scico.util import device_info - -""" -Create a ground truth image. -""" -N = 512 # phantom size -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, 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 svmbir 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()}") -x_rec, hist = {}, {} -for p in projectors.keys(): - print(f"\nSolving with {p} projector") - - # Set up ADMM solver object. - λ = 2e1 # L1 norm regularization parameter - ρ = 1e3 # ADMM penalty parameter - maxiter = 100 # number of ADMM iterations - cg_tol = 1e-4 # CG relative tolerance - 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() - A = projectors[p] - f = loss.SquaredL2Loss(y=y, A=A) - - # Set up the solver. - solver = ADMM( - f=f, - g_list=[g], - C_list=[C], - rho_list=[ρ], - x0=x0, - maxiter=maxiter, - subproblem_solver=LinearSubproblemSolver(cg_kwargs={"tol": cg_tol, "maxiter": cg_maxiter}), - itstat_options={"display": True, "period": 5}, - ) - - # Run the solver. - solver.solve() - hist[p] = solver.itstat_object.history(transpose=True) - 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") - - -""" -Display sinogram. -""" -fig, ax = plot.subplots(nrows=1, ncols=1, figsize=(15, 3)) -plot.imview(y, title="sinogram", fig=fig, ax=ax) -fig.show() - - -""" -Plot convergence statistics. -""" -fig, ax = plot.subplots(nrows=1, ncols=3, figsize=(12, 5)) -plot.plot( - np.vstack([hist[p].Objective for p in projectors.keys()]).T, - title="Objective function", - xlbl="Iteration", - ylbl="Functional value", - lgnd=projectors.keys(), - fig=fig, - ax=ax[0], -) -plot.plot( - np.vstack([hist[p].Prml_Rsdl for p in projectors.keys()]).T, - ptyp="semilogy", - title="Primal Residual", - xlbl="Iteration", - fig=fig, - ax=ax[1], -) -plot.plot( - np.vstack([hist[p].Dual_Rsdl for p in projectors.keys()]).T, - ptyp="semilogy", - title="Dual Residual", - xlbl="Iteration", - fig=fig, - ax=ax[2], -) -fig.show() - - -""" -Show the recovered images. -""" -fig, ax = plot.subplots(nrows=1, ncols=4, figsize=(15, 5)) -plot.imview(x_gt, title="Ground truth", fig=fig, ax=ax[0]) -for n, p in enumerate(projectors.keys()): - plot.imview( - x_rec[p], - title="%s SNR: %.2f (dB)" % (p, metric.snr(x_gt, x_rec[p])), - fig=fig, - ax=ax[n + 1], - ) -for ax in ax: - ax.get_images()[0].set_clim(-0.1, 1.1) -fig.show() - - -input("\nWaiting for input to close figures and exit") diff --git a/examples/scripts/index.rst b/examples/scripts/index.rst index 7a4f83abc..8100f06e6 100644 --- a/examples/scripts/index.rst +++ b/examples/scripts/index.rst @@ -25,7 +25,7 @@ Computed Tomography - ct_astra_odp_train_foam2.py - ct_astra_unet_train_foam2.py - ct_projector_comparison.py - - ct_multi_cs_tv_admm.py + - ct_multi_tv_admm.py Deconvolution ^^^^^^^^^^^^^ @@ -105,7 +105,7 @@ Total Variation - ct_astra_3d_tv_admm.py - ct_astra_3d_tv_padmm.py - ct_astra_weighted_tv_admm.py - - ct_multi_cs_tv_admm.py + - ct_multi_tv_admm.py - ct_svmbir_tv_multi.py - deconv_circ_tv_admm.py - deconv_tv_admm.py @@ -165,7 +165,7 @@ ADMM - ct_tv_admm.py - ct_astra_3d_tv_admm.py - ct_astra_weighted_tv_admm.py - - ct_multi_cs_tv_admm.py + - ct_multi_tv_admm.py - ct_svmbir_tv_multi.py - ct_svmbir_ppp_bm3d_admm_cg.py - ct_svmbir_ppp_bm3d_admm_prox.py From d4285d1ba9566e215cde87c0b6aa1ae1271292c8 Mon Sep 17 00:00:00 2001 From: Brendt Wohlberg Date: Tue, 9 Jul 2024 16:33:04 -0600 Subject: [PATCH 28/33] Edit script title and fix docs --- examples/scripts/ct_multi_tv_admm.py | 186 +++++++++++++++++++++++++++ 1 file changed, 186 insertions(+) create mode 100644 examples/scripts/ct_multi_tv_admm.py diff --git a/examples/scripts/ct_multi_tv_admm.py b/examples/scripts/ct_multi_tv_admm.py new file mode 100644 index 000000000..87fc86865 --- /dev/null +++ b/examples/scripts/ct_multi_tv_admm.py @@ -0,0 +1,186 @@ +#!/usr/bin/env python +# -*- coding: utf-8 -*- +# This file is part of the SCICO package. Details of the copyright +# and user license can be found in the 'LICENSE.txt' file distributed +# with the package. + +r""" +TV-Regularized Sparse-View CT Reconstruction (Multiple Projectors) +================================================================== + +This example demonstrates solution of a sparse-view CT reconstruction +problem with isotropic total variation (TV) regularization + + $$\mathrm{argmin}_{\mathbf{x}} \; (1/2) \| \mathbf{y} - A \mathbf{x} + \|_2^2 + \lambda \| C \mathbf{x} \|_{2,1} \;,$$ + +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, using a sinogram +computed with the astra projector. +""" + +import numpy as np + +from xdesign import Foam, discrete_phantom + +import scico.numpy as snp +from scico import functional, linop, loss, metric, plot +from scico.linop.xray import Parallel2dProjector, XRayTransform, astra, svmbir +from scico.optimize.admm import ADMM, LinearSubproblemSolver +from scico.util import device_info + +""" +Create a ground truth image. +""" +N = 512 # phantom size +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, 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()}") +x_rec, hist = {}, {} +for p in projectors.keys(): + print(f"\nSolving with {p} projector") + + # Set up ADMM solver object. + λ = 2e1 # L1 norm regularization parameter + ρ = 1e3 # ADMM penalty parameter + maxiter = 100 # number of ADMM iterations + cg_tol = 1e-4 # CG relative tolerance + 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() + A = projectors[p] + f = loss.SquaredL2Loss(y=y, A=A) + + # Set up the solver. + solver = ADMM( + f=f, + g_list=[g], + C_list=[C], + rho_list=[ρ], + x0=x0, + maxiter=maxiter, + subproblem_solver=LinearSubproblemSolver(cg_kwargs={"tol": cg_tol, "maxiter": cg_maxiter}), + itstat_options={"display": True, "period": 5}, + ) + + # Run the solver. + solver.solve() + hist[p] = solver.itstat_object.history(transpose=True) + 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") + + +""" +Display sinogram. +""" +fig, ax = plot.subplots(nrows=1, ncols=1, figsize=(15, 3)) +plot.imview(y, title="sinogram", fig=fig, ax=ax) +fig.show() + + +""" +Plot convergence statistics. +""" +fig, ax = plot.subplots(nrows=1, ncols=3, figsize=(12, 5)) +plot.plot( + np.vstack([hist[p].Objective for p in projectors.keys()]).T, + title="Objective function", + xlbl="Iteration", + ylbl="Functional value", + lgnd=projectors.keys(), + fig=fig, + ax=ax[0], +) +plot.plot( + np.vstack([hist[p].Prml_Rsdl for p in projectors.keys()]).T, + ptyp="semilogy", + title="Primal Residual", + xlbl="Iteration", + fig=fig, + ax=ax[1], +) +plot.plot( + np.vstack([hist[p].Dual_Rsdl for p in projectors.keys()]).T, + ptyp="semilogy", + title="Dual Residual", + xlbl="Iteration", + fig=fig, + ax=ax[2], +) +fig.show() + + +""" +Show the recovered images. +""" +fig, ax = plot.subplots(nrows=1, ncols=4, figsize=(15, 5)) +plot.imview(x_gt, title="Ground truth", fig=fig, ax=ax[0]) +for n, p in enumerate(projectors.keys()): + plot.imview( + x_rec[p], + title="%s SNR: %.2f (dB)" % (p, metric.snr(x_gt, x_rec[p])), + fig=fig, + ax=ax[n + 1], + ) +for ax in ax: + ax.get_images()[0].set_clim(-0.1, 1.1) +fig.show() + + +input("\nWaiting for input to close figures and exit") From f17ac18b1cfea3cc26639eea7119a26426d48b76 Mon Sep 17 00:00:00 2001 From: Brendt Wohlberg Date: Tue, 9 Jul 2024 16:49:59 -0600 Subject: [PATCH 29/33] Update example index --- docs/source/examples.rst | 6 +++--- examples/scripts/README.rst | 12 ++++++------ 2 files changed, 9 insertions(+), 9 deletions(-) diff --git a/docs/source/examples.rst b/docs/source/examples.rst index a5aa278e3..41a3df9e1 100644 --- a/docs/source/examples.rst +++ b/docs/source/examples.rst @@ -38,7 +38,7 @@ 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 ^^^^^^^^^^^^^ @@ -136,7 +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_cs_tv_admm + examples/ct_multi_tv_admm examples/ct_svmbir_tv_multi examples/deconv_circ_tv_admm examples/deconv_tv_admm @@ -208,7 +208,7 @@ ADMM examples/ct_tv_admm examples/ct_astra_3d_tv_admm examples/ct_astra_weighted_tv_admm - examples/ct_multi_cs_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 diff --git a/examples/scripts/README.rst b/examples/scripts/README.rst index 26a8f0a73..8c664c649 100644 --- a/examples/scripts/README.rst +++ b/examples/scripts/README.rst @@ -41,8 +41,8 @@ Computed Tomography CT Training and Reconstructions with UNet `ct_projector_comparison.py `_ X-ray Transform Comparison - `ct_multi_cs_tv_admm.py `_ - TV-Regularized Sparse-View CT Reconstruction (Multiple Projectors, Common Sinogram) + `ct_multi_tv_admm.py `_ + TV-Regularized Sparse-View CT Reconstruction (Multiple Projectors) Deconvolution ^^^^^^^^^^^^^ @@ -174,8 +174,8 @@ Total Variation 3D TV-Regularized Sparse-View CT Reconstruction (Proximal ADMM Solver) `ct_astra_weighted_tv_admm.py `_ TV-Regularized Low-Dose CT Reconstruction - `ct_multi_cs_tv_admm.py `_ - TV-Regularized Sparse-View CT Reconstruction (Multiple Projectors, Common Sinogram) + `ct_multi_tv_admm.py `_ + TV-Regularized Sparse-View CT Reconstruction (Multiple Projectors) `ct_svmbir_tv_multi.py `_ TV-Regularized CT Reconstruction (Multiple Algorithms) `deconv_circ_tv_admm.py `_ @@ -275,8 +275,8 @@ ADMM 3D TV-Regularized Sparse-View CT Reconstruction (ADMM Solver) `ct_astra_weighted_tv_admm.py `_ TV-Regularized Low-Dose CT Reconstruction - `ct_multi_cs_tv_admm.py `_ - TV-Regularized Sparse-View CT Reconstruction (Multiple Projectors, Common Sinogram) + `ct_multi_tv_admm.py `_ + TV-Regularized Sparse-View CT Reconstruction (Multiple Projectors) `ct_svmbir_tv_multi.py `_ TV-Regularized CT Reconstruction (Multiple Algorithms) `ct_svmbir_ppp_bm3d_admm_cg.py `_ From 275ef2c13fdc2873dc00edf943aef184664ea249 Mon Sep 17 00:00:00 2001 From: Brendt Wohlberg Date: Tue, 9 Jul 2024 16:59:50 -0600 Subject: [PATCH 30/33] Update submodule --- data | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/data b/data index 07fb707e3..866dfc544 160000 --- a/data +++ b/data @@ -1 +1 @@ -Subproject commit 07fb707e3c12f8a885778a4a06efafd93535c4b9 +Subproject commit 866dfc544deab6c08f141ea33b8fe4caaebf9d00 From 803b836c1879113f6273bd43fbde69a2235069bd Mon Sep 17 00:00:00 2001 From: Brendt Wohlberg Date: Fri, 19 Jul 2024 09:53:16 -0600 Subject: [PATCH 31/33] Update submodule --- data | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/data b/data index 866dfc544..cb97dd02c 160000 --- a/data +++ b/data @@ -1 +1 @@ -Subproject commit 866dfc544deab6c08f141ea33b8fe4caaebf9d00 +Subproject commit cb97dd02cfd96659804c74e7a5f9185440cec3ce From 2bf79a951109b3df7ab50748baf24404471041b2 Mon Sep 17 00:00:00 2001 From: Brendt Wohlberg Date: Fri, 19 Jul 2024 10:34:58 -0600 Subject: [PATCH 32/33] void doctest errors resulting from unimportable astra or svmbir --- .github/workflows/pytest_ubuntu.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.github/workflows/pytest_ubuntu.yml b/.github/workflows/pytest_ubuntu.yml index ffad0c04f..babfb1d9f 100644 --- a/.github/workflows/pytest_ubuntu.yml +++ b/.github/workflows/pytest_ubuntu.yml @@ -89,7 +89,7 @@ jobs: - name: Run doc tests if: matrix.group == 1 run: | - pytest --ignore-glob="*test_*.py" --doctest-modules scico + pytest --ignore-glob="*test_*.py" --ignore=scico/linop/xray --doctest-modules scico pytest --doctest-glob="*.rst" docs coverage: From 1c0cb426ac813f95c5d3d35a26d86e33cee4ba22 Mon Sep 17 00:00:00 2001 From: Brendt Wohlberg Date: Fri, 19 Jul 2024 10:46:23 -0600 Subject: [PATCH 33/33] Avoid doctest errors resulting from unimportable astra or svmbir --- docs/source/style.rst | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/source/style.rst b/docs/source/style.rst index f26bcb85a..80fe8fd27 100644 --- a/docs/source/style.rst +++ b/docs/source/style.rst @@ -439,7 +439,7 @@ classes, or method definitions. Comment explaining example 1. - >>> np.add(1, 2) + >>> int(np.add(1, 2)) 3 Comment explaining a new example.