PaperlayerView exact revision on arXiv

Source-linked AI summary

Convex optimization problem prototyping for image reconstruction in computed tomography with the Chambolle-Pock algorithm

Emil Y. Sidky, Jakob H. Jørgensen, Xiaochuan Pan

arXiv:1111.5632v2math.NAmath.OCphysics.med-ph

TL;DR

CT reconstruction formulations are difficult to prototype because imaging systems and tasks require varied optimization designs at very large scale. The paper applies the Chambolle–Pock primal-dual algorithm, deriving instances for many convex CT problems and demonstrating them in breast CT simulation. It reports solutions to difficult convex problems under challenging conditions in reasonable time without parameter tuning, while noting limitations of the prototyping approach.

  • Problem

    Large CT imaging models and varied optimization formulations make designing iterative reconstruction algorithms difficult to prototype.

  • Method

    The paper derives Chambolle–Pock algorithm instances for many convex CT reconstruction problems and applies them to a breast CT simulation.

  • Results

    The algorithm instances achieve solutions of difficult convex optimization problems under challenging conditions in reasonable time without parameter tuning.

  • Takeaways & Limitations

    The CP algorithm provides a tool for prototyping many CT optimization problems rather than advocating one formulation over another.

  • Takeaways & Limitations

    Optimization problem prototyping for CT image reconstruction has limitations, including in the breast CT application.

Abstract

from arXiv · show

The primal-dual optimization algorithm developed in Chambolle and Pock (CP), 2011 is applied to various convex optimization problems of interest in computed tomography (CT) image reconstruction. This algorithm allows for rapid prototyping of optimization problems for the purpose of designing iterative image reconstruction algorithms for CT. The primal-dual algorithm is briefly summarized in the article, and its potential for prototyping is demonstrated by explicitly deriving CP algorithm instances for many optimization problems relevant to CT. An example application modeling breast CT with low-intensity X-ray illumination is presented.

I. INTRODUCTION

CT reconstruction optimization problems vary with scanner properties, imaging tasks, representations, models, objectives, and constraints, making algorithm design difficult. The paper presents Chambolle–Pock as a convergent tool for rapidly prototyping many convex formulations.

  • CT optimization formulations depend on scanner type, imaging task, image representation, projection model, objective, and constraint design.
  • Changing an optimization formulation can require many weeks or months of algorithm development, motivating a mechanical prototyping tool.
  • The proposed tool generates a convergent algorithm from a specific CT reconstruction optimization problem, even if the resulting solver is not the most efficient.
  • The CP algorithm handles a broad class of convex problems, including non-smooth objectives and unconstrained or constrained formulations.
  • Large CT systems can involve millions of voxels and millions of X-ray transmission measurements, creating a major prototyping barrier.
  • The paper derives CP algorithm instances for CT problems and reports solutions to difficult convex problems under challenging conditions in reasonable time without parameter tuning.

A. Chambolle-Pock: Algorithm 1

The generic Chambolle–Pock method solves primal and dual convex problems through proximal mappings, and the paper turns CT formulations into concrete algorithm instances through five mechanical steps.

  • The CP algorithm simultaneously solves primal minimization and dual maximization problems for convex, possibly non-smooth functions linked by a linear transform.
  • In the paper’s usage, the CP setup has no free parameters, although convergence theory technically assumes L^2στ < 1.
  • The algorithm uses proximal mappings of F* and G; these mappings require minimizations that are simple enough to solve in closed form.
  • Prototyping maps the CT problem to the generic form, derives convex conjugates and proximal mappings, substitutes them into CP, and monitors the primal-dual gap.
  • Convex-function algebra, including non-smooth norms and indicator functions, permits many constrained and unconstrained CT formulations to be expressed in the generic form.

III. CHAMBOLLE-POCK ALGORITHM INSTANCES FOR CT

The paper formulates CT reconstruction problems within the Chambolle–Pock framework and derives corresponding primal-dual algorithm instances, convergence checks, and implementation requirements.

  • CT imaging model: CT reconstruction is modeled with a discrete linear system in which projection matrix A maps image coefficients u to line-integration data g.The model covers multiple expansion functions and both 2D fan-beam and 3D cone-beam configurations.
  • Primal-dual framework: The CP algorithm simultaneously solves primal and dual problems, with the primal-dual objective difference providing a convergence diagnostic in principle.At convergence, the primal and dual objectives become equal.
  • Convergence monitoring: Because the dual indicator imposes AT p = 0I, convergence monitoring uses a conditional primal-dual gap together with a separate check of AT p′ tending to 0I.The equality constraint is unlikely to be strictly satisfied during numerical computation, and the conditional gap should tend to zero.
  • Least-squares instance: The least-squares CP instance requires L = ∥A∥2 and an exact transpose AT of the projection matrix for guaranteed convergence.Obtaining an exact transpose is non-trivial when tomographic projection matrices are computed on the fly.
  • Underdetermined reconstruction: For underdetermined systems, CP returns a least-squares minimizer that can depend on the initial image, motivating objective augmentation to select a unique reconstruction.A quadratic term can favor an image close to 0I or a prior image.

