Investigation of the preconditioner-parameter in the preconditioned Chambolle-Pock algorithm applied to optimization-based image reconstruction

Abstract

The optimization-based image reconstruction methods have been thoroughly investigated in the field of medical imaging. The Chambolle-Pock (CP) algorithm may be employed to solve these convex optimization image reconstruction programs. The preconditioned CP (PCP) algorithm has been shown to have much higher convergence rate than the ordinary CP (OCP) algorithm. This algorithm utilizes a preconditioner-parameter to tune the implementation of the algorithm to the specific application, which ranges from 0 and 2, but is often set to 1. In this work, we investigated the impact of the preconditioner-parameter on the convergence rate of the PCP algorithm when it is applied to the TV constrained, data-divergence minimization (TVDM) optimization based image reconstruction. We performed the investigations in the context of 2D computed tomography (CT) and 3D electron paramagnetic resonance imaging (EPRI). For 2D CT, we used the Shepp-Logan and two FORBILD phantoms. For 3D EPRI, we used a simulated 6-spheres phantom and a physical phantom. Study results showed that the optimal preconditioner-parameter depends on the specific imaging conditions. Simply setting the parameter equal to 1 cannot guarantee a fast convergence rate. Thus, this study suggests that one should adaptively tune the preconditioner-parameter to obtain the optimal convergence rate of the PCP algorithm.

Keywords

Optimization-based image reconstruction total variation minimization Chambolle-Pock algorithm preconditioner convergence rate

1 Introduction

Optimization-based image reconstruction [1 –5] remains the focus of current research in the field of medical imaging. Evidence appears to show that appropriately designed optimization-based algorithms (OA) may outperform the analytic algorithms if the data is insufficient or inconsistent [6 –8]. The reason may be that OA can effectively model the sparsity of the reconstructed object (image or signal) or that of its transform in the optimization programs. If the sparsity is employed by minimizing the sparse transform, the class of OA is based on compressed sensing (CS) [9]. If the sparsity is employed by minimizing the rank of the reconstructed object, it is based on partial separability (PS) [10].

It appears that CS-based OA is more popular than PS-based OA because a sparse transform of the reconstructed object may always be obtained. An example of such a CS-based OA is to employ the discrete gradient transform (DGT) [11]. Total variation (TV) minimization is a popular and effective OA algorithm for constructing images accurately from sparse-view projections [1]. There are 3 frequently used algorithms to solve TV minimization: adaptive-steepest-decent-projection-onto-convex-sets (ASD-POCS) [2], alternating direction method of multipliers (ADMM) [12, 13] and Chambolle-Pock (CP) [14 –16] algorithms. The CP algorithm is a first-order primal-dual algorithm framework for solving convex problems [16]. It guarantees convergence and may solve both smooth and non-smooth convex problems [17], e.g. the ℓ₁ norm minimization problem. Thus, the CP algorithm is a powerful tool for solving the TV minimization program and has been applied more and more.

The work will focus on the CP algorithm applied to the TV constrained, data-divergence minimization (TVDM) [17 –20]. There are two types of CP algorithms, the ordinary CP (OCP) and preconditioned CP (PCP) [15]. Research has shown that PCP has a much faster convergence rate than OCP. In [21], the authors proposed the PCP algorithm and applied it to linear programming, TV-ℓ₁ image denoising, graph cuts and minimal partitions. Their work demonstrated that PCP leads to a significantly faster convergence rate. In [15], the authors applied PCP to unconstrained TV minimization and achieved convergence rate enhancement. In [22], the authors also applied PCP to TV minimization via simulated data and observed a similar acceleration effect. In all of these instances the preconditioner-parameter,, is set to be 1 (see Section 3 of [21], Equation E.1 and 2 of [15] and Algorithm 2 of [22]). However, the general form of the preconditioner allows for tuning of the parameter α within a range of [0,2]. It is not necessarily true that α = 1 is the optimal selection for achieving the fastest convergence. Thus, the impact of the preconditioner-parameter on the convergence rate of the PCP algorithm should be deeply investigated.

In the work, we use the PCP algorithm tailored to TVDM to perform such an investigation. To observe the behavior of different α values under different imaging conditions, we design TVDM-PCP algorithms for 2D computed tomography (CT) and 3D electron paramagnetic resonance imaging (EPRI). For CT, we use two FORBILD phantoms (Study 1 and 3) and the Shepp-Logan phantom (Study 2); for EPRI, we use a simulated 6-spheres phantom (Study 4) and a physical phantom (Study 5). For each investigation, we set up 3 groups of balance parameters and use 5 different α values to perform the reconstructions for each group of balance parameters. Thus, we evaluate the selection of α on the convergence rate for 15 separate image reconstruction cases.