1) Adding in non-negativity constraints:

The CP framework incorporates non-negativity and bound constraints by adding indicator functions to the primal objective, producing simple thresholding-based updates.

  • 1) Adding in non-negativity constraints:: Compared with conjugate gradients, adding constraints such as non-negativity is presented as considerably simpler in the CP formulation.For CP, the modification is achieved by introducing the appropriate indicator function.
  • 1) Adding in non-negativity constraints:: Adding an indicator for the set P of images with non-negative components extends the least-squares problem to enforce non-negativity.The constrained primal problem differs from the unconstrained case through the function G(x).
  • 1) Adding in non-negativity constraints:: The constrained formulation has dual feasibility checks requiring both AT p and u to be non-negative.These checks accompany the conditional primal-dual gap.
  • 1) Adding in non-negativity constraints:: The proximal mapping for the non-negativity indicator reduces to component-wise thresholding, so each iterate satisfies the non-negativity constraint automatically.The relevant update places the right-hand side inside the pos(·) operator.

B. Optimization problems based on the Total Variation (TV) semi-norm

The paper extends CP prototyping to TV-regularized CT objectives, including an explicitly derived ℓ2-TV instance and mechanisms for handling TV’s non-smoothness.

  • B. Optimization problems based on the Total Variation (TV) semi-norm: TV regularization is motivated by its use in CT reconstruction and its potential to provide high-quality images from sparse-view sampling.The TV semi-norm is associated with edge-preserving regularization.
  • B. Optimization problems based on the Total Variation (TV) semi-norm: The TV semi-norm is convex but non-linear, non-quadratic, and non-differentiable, preventing direct use of standard gradient-based optimization algorithms.The paper derives a CP instance for a TV-regularized least-squares objective.
  • B. Optimization problems based on the Total Variation (TV) semi-norm: The ℓ2-TV objective is expressed as F(Kx), where K maps an image to projection data and a spatial gradient, while KT combines AT and −div.The gradient ∇u is represented as a vector-valued image array in V.
  • B. Optimization problems based on the Total Variation (TV) semi-norm: The derived ℓ2-TV CP instance constrains the dual gradient variable through |q′| ∈ Box(λ) and enforces AT p′ − div q′ = 0I.The TV proximal operation rescales spatial vectors whose magnitudes exceed λ.
  • B. Optimization problems based on the Total Variation (TV) semi-norm: The CP framework can also incorporate positivity into the TV formulation by adding the indicator δP(u) to the primal objective.This introduces the pos(·) operator in the corresponding update.

1) Alternate data divergences:

The paper derives CP instances for alternate CT data divergences, combining TV with KL or ℓ1 data terms and preserving direct handling of non-smooth objectives.

  • 1) Alternate data divergences:: KL divergence is motivated for CT settings with significant noise modeled by a multivariate Poisson distribution and can be combined with TV for sparse-view sampling.The resulting KL-TV objective includes a non-negativity indicator.
  • 1) Alternate data divergences:: The paper derives an ℓ1-TV CP instance in which the non-smooth objective can be handled without smoothing parameters for the ℓ1 and TV terms.The corresponding pseudocode differs from the preceding instances at the data-term update.
  • 1) Alternate data divergences:: The KL formulation requires a positive linear system A and non-negative image pixels u and data g.An indicator function restricts the domain so the objective remains defined over the whole space.
  • 1) Alternate data divergences:: The KL-TV CP instance differs from the preceding ℓ2-TV instance only in the update for the data term.Intermediate image estimates may be negative even though the converged solution is non-negative.
  • 1) Alternate data divergences:: The ℓ1 data-error term is presented as a robust fit that weakly penalizes outliers and may reduce sensitivity to non-uniform projection-domain errors.The paper relates these errors to discretization and metal objects, which can otherwise contribute to streak artifacts under quadratic fitting.

2) Constrained, TV-minimization:

This section reformulates least-squares data fidelity with TV regularization as a constrained TV-minimization problem, using a data-error bound ε instead of the less physically interpretable λ. It derives the corresponding Chambolle–Pock algorithm instance and its key proximal mappings.

  • Problem formulation: Constrained TV-minimization replaces the TV-penalized formulation with an ℓ2 data-error constraint.The reformulation uses δBall(ε)(Au − g), which is zero when ∥Au − g∥2 ≤ ε.
  • Problem formulation: For every positive ε, a corresponding λ yields the same solution as the unconstrained formulation.This establishes equivalence between the constrained and unconstrained problems in the sense of matching solutions.
  • CP derivation: The constrained formulation produces both primal and dual problems through the CP derivation.The section defines the primal function and its conjugate before presenting the dual problem.
  • CP derivation: The algorithm instance requires proximal mappings for the conjugate data term and the primal regularizer.These mappings are the resolvent operators used in the CP framework.
  • Algorithm instance: The resulting pseudocode solves the ℓ2-constrained, TV-minimization problem, with ε interpreted as the data-error bound.The listed instance initializes u0, p0, and q0 to zero values.

IV. DEMONSTRATION OF CP ALGORITHM INSTANCES FOR TOMOGRAPHIC IMAGE RECONSTRUCTION

The paper demonstrates that CP algorithm instances can be derived for many convex CT reconstruction problems and used to prototype alternatives without selecting one formulation as universally preferable. A simulated sparse-view breast CT study uses Poisson noise to compare KL-TV and least-squares-TV reconstructions, finding perceptibly lower ROI noise for KL-TV while noting that firmer conclusions require multiple noise realizations.

  • Demonstration scope: CP instances are guaranteed to converge across practically infinite combinations of system matrix, data, and optimization parameters.The paper presents this convergence as applying to the many optimization problems covered by the derived instances.
  • Demonstration scope: The results section demonstrates CP's utility for optimization-problem prototyping rather than advocating one formulation over another.Example reconstructions are used to illustrate the framework's practical utility.
  • Sparse-view experiment: The sparse-view experiment models breast CT with 60 projections over a full 360° range and Poisson-distributed measurement noise.The simulation uses a digital 256 × 256 breast phantom and selects Poisson noise to examine matching the data-error term to the noise model.
  • Sparse-view experiment: The experiment compares KL-TV, motivated by maximum likelihood under the Poisson model, with least-squares-TV reconstruction.Both objectives are readily prototyped within CP without smoothing the TV semi-norm, setting algorithm parameters, or proving convergence.
  • Sparse-view results: Perceptible noise reduction appears in the ROIs of KL-TV images, while the two optimization problems show other visual differences.The images address fibro-glandular morphology and micro-calcification detection or classification tasks.
  • Sparse-view results: A firm conclusion about the reconstruction comparison awaits a more complete study with multiple noise realizations.The authors identify multiple realizations of the random data as a necessary next step for the CT simulation conclusion.
  • Broader applicability: The CP framework supports prototyping other convex formulations, including sparsity-based penalties, TV constraints, alternative data fidelities, and physical image constraints.The paper gives least-squares, ℓ1, and TV combinations for sparse-view phase-contrast imaging as an example.

C. Iteration dependence of the CP algorithm

The CP algorithm is used to prototype convex CT reconstruction formulations, with convergence checked through the conditional primal-dual gap and problem-specific indicator conditions. In a KL-TV example, smaller λ values require more iterations, while the broader workflow separates formulation design from later efficiency optimization.

  • Iteration dependence: For KL-TV, condition 4 is automatically enforced, condition 3 is satisfied early, and condition 1 is illustrated as converging in the plotted diagnostic.Condition 2 has a similar dependence to condition 1.
  • Iteration dependence: The KL-TV experiment evaluates conditional primal-dual-gap convergence for different λ values and checks agreement with indicator-function conditions.The plotted checks include div q − AT p = 0I, Au ≥ 0D, p ≤ 1D, and |q| < λ.
  • Discussion: The CP framework derives convergent algorithm instances for many convex CT image-reconstruction optimization problems.It is intended for rapid prototyping rather than necessarily being the most efficient solver.
  • Discussion: The CP algorithm is most useful during early formulation design, after which attention shifts from convergence and optimization parameters to efficiency.This separates optimization parameters from algorithm parameters during prototyping.
  • Limitations: Clinical iterative reconstruction may use only one to ten iterations, potentially too few to claim an accurate solution to the designed optimization problem.The paper presents this as a limitation of applying prototyping insights directly to severely truncated iteration counts.