We organize the paper as follows. In Section 2, we show the image reconstruction chain, from imaging model, optimization program to PCP algorithm. In Section 3, we show the reconstruction results and investigate the dependence on α. In Section 4, we provide a theoretical analysis and present a brief discussion and our conclusions.

2 Methods

In this section, we describe the entire image reconstruction chain: imaging system modeling, optimization program, PCP algorithm and parameter selection. Because we aim to study the convergence property of the PCP algorithm applied to TVDM in both CT and EPRI, the reconstruction-chain description will be more generalized rather than focused on any specific imaging modality.

2.1 Imaging system modelling

Iterative algorithms are, in fact, based on discrete-to-discrete (D2D) models. D2D means that the projection and the image are both discrete. For 2D parallel CT, the projection set is the 2D radon transform of the 2D object; for 3D EPRI, the projection set is the 3D radon transform of the 3D object. To describe the whole reconstruction chain generally, we will consider the generalized radon transform rather than its specific form applied to either 2D CT or 3D EPRI. The Radon transform is a forward process to convert the object into projections, whereas the inverse radon transform is a process to convert the projections into the object.

For the 2D radon transform, a measurement on a projection signal is the line integral of the image (object) function on the line specified by the measurement. For the 3D radon transform, a measurement on the projection signal is the plane integral of the object function on the plane specified by the measurement. The D2D model of the imaging system based on Radon transform may be modeled as $g = Au$ (1) where, g of size M denotes the discrete projection data, u of size N the discrete image, and A the system matrix of size M × N, corresponding to the Radon transform, meaning that A_m,n is the contribution of the nth pixel (voxel) to the mth measurement.

For the 2D radon transform, A_m,n is the intersecting length of the nth pixel with the mth line. For the 3D radon transform, A_m,n is the interesting area of the nth voxel with the mth plane. A pixel is regarded as a square and a voxel as a cube; hence the key for calculating A_m,n is calculating the intersecting length of a line and a square or the intersecting area of a plane and a cube. In fact, the method used to calculate these values for A_m,n is the so-called projection method. The direct length calculation forms the Siddon ray-driven method [23] and area calculation may borrow and extend the method. By the use of interpolation in the projection generation process in the 2D case, the Joseph ray-driven [24] and distance-driven projection [25] methods were proposed. The projection method for the 3D case may also be designed by extending the two methods. However, these methods are mathematically complex, especially in 3D case. In this work, we use the pixel-driven projection and backprojection methods.

Usually, the linear system is large-scale, often underdetermined and ill-posed, so it is not practical to solve the system by simply performing an inversion to the system matrix. It is possible that some useful and desired solutions of the linear system may be obtained by further modelling the D2D data model into an optimization program.

2.2 Optimization program

The solution to the imaging model discussed above can be obtained via an optimization program. Here, we consider the TVDM program shown as $u^{*} = argmin {∥ g - Au ∥}_{2}^{2} s . t . {∥ u ∥}_{TV} ⩽ t_{1}$ (2) where u^* denotes the solution, i.e. the reconstructed image, ${∥ \cdot ∥}_{2}^{2}$ the square of the ℓ₂-norm of the difference between measured data g and modelled data Au, and image-TV ∥u ∥ _TV the ℓ₁-norm of the gradient magnitude image |Du|_mag is given by ${∥ u ∥}_{TV} = {∥ {| Du |}_{mag} ∥}_{1}$ (3) where D is a matrix of size 3N × N (for the 3D case only, the 2D case may be obtained by omitting one dimension), expressed as $D = (\begin{matrix} D_{1} \\ D_{2} \\ D_{3} \end{matrix})$ (4) where, D₁, D₂, and D₃, matrices of size N × N, denote the discrete gradient transforms along the directions of the x, y and z, respectively.