APPENDIX A

The appendix supplies computational details for CP implementations, including operator-norm estimation, convex-conjugate derivations, and proximal mappings. These derivations use case analysis, geometric reasoning, and componentwise constraints to obtain the required updates.

  • Operator norms: The CP parameter L uses the largest singular value of K, estimated with a standard power method.For a discrete X-ray transform, the power method converges to numerical precision in twenty iterations or less in the authors’ experience.
  • Operator norms: TV-penalized minimization requires the norm of the combined transform (A, ∇), whose power-method implementation uses AT A x_n − div ∇ x_n.The identity −div = ∇T is used in this computation.
  • Convex conjugates: The appendix derives convex conjugates of one-dimensional indicator functions by analyzing separate regions and generalizing the resulting cases to higher dimensions.Figures 4–6 visualize the objectives for two- or three-case maximizations.
  • Convex conjugates: For δP, the conjugate is 0 when x ≤ 0 and unbounded when x > 0; the x ≤ 0 maximum occurs at x′ = 0.The two cases are illustrated in Fig. 4.
  • Convex conjugates: For δBox(1), the one-dimensional conjugate has maxima at x′ = −1 or x′ = 1 outside the central case, with values determined by −x or x.The appendix distinguishes x ≤ 0 and x > 0 in one dimension and notes that multidimensional box and ball indicators differ.
  • Proximal mappings: The proximal mappings are obtained by enforcing magnitude constraints componentwise, selecting the closest admissible vector and scaling when |z_i| exceeds λ.The scaling is λz_i/|z_i|, and σ does not enter this calculation.
  • Proximal mappings: The KL data-term proximal mapping is derived by differentiating a smooth objective and selecting the negative quadratic root to preserve 1D − p′ ≥ 0.The non-negativity of the data is used in checking the admissible root.
  • Proximal mappings: For the constrained TV problem, completing the square and geometric reasoning place the minimizer on the line segment between 0D and p − σg.Analyzing the resulting one-dimensional minimization yields Eq. (59).

THE FINITE DIFFERENCING FORM OF THE IMAGE GRADIENT AND DIVERGENCE

The paper specifies discrete gradient and divergence operators for CT images, emphasizing their transpose relationship, then describes a preconditioned CP variant for KL-TV problems that can improve efficiency for small λ.

  • THE FINITE DIFFERENCING FORM OF THE IMAGE GRADIENT AND DIVERGENCE: The discrete image domain uses an M × M image x ∈ I, with x_i,j denoting its (i, j)th pixel.
  • THE FINITE DIFFERENCING FORM OF THE IMAGE GRADIENT AND DIVERGENCE: The two-dimensional gradient ∇ is defined through horizontal and vertical differencing images Δ_sx and Δ_tx.
  • THE FINITE DIFFERENCING FORM OF THE IMAGE GRADIENT AND DIVERGENCE: Pixels referenced outside the image border are set to zero in the finite-difference definitions.
  • THE FINITE DIFFERENCING FORM OF THE IMAGE GRADIENT AND DIVERGENCE: The discrete operators ∇ and −div must be transposes of each other, even though the particular finite-difference form is not essential.
  • PRECONDITIONED CHAMBOLLE-POCK ALGORITHM DEMONSTRATED ON THE KL-TV OPTIMIZATION PROBLEM: The preconditioned CP algorithm replaces scalar σ and τ with vector quantities computed directly from the system matrix K.
  • PRECONDITIONED CHAMBOLLE-POCK ALGORITHM DEMONSTRATED ON THE KL-TV OPTIMIZATION PROBLEM: A diagonal preconditioner uses weights involving Σ ∈ Y, T ∈ X, and the elementwise absolute-value matrix |K|.
  • PRECONDITIONED CHAMBOLLE-POCK ALGORITHM DEMONSTRATED ON THE KL-TV OPTIMIZATION PROBLEM: Preconditioning modifies the proximal mapping, whose quadratic term is no longer spherically symmetric and can therefore make derivation more involved.
  • PRECONDITIONED CHAMBOLLE-POCK ALGORITHM DEMONSTRATED ON THE KL-TV OPTIMIZATION PROBLEM: For 2-TV and KL-TV with small λ, absorbing λ into the TV term lets Σ address mismatch between TV and data-agreement terms, while the resulting algorithm remains an alternative prototyping method.
Loading 1111.5632v2…