Let u_x,y,z be the voxel at index (x, y, z) x ∈ [1, N_x], y ∈ [1, N_y], z ∈ [1, N_z]. Transform D₁ can be expressed as $(D_{1} u)_{x, y, z} = {\begin{matrix} u_{x, y, z} - u_{x - 1, y, z} x \in [2, N] \\ 0 x = 1 \end{matrix}$ (5) and similar expressions can readily be obtained for D₂ and D₃. Finally, we obtain the gradient magnitude image in terms of image vector u. ${({| Du |}_{mag})}_{i} = \sqrt{{(Du)}_{i}^{2} + {(Du)}_{i + N}^{2} + {(Du)}_{i + 2 N}^{2}} i \in [1, N]$ (6)

To improve the convergence rate, two algorithm parameters λ and ν are introduced to control and adjust the relative amplitudes of the two convex functions, data-ℓ₂ norm and image-TV norm. Thus we get the new optimization program $u^{*} = argmin \frac{λ}{2} {∥ g - Au ∥}_{2}^{2} s . t . ν {∥ u ∥}_{TV} ⩽ ν t_{1},$ (7)

The two parameters do not impact the solution of the optimization program but may impact the convergence rate of the algorithm designed to solve the optimization program. Since the two parameters appear to aim at achieving the balance point between the two convex functions, they will be referred to as balance parameters.

2.3 The general PCP algorithm

The TVDM program may be solved by the Chambolle-Pock algorithm. The algorithm tailored to solve Equation (7) has been similarly designed in [17 –20]. However, in fact, these algorithms are all using a special form of the general PCP algorithm. In [21], the authors show the general PCP algorithm frame in Lemma 2. Applying the general form to the CP algorithm instance tailored to TVDM, we get the pseudo-code of the general TVDM-PCP algorithm, shown in Table 1.

Table 1
Pseudo code of the general PCP algorithm for solving Equation (7)

INPUT: g, A, and t₁

1: input λ, b

2: $v = {bv}_{A} = b \frac{| | A | |_{sv}}{| | D | |_{sv}}$

3: $\sum_{1} = \frac{1_{D}}{| A |^{α} 1_{I}}; \sum_{2} = \frac{1_{V}}{| vD |^{α} 1_{I}}; T = \frac{1_{I}}{| A^{T} |^{2 - α} 1_{D} + | (vD)^{T} |^{2 - α} 1_{ν}}; θ = 1; n = 0$

4: $u_{0} = 0; {\bar{u}}_{0} = 0; p_{0} = 0; q_{0} = 0$

5: repeat

6:        $p_{n + 1} = (p_{n} + \sum_{1} (A {\bar{u}}_{n} - g)) / (1 + \sum_{1} / λ)$

7: 7.1 $a = q_{n} + \sum_{2} vD {\bar{u}}_{n}$

7: 7.2 s = ProjectOnto ℓ ₁ Ball_{vt
₁} (|a|_mag/Σ₂)

7: 7.3 q_n+1 = a (1_v - Σ₂s/|a|_mag)

8:       u_n+1 = u_n - T (A^T p_n+1 + vD^T q_n+1)

9:        ${\bar{u}}_{n + 1} = u_{n + 1} + θ (u_{n + 1} - u_{n})$

10:     n = n + 1

11: until the practical convergence conditions are satisfied

OUTPUT: image u_conv = u_N

The 3D case is used to illustrate the code, however the 2D case is simpler than the 3D case and easily derived from the 3D case. In Table 1, ∥ · ∥ _sv denotes the largest singular value of a matrix, which can be calculated by use of Algorithm 8 in [15]. Let I = R^N and D = R^M denote an image space and a data space, respectively. Therefore, we have gradient-image space, V = I³ = R^3N. $u, \bar{u}, s \in I; p \in D; and q, a \in V$ . In line 7.3, 1_V is a unit vector in space V. In line 7.2 the projection operation, ProjectOnto ℓ₁ Ball, can be computed with the pseudo code shown in [18].

According to the specific form of |D| and |D^T|, Σ₂ and T may be simplified. $Σ_{2} = \frac{1_{V}}{{| ν D |}^{α} 1_{I}} = \frac{1_{V}}{ν^{α} {| D |}^{α} 1_{I}} = \frac{1_{V}}{ν^{α} | D | 1_{I}} = \frac{1_{V}}{ν^{α} \times 21_{V}} = \frac{1}{2 ν^{α}}$ (8) for | (νD) ^T|^2-α = ν^2-α|D^T|^2-α = ν^2-α|D^T| $T = \frac{1_{I}}{| A^{T} |^{2 - α} 1_{D} + | (ν D)^{T} |^{2 - α} 1_{ν}} = \frac{1_{I}}{| A^{T} |^{2 - α} 1_{D} + ν^{2 - α} | D^{T} | 1_{ν}}$ and since |D^T|1_V = 61_I $T = \frac{1_{I}}{| A^{T} |^{2 - α} 1_{D} + v^{2 - α} | D^{T} | 1_{v}} = \frac{1_{I}}{| A^{T} |^{2 - α} 1_{D} + v^{2 - α} \times 61_{I}} = \frac{1_{I}}{| A^{T} |^{2 - α} 1_{D} + 6 v^{2 - α}}$ (9)

2.4 Parameter selection

In the PCP algorithm for solving TVDM, there are mainly two groups of parameters: program parameters and algorithm parameters. Program parameters define the final solution to the optimization program, including TV constraint value t₁, the projection-method selection, the definition of TV norm etc. Algorithm parameters, however, may not impact the final solution but impact the convergence rate, including λ, ν, Σ₁, Σ₂, T and θ. λ and ν are the two balance parameters, which have been shown to impact the convergence rate [19]. Σ₁, Σ₂ and T are the inherent algorithm parameters in the PCP algorithm. They are termed preconditioner in [21]. These three preconditioners may be tuned by choosing the value of α. In this work, the main task is to evaluate how, if at all, selection of α impacts the convergence rate.

3 Results

To extensively explore the impact of the preconditioner-parameter, α, on convergence rate, we perform 5 studies: (1) 2D complete FORBILD phantom simulation CT-study; (2) 2D Shepp-Logan phantom simulation CT-study; (3) 2D simple FORBILD phantom simulation CT-study; (4) 6-spheres phantom simulation EPRI-study and (5) Experimental 3D data EPRI-study.

3.1 Evaluations of the CT reconstructions of 2D complete FORBILD phantom with two ears

The size of the complete FORBILD phantom is 500 × 500, which is scanned using parallel beam X rays. The size of the projection-data set is also 500 × 500, meaning that the number of views is 500 and the number of measurements on each projection is also 500. The pixel size and detector-cell size are both normalized to be 1. The projection views are distributed uniformly in the range [0, π].

We set up λ and ν into three combinations: (1) λ = 1, ν = 0.1ν_A; (2) λ = 1, ν = ν_A; (3) λ = 1, ν = 10ν_A. We set up α to be 0, 0.5, 1, 1.5 and 2 and evaluate their impact on convergence rate. We let each reconstruction stop at iteration 2000 and compare the root-mean-square-error (RMSE) of the reconstructed objects to perform the evaluation of the convergence rate. Lower RMSE indicates faster convergence, whereas higher RMSE indicate slower convergence. Figure 1 shows the reconstructed images; Fig. 2 shows the iteration trend of the RMSE for each reconstruction and Table 2 shows the RMSE of each reconstruction at iteration 2000.

Fig.1

The reconstructed complete FORBILD-images via the TVDM algorithm. The text above the images shows the values of the two balance parameters. The text on the left of the images shows the values of α. The display window is [0, 1]. We use image(i, j) to indicate the image at the ith row and the jth column (the same below).

Fig.2

Comparison of iteration trend of the RMSE of the reconstructed objects. (a) is the RMSE comparison under the condition that λ = 1, ν = 0.1ν_A; (b) under the condition that λ = 1, ν = ν_A; and (c) under the condition that λ = 1, ν = 10ν_A.

Table 2

RMSE of the reconstructed complete FORBILD-images via different α

RMSE (×10^-4)	α = 0	α = 0.5	α = 1	α = 1.5	α = 2
λ = 1, ν = 0.1v_A	5.25	3.14	3.28	4.55	12.00
λ = 1, ν = v_A	310.8	8.82	3.26	3.35	12.00
λ = 1, ν = 10v_A	1292.9	309.7	3.49	2.51	12.00

From Fig. 1, we may see that α can have an appreciable impact on the convergence rate, observing that image (1, 3) and image (2, 3) still have large RMSE at iteration 2000 compared to the truth. From Table 2 and Fig. 2, we see that α = 1 must not always achieve the fastest convergence. For the case of λ = 1, ν = 0.1ν_A, α = 0.5 achieves the fastest convergence as its RMSE is the lowest; for the other two cases, α = 1 and 1.5 achieves the fastest convergence, respectively. For the curve comparison and value comparison of RMSE are consistent, we will just show the RMSE-value comparisons in the following studies.

3.2 Evaluations of the CT reconstructions of 2D Shepp-Logan phantom

We continue to perform investigation on a small size phantom relative to the last study for maintaining the extensiveness of the whole investigation. The size of the digital Shepp-Logan phantom is 64 × 64, which is scanned using parallel beam X rays. The size of the projection-data set is also 64 × 64, meaning that the number of views is 64 and the number of measurements on each projection is also 64. The pixel size and detector-cell size are both normalized to be 1. The projection views are distributed uniformly in the range [0, π].

The settings of λ, ν and α are the same with those in study 1 (Section 3.1). Figure 3 shows the reconstructed images and Table 3 shows the RMSE of each reconstruction at iteration 2000.

Fig.3

The reconstructed Shepp-Logan-images via the TVDM algorithm. The text above the images shows the values of the two balance parameters (λ, ν). The text on the left of the images shows the values of the preconditioner-parameter, (α). The display window is [0, 1].

Table 3

RMSE of the reconstructed Shepp-Logan-images via different α

RMSE (×10^-6)	α = 0	α = 0.5	α = 1	α = 1.5	α = 2
λ = 1, ν = 0.1v_A	5.80	6.30	7.46	7.24	7.43
λ = 1, ν = v_A	10000	1100	7.16	4.45	5.65
λ = 1, ν = 10v_A	113100	43800	4000	8.89	10.03

From Fig. 3, we may see that α can have an appreciable impact on the convergence rate, observing that image (1, 3) and image (2, 3) still have large RMSE at iteration 2000 compared to the truth. From Table 3, we see that α = 1 does not achieve the fastest convergence. For the case of λ = 1, ν = 0.1ν_A, α = 0 achieves the fastest convergence as its RMSE is the lowest; for the other two cases, α = 1.5 achieves the fastest convergence.

3.3 Evaluations of the CT reconstructions of a 2D simple FORBILD phantom

For maintaining the extensiveness of the investigation, we perform reconstructions on a medium size FORBILD phantom. The size of the FORBILD phantom is 100 × 100, which is scanned using parallel beam X rays. The size of the projection-data set is also 100 × 100, meaning that the number of views is 100 and the number of measurements on each projection is also 100. The pixel size and detector-cell size are both normalized to be 1. The projection views are distributed uniformly in the range [0, π].

The settings of λ, ν and α are the same with those in study 1 (Section 3.1). Figure 4 shows the reconstructed images and Table 4 shows the RMSE of each reconstruction at iteration 2000.

Fig.4

The reconstructed FORBILD-images via the TVDM algorithm. The text above the images shows the values of the two balance parameters. The text on the left of the images shows the values of α. The display window is [0, 1].

Table 4

RMSE of the reconstructed FORBILD-images via different α

RMSE (×10^-5)	α = 0	α = 0.5	α = 1	α = 1.5	α = 2
λ = 1, ν = 0.1v_A	1.36	1.05	1.45	2.36	3.66
λ = 1, ν = v_A	440	48.9	1.13	1.46	4.64
λ = 1, ν = 10v_A	9580	2320	150	1.12	4.10

From Fig. 4, we may see once more that α can have an appreciable impact on convergence rate, observing that image (1,3) and image (2,3) still have large errors compared to the truth at iteration 2000. From Table 4, we see that α = 1 must not always achieve the fastest convergence. For case of λ = 1, ν = 0.1ν_A, α = 0.5 achieves the fastest convergence for its RMSE is the lowest; for the other two cases, α = 1 and 1.5 achieves the fastest convergence, respectively.

3.4 Evaluations of the EPRI reconstructions of a simulated 3D 6-spheres phantom

Further, we evaluate the impact of the preconditioner parameter on convergence rate via a 3D 6-spheres phantom. The phantom consists of 6 spheres with each one having a different uniform value. Five small spheres are embedded within a larger sphere.

The size of the digital phantom is 32 × 32 × 32. The size of the projection-data set is 32 × 828. There are 828 projections which are distributed in the full range using an equal-solid-angle pattern [28]. The number of measurements on each projection is 32. The size of the voxels and the visual detector cells are both normalized to be 1.

The settings of λ, ν and α are the same with those in study 1 (Section 3.1). Figure 5 shows the reconstructed images and Table 5 shows the RMSE of each reconstruction at iteration 2000. From Fig. 5, again we may see that α can have an impact on convergence rate, observing that image (1,:), image (2,2:3) and image (3,3) still have big error from the truth at iteration 2000. From Table 5, we see that α = 1 does not necessarily achieve the fastest convergence. For the case of λ = 1, ν = 0.1ν_A, α = 0.5 achieves the fastest convergence as its RMSE is the lowest; for the other two cases, α = 1 and 1.5 achieves the fastest convergence, respectively.

Fig.5

The slice-images of the reconstructed 6-spheres-objects via the TVDM algorithm. The text above the images shows the values of the two balance parameters. The text on the left of the images shows the values of α. The display window is [0, 1].

Table 5

RMSE of the reconstructed 6-spheres-phantoms via different α

RMSE (×10^-4)	α = 0	α = 0.5	α = 1	α = 1.5	α = 2
λ = 1, ν = 0.1v_A	24	4.25	6.88	28	67
λ = 1, ν = v_A	411	72	3.02	8.92	67
λ = 1, ν = 10v_A	1325	506	38	3.52	67

3.5 Evaluation of the EPRI reconstructions of a physical phantom

Finally, we use real experimentally acquired data to evaluate the role that the selection of α value plays in the convergence of the reconstructed images. We made a physical phantom consisting of 3 small bottles and 3 tubes, each filled with uniform solutions containing different concentrations of triarylmethyl, the EPRI spin probe imaging agent.

A pulse 250 MHz imager [26] controlled by SpecMan4EPR v2.1 (Femi instruments, Chicago, IL) [27] is used. The gradient strength is 0.75G/cm and the scanning time for the full projection set, i.e. 828 projections, is 5.5 minutes. These setting are standard and have been found to achieve a good balance between spatial resolution, which is impacted by gradient strength, and signal to noise radio (SNR), which is impacted by scanning time.

The object is of size 80 × 80 × 80; the field of view (FOV) is 7cm × 7cm × 7cm so the voxel size is 0.0875cm × 0.0875cm × 0.0875cm; the number of spatial projections is 828; the number of measurement points on a spatial projection is 80; the length of the spatial projection is 7 cm so the sampling interval of the spatial projection is also 0.0875cm; the sampling pattern is the equal solid angle pattern and the sampling order is MSPS [28].

The TV constraint value, t₁ = 0.275t₀ is used. Here, t₀ is the TV value of the reconstructed object when the analytic filtered backprojection (FBP) algorithm is used [29]. We use three combinations of λ and ν: (1) λ = 1, ν = 0.1ν_A; (2) λ = 1, ν = ν_A; (3) λ = 1, ν = 10ν_A. We set up α to be 0, 0.5, 1, 1.5 and 2 and evaluate their impact on convergence rate. We use a different method here with discussed above to evaluate the convergence rate of the reconstructions when varying α. If the normalized RMSE of the reconstructed object is less than or equal to 3, then we stop the iteration and record the iteration index. Smaller iteration index means faster convergence. Figure 6 shows the reconstructed images and the photo of the phantom and Table 6 shows the iteration index achieving the stopping condition.

Table 6
Iteration index achieving the stopping condition for different α in experimental data study

Iteration index α = 0 α = 0.5 α = 1 α = 1.5 α = 2

λ = 1, ν = 0.1v_A 110 93 84 82 87

λ = 1, ν = v_A 286 230 199 179 168

λ = 1, ν = 10v_A 909 738 637 568 518

Iteration index	α = 0	α = 0.5	α = 1	α = 1.5	α = 2
λ = 1, ν = 0.1v_A	110	93	84	82	87
λ = 1, ν = v_A	286	230	199	179	168
λ = 1, ν = 10v_A	909	738	637	568	518

From Fig. 6, we may see that the unpaired-electron spin-density (UESD) image can be clearly reconstructed. Red color indicates a high UESD value, whereas yellow and blue colors indicate a low UESD value. The image is consistent with the fact that phantom is comprised of one high density bottle, two lower density bottles, one high density tube and two lower density tubes.

Fig.6

(a) is the physical phantom consisting of 3 small bottles and 3 tubes. (b) is the slice-image at y = 44 of the reconstructed object with the setting that λ = 1, v = 0.1v_A and α = 2. The display window is [0, 0.0036]. The other images reconstructed using other settings are visually similar with (b) for they are all obtained at the stop condition that NRMSE is less than or equal to 3.

From Table 6, we see that α = 1 cannot achieve the fastest convergence. For the case of λ = 1, ν = 0.1ν_A, the fastest convergence is achieved with α = 1.5 because the iteration stopping condition is reached at iteration 82. For the other two cases, α = 2 achieves the fastest convergence. In the experimental real-data study, the truth object is not known, so we use the FBP-reconstructed object as a reference object to calculate the NRMSE between the TVDM-reconstructed object and the reference object. According to our experience, when the NRMSE goes down to 3, the image will be the desired image of practical utility. In fact, if we let the iteration number be large enough, we may find that the final image is visually almost the same with the image achieved at the stop condition that NRMSE is less than or equal to 3.

4 Discussions and conclusions

We have investigated the impact of the preconditioner-parameter on the convergence rate of the preconditioned Chambolle-Pock algorithm applied to TV constrained data-divergence minimization. The investigations are performed in the context of 2D CT and 3D EPRI, corresponding to 2D and 3D inverse radon transforms, respectively. For 2D CT, we use the Shepp-Logan and two FORBILD phantoms; for 3D EPRI, we use a simulated 6-spheres phantom and a physical phantom. For each phantom evaluation, we set up 3 groups of balance parameters and then evaluate the preconditioner-parameter by tuning it from 0 to 2. Evaluations for all the 15 cases show that setting the preconditioner-parameter (α) to be 1 cannot necessarily achieve the fastest convergence and appear to suggest that one should tune the parameter in the range between 0 and 2 to achieve the optimal value under a specific imaging condition rather than simply set up it to be 1 as has been done in many previous implementations of the preconditioned Chambolle-Pock algorithm.

However, we find that, while setting α = 1 may not provide the fastest convergence, it always achieves an acceptable convergence rate. Selecting a value of α that is too low will likely result in a very slow convergence rate. The reason for this dichotomy may be that a low α leads to T that is too small and subsequently leads to image updates that are too small (see Equation (9) and line 8 of Table 1). Though α = 2 achieves the fastest convergence in several cases, it does not necessarily guarantee the fastest convergence. If the convergence-rate difference between the optimal α and α = 1 is small enough, one may always set α to be 1 for simplicity. One efficient and practical parameter-tuning method may be that we tune the parameters among three representative values (e.g., 1, 1.5 and 2), assuming that the value providing the fastest convergence may be safely implemented as the approximately optimal selection.

The ordinary CP (OCP) algorithm has convergence rate of O (1/N) [14, 16]. In the algorithm, δ is the step-size for increasing the dual objective, whereas, τ is the step-size for decreasing the primal objective, aiming at converging to the saddle-point, i.e. the minimal point of the primal objective. If the primal or dual objective is uniformly convex, then the convergence rate of CP algorithm may be improved to O (1/N²) by updating the step-size in each iteration. If the primal and dual objectives are both uniformly convex, then the algorithm may achieve linear convergence of O (ω^N) for ω < 1 by superior selection of the two step-sizes [14]. Actually, the PCP algorithm has much faster convergence rate than OCP algorithm by use of vector-step-size other than scalar-step-size [14]. Clearly, superior step-size selection may lead to faster convergence rate. In the work, we found that the preconditioner parameter has significant impact on convergence rate. Theoretically speaking, the further acceleration of the PCP algorithm by selecting superior preconditioner parameter also depends on the superior step-size selections. But we must stress that the optimal selection of the preconditioner parameter is object and imaging conditions dependent. For a specific imaging task, we may always achieve the optimal parameter by the determination method described in the work.

Though the investigations presented in this work are within the context of TVDM for image reconstruction, the observations may be extended to other optimization-based problems solved by PCP algorithm, such as image denoising, graph cuts, and etc. It may be suggested that one should tune the preconditioner-parameter so as to achieve the fastest convergence rate of the preconditioned Chambolle-Pock algorithm for solving the specific optimization program to which it is applied.

Footnotes

Acknowledgments

This work was supported in part by the Natural Science Foundation of Shanxi Province under grant 201601D011041, by the China Scholarship Council under grant 201608140014 and by the National Institutes of Health (NIH) under grants P41 EB002034, R01 s CA98575 and CA182264. The authors also would like to thank Prof. Xiaochuan Pan and Dr. Emil Sidky of The University of Chicago for their heuristic discussions.

References

Bian

, Siewerdsen

J.H.

, Han

et al., Evaluation of sparse-view reconstruction from flat-panel-detector cone-beam CT, Physics in Medicine and Biology 55 (2010), 6575–6599.

Sidky

E.Y.

and Pan

, Image reconstruction in circular cone-beam computed tomography by constrained, total-variation minimization, Physics in Medicine and Biology 53 (2008), 4777–4807.

Qiao

, Redler

, Epel

, Qian

and Halpern

, 3D pulse EPR imaging from sparse-view projections via constrained, total variation minimization, Journal of Magnetic Resonance 258 (2015), 49–57.

Omotavo

and Elbakri

, Objective performance assessment of five computed tomography iterative reconstruction algorithms, Journal of X-ray Science and Technology 24 (2016), 913–930.

Cai

, Wang

, Li

et al., Optimization-based image reconstruction in computed tomography by alternating direction method with ordered subsets, Journal of X-ray Science and Technology 25 (2017), 429–464.

Pan

, Sidky

E.Y.

and Vannier

, Why do commercial CT scanners still employ traditional, filtered back-projection for image reconstruction? Inverse Problems 25 (2009), 123009.

Chen

, Ma

, Feng

et al., Nonlocal prior bayesian tomographic reconstruction, Journal of Mathematical Imaging and Vision 30 (2008), 133–146.

, Huang

, Feng

et al., Low-dose computed tomography image restoration using previous normal-dose scan, Medical Physics 38 (2011), 5713–5731.

Donoho

D.L.

, Compressed sensing, IEEE Transactions on Information Theory 52 (2006), 1289–1306.

10.

Zhao

, Haldar

J.P.

, Christodoulou

A.G.

and Liang

Z.P.

, Image reconstruction from highly undersampled (k, t)-space data with joint partial separability and sparsity constraints, IEEE Transactions on Medical Imaging 31 (2012), 1809–1820.

11.

Sidky

E.Y.

, Kao

C.M.

and Pan

, Accurate image reconstruction from few-views and limited-angle data in divergent-beam CT, Journal of X-ray Science and Technology 14 (2006), 119–139.

12.

Esser

, Applications of Lagrangian-based alternating direction methods and connections to split Bregman, CAM Report 9 (2009), 31.

13.

Esser

, Zhang

and Chan

, A general framework for a class of first order primal-dual algorithms for TV minimization, UCLA CAM Report (2009), 9–67.

14.

Chambolle

and Pock

, A first-order primal-dual algorithm for convex problems with applications to imaging, Journal of Mathematical Imaging and Vision 40 (2011), 120–145.

15.

Sidky

E.Y.

, Jørgensen

J.H.

and Pan

, Convex optimization problem prototyping for image reconstruction in computed tomography with the Chambolle–Pock algorithm, Physics in Medicine and Biology 57 (2012), 3065.

16.

Chambolle

and Pock

, An introduction to continuous optimization for imaging, Acta Numerica 25 (2016), 161–319.

17.

Zhang

, Han

, Pearson

et al., Artifact reduction in short-scan CBCT by use of optimization-based reconstruction, Physics in Medicine and Biology 61 (2016), 3387.

18.

Sidky

E.Y.

, Kraemer

D.N.

, Roth

E.G.

et al., Analysis of iterative region-of-interest image reconstruction for x-ray computed tomography, Journal of Medical Imaging 1 (2014), 031007.

19.

Xia

, Langan

D.A.

, Solomon

S.B.

et al., Optimization-based image reconstruction with artifact reduction in C-arm CBCT, Physics in Medicine and Biology 61 (2016), 7300.

20.

Zhang

, Ye

, Chen

et al., Investigation of optimization-based reconstruction with an image-total-variation constraint in PET, Physics in Medicine and Biology 61 (2016), 6055.

21.

Pock

and Chambolle

, Computer Vision (ICCV), vol. Series: IEEE, 2011 IEEE International Conference on, 2011, pp. 1762–1769.

22.

Tang

, Cai

, Wang

and Peng

, Medical Imaging Physics and Engineering (ICMIPE), vol. Series: IEEE, 2013 IEEE International Conference on, 2013, pp. 6–10.

23.

Siddon

R.L.

, Fast calculation of the exact radiological path for a three-dimensional CT array, Medical Physics 12 (1985), 252–255.

24.

Joseph

P.M.

, An improved algorithm for reprojecting rays through pixel images, IEEE Transactions on Medical Imaging 1 (1982), 192–196.

25.

De Man

and Basu

, Distance-driven projection and backprojection in three dimensions, Physics in Medicine and Biology 49 (2004), 2463.

26.

Epel

, Sundramoorthy

S.V.

, Mailer

and Halpern

H.J.

, A versatile high speed 250-MHz pulse imager for biomedical applications, Concepts in Magnetic Resonance Part B: Magnetic Resonance Engineering 33 (2008), 163–176.

27.

Epel

, Gromov

, Stoll

et al., Spectrometer manager: A versatile control software for pulse EPR spectrometers, Concepts in Magnetic Resonance Part B: Magnetic Resonance Engineering 26 (2005), 36–45.

28.

Redler

, Epel

and Halpern

H.J.

, Maximally spaced projection sequencing in electron paramagnetic resonance imaging, Concepts in Magnetic Resonance Part B: Magnetic Resonance Engineering 45 (2015), 33–45.

29.

Qiao

, Redler

, Epel

and Halpern

H.J.

, Comparison of parabolic filtration methods for 3D filtered back projection in pulsed EPR imaging, Journal of Magnetic Resonance 248 (2014), 42–53.

Investigation of the preconditioner-parameter in the preconditioned Chambolle-Pock algorithm applied to optimization-based image reconstruction

Abstract

Keywords

1 Introduction

2 Methods

2.1 Imaging system modelling

3 Results

3.1 Evaluations of the CT reconstructions of 2D complete FORBILD phantom with two ears

Table 6 Iteration index achieving the stopping condition for different α in experimental data study Iteration index α = 0 α = 0.5 α = 1 α = 1.5 α = 2 λ = 1, ν = 0.1v A 110 93 84 82 87 λ = 1, ν = v A 286 230 199 179 168 λ = 1, ν = 10v A 909 738 637 568 518

Footnotes

Acknowledgments

References

Table 6
Iteration index achieving the stopping condition for different α in experimental data study

Iteration index α = 0 α = 0.5 α = 1 α = 1.5 α = 2

λ = 1, ν = 0.1v_A 110 93 84 82 87

λ = 1, ν = v_A 286 230 199 179 168

λ = 1, ν = 10v_A 909 738 637 568 518