Optimization methods for very accurate Digital Breast Tomosynthesis image reconstruction

07/20/2020 ∙ by Elena Morotti, et al. ∙ University of Bologna 0

Digital Breast Tomosynthesis is an X-ray imaging technique that allows a volumetric reconstruction of the breast, from a small number of low-dose two-dimensional projections. Although it is already used in clinical setting, enhancing the quality of the recovered images is still a subject of research. Aim of this paper is to propose, in a general optimization framework, very accurate iterative algorithms for Digital Breast Tomosynthesis image reconstruction, characterized by a convergent behaviour. They are able to detect the cancer object of interest, i.e. masses and microcalcifications, in the early iterations and to enhance the image quality in a prolonged execution. The suggested model-based implementations are specifically aligned to Digital Breast Tomosynthesis clinical requirements and take advantage of a Total Variation regularizer. We also tune a fully-automatic strategy to set a proper regularization parameter. We assess our proposals on real data, acquired from a breast accreditation phantom and a clinical case. The results confirm the effectiveness of the presented solutions in reconstructing breast volumes with particular focus on the masses and microcalcifications.



There are no comments yet.


page 10

page 12

page 14

This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

1 Introduction

Digital Breast Tomosynthesis (DBT) is a 3D X-ray cone-beam Computed Tomography (CT) technique for the early detection of breast tumors Wu et al. (2003); Andersson et al. (2008).

While the traditional digital mammography provides a unique 2D breast image, DBT reconstructs the breast as a stack of 2D images by using a comparable radiation dose. Hence DBT is also used in screening programs, because the volumetric reconstruction reduces the tissue overlaps allowing for a better visibility of malignant structures. DBT is characterized by a limited-angle geometry: since the object is scanned only from a narrow angular range, the DBT projection data is incomplete if compared to classical CT cases.

The reconstruction algorithm plays an important role, influencing the accuracy of the recovered breast images. It is well known that traditional fast analytic reconstruction methods, such as Feldkamp Feldkamp et al. (1984), produce poor noisy images in limited-angle tomography, hence they have been left in favour of Iterative Reconstruction (IR) algorithms Wu et al. (2004); Bliznakova et al. (2012); Reiser et al. (2009). IR solvers provide a sequence of solutions, by computing an improved reconstructed volume at each iteration. Many iterative reconstruction solvers have been proposed in literature. An overview of the IR methods is discussed in 2 and a good paper reviewing IR methods is Beister et al. (2012).

In this work we consider IR algorithms as solvers of a model-based formulation through an unconstrained optimization problem, where the objective function both describes the CT process by modelling the physics of the system (including the presence of noise on the projection data) and introduces some image priors. Such a mathematical approach is quite uncommon in 3D tomographic imaging, where a constrained formulation is preferred Sidky et al. (2009a); Choi et al. (2010); Ritschl et al. (2011); Luo et al. (2017). In particular we consider the objective function as the sum of the Least Squares (LS) data fitting term and the Total Variation (TV) regularization function. The TV regularizer is chosen by many authors because of its excellent shape recovering and denoising properties, even if it is known that it can produce staircasing effects when the regularization parameter is too high Liu (2014); Sidky et al. (2009a); Choi et al. (2010); Ritschl et al. (2011); Graff and Sidky (2015); Luo et al. (2017). Hence the choice of the regularization parameter plays a fundamental role in the model-based formulation.

Figure 1 shows how we approach the entire DBT imaging process, from the numerical modeling of the projection step during the breast scanning, to the reconstructed volume inspection looking for breast cancer objects, via the implementation of an iterative solver for the model-based minimization problem.

Aim of the paper is to propose both a TV-based optimization framework and three accurate iterative solvers which use accelerated first order strategies, for DBT image reconstruction. We are also interested in finding an automatic strategy to set a satisfactory regularization parameter and thus avoid its manually tuning which is infeasible in a clinical setting.

The contribution of this work can be summarized as follows.

  • We present three IR solvers in a unique optimization framework which can reconstruct clinically usable DBT images in few iterations as well as very accurate reconstructions if more iterations are allowed. Even if in clinical routine almost real time reconstructions are required, we remark the importance of improving the image quality with ongoing iterations in longer execution times, for two main reasons: first, having more reliable images can be crucial in difficult diagnosable cases to avoid false responds; second, the fast evolution of multiprocessor boards, such as GPUs, is drastically reducing the time per iteration of the methods, hence we can suppose that more iterations could be performed in clinical reconstructions in the next future.

  • We propose a user independent and computationally effortless rule to set and adapt the regularization parameter at each iteration of the algorithms.

  • In order to assess our proposals, we implement the methods and test them on real projection data of both a breast accreditation phantom and a human patient. We analyse the algorithms performance in recovering the breast tumor objects of interest, by means of measures of merits and visual inspection, at different stages of the iterative reconstruction process. We analyse the volume via its recovered slices, both perpendicularly and along the direction (see Figure 1).

The paper is organized as follows. We present an overview of IR methods in Section 2. In Section 3 we state the optimization framework for the image reconstruction, thus we illustrate the three proposed IR solvers in Section 4. Sections 5 and 6 present the data sets and the experimental results, respectively. Finally, Section 7 contains some conclusions.

Figure 1: Scheme of the DBT reconstruction process. On the left, a draft of the frontal (coronal) section of a DBT system acquiring projection images of the breast; in the centre, a chart representing the -th iteration of the algorithm computing the sequence of approximate solutions by solving the model-based minimization problem; on the right, the evaluation of reconstructed volumes by inspection of cancer objects of interest.

2 State of art

Iterative approaches have been introduced since the first years of CT, but they have not been used for long time due to their high computational time request. Recently, IR methods got a renewed interest in scientific communities and among the major vendors, due to the advent of more performing processors Beister et al. (2012). As a consequence, a wide amount of IR methods has been proposed to reconstruct tomographic images and an exhaustive analysis can be found in Graff and Sidky (2015).

Initial efforts to solve tomographic imaging with IR methods took an algebraic approach. Algorithms such as ART, SIRT, SART and their modifications iteratively solve a linear system of equations by sequentially projecting a solution onto different hyperplanes

Kak and Slaney (1988).

On the other hand, the worldwide increasing interest in Compressive Sensing (CS) Donoho (2006) promoted a novel model-based iterative approach, which uses an optimization framework to exploit CS theory. Among the wide class of model-based IR methods, the so called Sparsity-Exploiting Image Reconstruction (SEIR) methods have produced significant improvement to the image quality in all the low-dose CT applications (see Liu (2014); Graff and Sidky (2015) and references therein). In particular, many authors introduce the TV function to take advantage of the sparsity in the image gradient domain for edge detection Wu et al. (2004); Bliznakova et al. (2012); Sidky et al. (2006); Sidky and Pan (2008); Sidky et al. (2008, 2009b); Lu et al. (2010); Chen et al. (2013); Ertas et al. (2013); Rose et al. (2015). This property turns into practise as a noise smoothing effect and as a reliable detection of shape and size of anatomical objects (such as microcalcifications and masses), which are fundamental tasks of DBT imaging.

It is possible to distinguish two main categories of algorithms in the class of SEIR methods: the approximate solvers and the accurate solvers. The first one contains algorithms which use, at each step, an algebraic approach (such as SART and SIRT) sequentially and then decrease the TV of the just calculated solution. Examples are the well-known POCS algorithm and its developments Sidky and Pan (2008); Chen et al. (2008); Ramirez-Giraldo et al. (2011). They provide reliable reconstructions in few iterations, but the quality of the recovered images strongly depends on the tuning of many inner parameters and the algorithm convergence is not guaranteed.

On the other hand, the accurate solvers are optimization methods which minimize an objective function defined as a sum of a fit-to-data term and a regularization function. The two quantities are typically weighted by a regularization parameter. This class is represented by classical optimization methods adapted to the huge size 3D tomographic reconstruction problems. Their solution is proved to converge to the exact solution of the minimization problem.
Nowadays, only preliminary investigations on simulations or phantoms have been performed to analyse the results of accurate solvers for few-views CT applications Jensen et al. (2012). In our previous works we have investigated a Fixed Point (FP) algorithm and a Scaled Gradient Projection (SGP) method in Loli Piccolomini and Morotti (2016) and Loli Piccolomini et al. (2018) respectively, and we have applied them to small simulated data sets. A parallel SGP implementation on GPUs has been presented in Cavicchioli et al. (2020), where we were interested in showing the computational efficiency of the algorithm in terms of execution time. We also remark that in Park et al. (2012) an accelerated Gradient Projection method outperformed a sequential reconstruction algorithm on data acquired in a sub-sampled 2D circular geometry. A further example is the Chambolle-Pock (CP) algorithm which has been applied in Sidky et al. (2014) onto 2D circular geometry breast CT, to solve a TV-based convex optimization problem. In this work we consider the three iterative optimization solvers, namely the SGP, FP and CP, to further evaluate their feasibility in reconstructing DBT real volumes and recovering breast tumor objects, like masses and microcalcifications.

Concerning existing rules for the regularization parameter choice in tomography, in Niinimäki et al. (2016) the authors propose a strategy based on multiresolution and apply it to 2D reconstructions. The proposed rule is very promising, but it is quite expensive for a very large size 3D application, such as DBT image reconstruction. An exhaustive list of existing rules for the selection of the regularization parameter is reported in Niinimäki et al. (2016).

3 The optimization framework in model-based formulation

Mathematically, tomographic image reconstruction is an inverse ill-posed problem whose solution can be obtained by minimizing a suitable objective function related to the physical process. To define the model describing image reconstruction is therefore crucial a deep understanding of the acquisition steps characterizing the DBT technique. A schematic example of a DBT system is shown on the left of Figure 1. In DBT routine, the breast is first compressed along the -axis, over the flat detector plane. The source moves along an arc trajectory and emits low-dose radiations from a discrete number of angles. Once the X-ray cone-beam has passed through the body, the detector records its attenuation: the set of the resulting projection images constitutes the raw tomographic data set. The breast volume to be recovered is composed by a stack of high resolution images, parallel to the detector plane along the vertical direction.
In order to define the numerical model of tomographic image formation, we discretize the 3D object into voxels, whereas the 2D detector panel is made of recording units. For each fixed projection angle and -th detector recording unit, the Lambert-Beer law relates the projections , along a ray , to the attenuation coefficient function of the voxel crossed by Epstein (2007) as:


where represents the intensity of the energy emitted by the X-ray source. The discretization of the integral in (1) for all the scanning angles arises the following linear system:


In equation (2) we denote with the

dimensional vector stacking the attenuation coefficients of all the voxels, while

is the vector of size storing all the projections (i.e. the right hand sides of (1)) and is the matrix of size , built according to the DBT device geometry and representing the projection process onto the detector.

Some issues arise when solving the linear system (2) as an inverse problem, such as the existence of infinite solutions (since ) and the presence of high noise in the reconstructed images (due to the ill-posedness of the problem). The model-based approach is introduced to overcome these numerical controversies, by adding some a priori information. The resulting formulation can be stated as an unconstrained or constrained minimization problem Graff and Sidky (2015). We consider here the former problem and express it as:


where is a fit-to-data function, is the prior function (acting here as a regularizer) and is the regularization parameter.

To such DBT mathematical formulation, we can add the box constraint reflecting the non-negativity property of the linear attenuation coefficient in (1).

In particular, in this work we settle as the Least Squares (LS) function


and as the Total Variation (TV) operator defined as Vogel (2002):


Since TV is not differentiable in the origin, in the algorithms requiring the computation of the gradient, we consider its smoothed version:


where is a small positive parameter Vogel (2002). Exploiting the linearity of (3), the objective function gradient can be evaluated by separately computing as


and through finite forward differences.

4 Iterative optimization methods

To solve the minimization problem (3), we propose three accurate solvers the Scaled Gradient Projection (SGP), the Chambolle-Pock (CP) and the Fixed Point (FP) methods. For all these methods the convergence to the solution of the model-based minimization problem d Among the wide class of optimization methods they have been chosen since they satisfy the requirements necessary to be usable on DBT devices:

  • a fast error decreasing in the initial algorithm execution, in order to obtain a good image in few iterations;

  • a low computational cost per iteration (which is mainly determined by the number of matrix-vector products), to efficiently run the solver in short time;

  • a limited request of memory, to solve real-size problems on commercially affordable hardware.

A challenging issue (common to the implementation of the three algorithms) is the computation of the projection matrix : since it can not be stored due to its huge dimensions, it must be recalculated at each call. Thus, this section ends with a focus on the algorithm we use to generate .

4.1 Scaled Gradient Projection algorithm

The SGP algorithm is a first order accelerated method. We apply it to solve the non-negative constrained optimization problem:


Algorithm 1 reports the main steps of the SGP algorithm.

At each -th iteration, the new solution is computed by moving along a descent direction of a quantity , as:


The direction is obtained through a projection onto the non-negative orthant:


where is the step length and is the scaling matrix (step 7 in Algorithm 1).
Essentially, the method follows a Gradient Projection approach accelerated by choosing the step length with Barzilai-Borwein techniques and by introducing a suitable scaling matrix improving the matrix conditioning Loli Piccolomini et al. (2018). In particular, the scaling matrix is a diagonal matrix with entries in a limited interval. To update (line 5 of Algorithm 1), we compute a splitting of the objective function gradient into its positive and negative parts, as:


where and . The diagonal elements of are updated, for as:


where is a decreasing positive sequence.

Regarding the convergence, it is proved in Bonettini and Prato (2015) that the SGP algorithm converges without any further restriction on the step length and on the scaling matrix to the unique minimum of (8). In Bonettini et al. (2016), the authors proved that the theoretical convergence rate of the SGP method is (1/k).

4:while not convergence do
5:     Compute
6:     Compute
7:     Define with alternate BB rules
10:     while   do
13:     k = k+1
Algorithm 1 Scaled Gradient Projection algorithm (SGP)

4.2 The Fixed Point algorithm

The FP algorithm for the solution of the minimization problem:


has been firstly proposed for image denoising by Rudin, Osher and Fatemi in Rudin et al. (1992). Starting from this approach we derived the lagged diffusivity FP Algorithm 2 for 3D tomographic image reconstruction.

3:for  to  do
4:     Compute
5:     Solve the linear system , where , with the Conjugate Gradient method.
Algorithm 2 Lagged diffusivity Fixed Point algorithm (FP)

At each -th iteration, the FP algorithm updates the solution with the following rule:


where the descent direction is computed by solving a linear system (line 4 of Algorithm 2). The matrix in line 4 approximates the Hessian matrix. It contains in fact the seven diagonals banded matrix which is the discretization matrix of the diffusion operator so that Vogel (2002). We solve the linear system with very few iterations of a Conjugate Gradient (CG) algorithm Hestenes et al. (1952): we stop it far before convergence, both to limit the computational time and to prevent noise from affecting the solution. We remark that each CG iteration requires a matrix-vector product involving and that, to save memory space, we perform it without storing the matrix : we only store and re-compute and at run time. At the end, we project the last computed solution onto the non-negative orthant. For more details on the FP method applied to tomographic image reconstruction and its convergence, see Loli Piccolomini and Morotti (2016) and Chan and Mulet (1999) respectively.

4.3 The Chambolle-Pock algorithm

The CP algorithm has been firstly proposed in Chambolle and Pock (2011) for the solution of the general minimization problem:


where is a convex, lower-semicontinuous and proper function, is convex and lower-semicontinuous and K is a continuous linear operator. To fit the problem statement (15) and fully exploit the linearity of , we assign:


which results in defining the operator as a matrix composed by the four following blocks:


where and are the forward differences operators acting along the and axes respectively. In order to include the non-negative constraints, we fixed as the indicator function of the convex set , i.e.


Considering the convex conjugate of , defined as , and the proximal mappings of and , i.e.


the -th CP iteration can be described with the following three steps:

  1. compute as ;

  2. compute as ;

  3. define with an extrapolation step: and .

In particular, the proximal mapping can be computed as sum the two independent blocks, as in lines 5-6 and 7-8 of the Algorithm 3; its detailed derivation can be found in Sidky et al. (2014). The proximal mapping of is defined as:


hence it is exactly the projection of onto the feasible set (lines 9-10 of Algorithm 3). The updated iterate is computed with a FISTA strategy as in line 11 of Algorithm 3. The algorithm convergence is demonstrated in Chambolle and Pock (2011).
We finally observe that the algorithm needs to compute the value (line 1 of Algorithm 3

): to estimate the matrix 2-norm as


is the spectral radius of a matrix), we perform two iterations of the power method for the maximum eigenvalue computation

Golub and Van Loan (2012).

2:Compute: as an approximation of
3:Initialize: ,
4:Initialize: and to zeros-vectors
5:for  to  do
Algorithm 3 Chambolle Pock algorithm (CP)

4.4 User-independent choice of the regularization parameter

In model-based optimization approach (3), the choice of the regularization parameter plays a key role for the quality of the reconstruction and it represents a crucial challenge in a clinical setting, where the trial-and-error approach is not doable for each reconstruction. Moreover, experimental results show that a strong regularization is required in the first iterations, to avoid noise propagation and force the algorithm towards a good solution, whereas a weaker regularization in the last iterations can prevent the TV staircaising effects on final reconstructions. Hence, we propose to reduce the regularization weight along the iterations, by choosing the values with a decreasing updating rule. Interestingly, state of art studies have already proposed semi-automatic rules for the selection of a decreasing sequence of regularization parameters defining a sequence of minimization problems stated as (13), whose solutions converge to a good reconstructed image. See Lazzaro et al. (2019) for more details and the convergence proof.

We propose the following fully-automatic strategy to compute a decreasing sequence . At the beginning of our algorithm, we leave out the regularization by setting the first parameter : we are in fact interested in a very good data fitting, to recover as many image features as possible. Next, the starting value is set to balance the residual norm and the amount of TV of the first iterate. Afterward, we propose to decrease of a constant factor at each -th iteration, since we need a very simple and computationally cheap rule, reducing the regularization weight slightly. The resulting strategy is summarized in the following scheme and it can be introduced in each of the previously considered algorithms.

  • Set to initialize the algorithm and run the first iteration (labelled with k=0) to compute ;

  • Set and use it to compute ;

  • For each , set


    and use it to compute .

4.5 The projection matrix algorithm

Besides the choice of the model parameter and the solver, in optimization approach a key point consists in numerical modeling the geometric projection process, schematically displayed from a frontal view in Figure 1, through a matrix.

The coefficient matrix of the linear system (2) is commonly called projection operator in tomography, since it represents the action of the tomographic system in projecting an object onto the detector, whereas the matrix modeling the backprojection of the tomographic data onto a volume is called backprojection operator. In the proposed optimization algorithms, the backprojection coincides with the transpose matrix .
Different algorithms have been proposed in literature for the computation of the matrix . We have adopted the Distance Driven (DD), which accurately models the discretization of the Lambert-Beer’s law (1) for cone-beam projections De Man and Basu (2004). In DD, of size is constituted by submatrices of size . Each element represents the contribution of the -th voxel (for ) to the projection onto the -th detector pixel (for ), for a projection angle . Images in Figure 2 help in understanding the DD procedure. In Figure 2 (a), for a scanning angle, we consider the X-ray cone-beam projecting onto the -th blue pixel and intersecting the voxels with bold contours (in the magenta coloured area). Only these voxels contribute to the value of the projection in the considered pixel. In Figure 2 (b) we highlight the i-th cell of the detector (the blue area) and its backward footprint on a plane parallel to the detector (the magenta area). The ratio between the magenta area inside the -th voxel and the whole magenta extension is proportional to the value of the matrix. For all the voxels not contributing to the -th projection the corresponding matrix element ; hence is extremely sparse. However, despite the huge number of nonzero elements, for its very large size, in real applications cannot be stored and it must be recomputed whenever a matrix-vector product is needed.

We finally remark that we have modified the general approach presented in De Man and Basu (2004) by efficiently exploiting the characteristics of our specific mammographic setting. Really, since the DBT detector is a stationary flat panel and it is parallel to the compression plane of the breast, the footprints can be directly projected onto the detector plane, thus avoiding the use of an intermediate projection plane and further computational costs.

Figure 2: Schematic draw representing the Distance Driven approach to compute the system matrix. (a) View on the plane of an X-ray projection onto a single pixel, from a fixed angle. The intersection of the X-ray beam with the volume is highlighted in magenta. (b) The magenta area represents the backward projection of the blue recording unit onto a volume slice parallel to the plane.

5 Materials

5.1 DBT system configuration

Our tests are performed on the digital system Giotto Class of the Italian I.M.S. Giotto Spa company in Bologna 21. The source executes scans from equally spaced angles in approximately degrees range; in the highest vertical position, the source is about over the detector. The stationary digital detector has a sensitive area of and squared pixel pitch of 0.085 mm; the reconstructed voxel dimensions along the three cartesian axes are and respectively.
The system uses a polychromatic ray with energies in a narrow range around 20 keV to avoid the photon scattering. As always happens in CT reconstruction algorithms, we approximate the polychromatic beam with a monochromatic one.

5.2 Data sets

We consider two data sets in our experiments: a breast 3D phantom and a clinical acquisition from a human subject. Both volumes contain the objects of interest for breast cancer detection, i.e. small high contrast microcalcifications and larger but lower contrasted masses.

The phantom is the model 020 of BR3D breast imaging phantom, produced by CIRS Tissue Simulation and Phantom company 12. It is characterized by a heterogeneous background, where adipose-like and gland-like tissues are mixed in about 50/50 ratio and it is made of six slabs that may be arranged to create multiple anatomical backgrounds. Each slab has a semicircular shape and its size is . Inside one of them, we find acrylic spheres simulating breast masses (MSs), 1 length fibers and many clusters of calcium carbonate specks simulating microcalcifications (MCs). We report in Table 1 the length of the diameters of all the MSs and of each sphere of a MC cluster. In particular, we reconstruct a volume of 50 slices of and we analyze in the reconstructed images objects with different diameters, such as the microcalcifications in clusters 3, 5 and 6 (having 230, 165 and 130 diameter, respectively) and the second and fourth mass (with diameter 4.7 and 3.1 , respectively). All such objects lye on the same slice.

As human DBT data set, we have chosen a case containing microcalcifications, circular and spiculated masses. The clinical volume is constituted of 55 slices of .

1 2 3 4 5 6
MC 400 290 230 196 165 130
MS 6300 4700 3900 3100 2300 1800
Table 1: Diameters of a microcalcification (MC) in a cluster and of a mass (MS) in the BR3D phantom as reported in 12. Measures are in micrometers ().

5.3 Measure and graphics of merits

In order to quantitatively evaluate the reconstructed objects of interest in the volumes, we compute two widely used measure of merits: the Contrast-to-Noise Ratio (CNR) and the Full-Width at Half Maximum (FWHM).

The CNR measure on a mass is calculated as:


where and

are the mean and standard deviation computed on the reconstructed volume, in small regions located inside the mass (MS) or in the background (BG). Similarly, we define the CNR measure on a microcalcification as:


where is the maximum intensity inside the considered microcalcification (MC). Higher values of the CNR indices reflect a better detection of an object from the background.

To compute the FWHM parameter, we consider the transverse slice (parallel to the plane) where the microcalcification lies and we extract the Plane Profile (PP) along the axes. The FWMH index is computed as:


where is the standard deviation of the gaussian curve fitting the PP. We remark that


approximates the width of the examined microcalcification. The Plane Profiles are also useful tools to evaluate the reconstruction accuracy on the transverse plane.

To estimate the solver effectiveness along the direction, which is the most challenging purpose in DBT imaging, we plot the Artifact Spread Function (ASF) vector, whose components are computed on a microcalcification as:


where is the mean of the reconstructed values inside a circular region of three pixels diameter inside the considered MC and in the background, corresponds to the slice where the object is on focus and is the total number of discrete slices. Similarly, we compute the ASF for the masses.

6 Numerical results and discussion

In this section we present the results obtained with the proposed optimization approach and the accurate solvers described in section 4. At first we compare the BR3D phantom reconstructions produced by the SGP, FP and CP solvers in a similar computational time. Then, to analyse the best obtainable image quality we have run the SGP up to convergence both on the phantom and on the clinical data set. In all the previous tests we used a constant value of , set by trial and error. At last, we test the automatic rule proposed in section 4.4 to decrease the values along the iterations.

6.1 Methods comparison for early reconstructions

Aim of this paragraph is to show the behaviour of the proposed solvers at different stages of their executions. We fixed 5 and 15 iterations: the workload of 5 iterations is compatible with the execution of a reconstruction on a commercial hardware in a clinical setting, while in 15 iterations we get fairly accurate reconstructions with all the three methods, reflecting that they are sufficiently close to the convergence solution. Each SGP and CP iteration requires approximately the same time, whereas in the special case of FP solver, the number of allowed iterations corresponds to the sum of the external and CG iterations. Since we perform 4 CG iterations, we have stopped the FP algorithm after one or three outer iterations (loop in Algorithm 2), respectively.

The value of in (6) has been fixed as . Since the three considered algorithms solve a slightly different optimization problem, the three parameters have been chosen independently for each method to achieve the best reconstruction in 5 iterations: we have set for both SGP and CP methods and for the FP algorithm.

In the following analysis, we focus on the reconstruction of MC cluster number 3 in the BR3D phantom. For each solver, Figure 3 reports a pixels crop taken from the fifth slice of the reconstructions in 5 and 15 iterations. The images are represented by automatically enhancing the gray level contrast computed on the same considered region. In Figure 4 we compare the PP and ASF curves taken on one MC.

Looking at Figure 3, we observe that the detection of the MC cluster is comparable at equal iterations whereas the background appears slightly different for the three methods. For example, in the case of CP reconstruction in 15 iterations it looks smoother and more blurred. Focusing on the objects of interest, we notice that in 5 iterations the MCs are perfectly visible; moreover, in 15 iterations the MC edges are sharper as confirmed by the plots (a) and (c) in Figure 4. From plots (b) and (d) of Figure 4 we observe that in all the three reconstructions the object is placed in the correct slice and it is not diffused in the adjacent layers. Hence we can conclude that the proposed model-based optimization framework yields good quality images in early reconstructions, regardless the applied solver.

(a) SGP
(b) SGP
(c) FP
(d) FP
(e) CP
(f) CP
Figure 3: Reconstructions of microcalcification cluster number 3 in BR3D phantom obtained with SGP, FP and CP methods. On the left column, reconstructions in 5 iterations; on the right, reconstructions in 15 iterations.
(a) PP in 5 iterations
(b) ASF in 5 iterations
(c) PP in 15 iterations
(d) ASF in 15 iterations
Figure 4: Plots of the Plane Profile on the left and of the ASF vectors on the right, taken over one microcalcification of cluster number 3 in BR3D phantom obtained. In all the plots: the red line corresponds to SGP method, the blue line to FP method and the green line to CP method.

6.2 SGP algorithm insights

Figure 5: Objective function values vs. iteration number for the SGP execution on the phantom test. The convergence has been reached after 44 iterations by satisfying condition (27). The red labels outline the function values at 5, 15 and 30 iterations.

In the following, we raise the SGP as the representative solver in the proposed optimization framework. We explore the performance of the optimization approach on many different objects of the BR3D phantom (such as microcalcifications of very small diameter and masses with a low contrast with the background tissue) and we analyse the quality of the reconstructions also after 15 iterations, i.e. approaching convergence.

In fact, we run the SGP solver on the BR3D phantom until the stopping condition


is satisfied. It occurs after 44 iterations. In Figure 5, we plot the objective function values vs. the number of iterations: we observe that the objective function fast decreases in the first 5 iterations, whereas it exhibits a very flat trend from 10 iterations on, as it is confirmed by the red labelled values. We have seen, in fact, that the reconstructed images are visually almost indistinguishable after 30 iterations.
In Figure 6 we exhibit the reconstructions of the 165 m MCs of cluster 5, and the 4.7 mm mass (MS 2), obtained by the SGP algorithm after 5, 15 and 30 iterations. In Figure 7 we report the corresponding PP and ASF plots. Figure 6 (a) shows that the MC of cluster 5 can be clearly visible after only 5 iterations and the PP plots of Figure 7 (a) confirms that the it gets more and more enhanced from the background. The ASF plot in Figure 6 (c) shows an improvement in the object detection along the direction. As visible in Figure 6 and from the PP plot of Figure 7 (b), MS 2 is out of focus at 5 iterations but its contours are more and more defined when the algorithm approaches to convergence. The previous plots confirm that the proposed model with TV regularization is more effective in recovering high contrast objects such as microcalcifications than low absorbing structures such as masses.

In Table 2 we report the values of the CNR parameter on the examined reconstructed microcalcifications and masses. In particular, recalling the CNR definition (22) for the masses, the background area is a circle with diameter of 80 voxels, whereas we have considered circles of diameter 40 and 25 voxels inside the masses 2 and 4, respectively. When we compute the CNR value on a MC with equation (23) we consider the background as a circle of 20 voxels diameter and we compute on a small circle of diameter 5 voxels containing the microcalcification. Table 3 shows the values of the FWHM index defined in (24) and computed on one of the reconstructed microcalcification in each cluster. The corresponding MCs width computed as in (25) in micrometers are reported to be compared to the values of the diameters of the actual objects, shown in Table 1. Both tables demonstrate that we can get improved and more accurate reconstructions as the SGP approaches to the convergence: the increasing CNR indexes exhibit good denoising effects whereas the object enhancement is confirmed by the FWHM decreasing values. We remark that the MCs of cluster 6 are not discernible from the background in only 5 iterations (the FWHM is not measurable on the sixth MC cluster), because they are 130 m width and they should approximately fill inside only two voxels. However, they can be well recovered after more iterations with a good approximation of their real size.

(a) 5 iterations
(b) 5 iterations
(c) 15 iterations
(d) 15 iterations
(e) 30 iterations
(f) 30 iterations
Figure 6: SGP results on BR3D phantom. (a)-(c) Reconstructions of MC cluster number 5 obtained after 5, 15 and 30 iterations. (f)-(h) Reconstructions of mass number 2 obtained after 5, 15 and 30 iterations.
(a) Plane Profile
(b) Plane Profile
(c) ASF
(d) ASF
Figure 7: SGP results on BR3D phantom. (d)-(e) Plane and Depth profile on one microcalcification of cluster 5. (i)-(j) Plane Profiles and ASF profile on the mass. In all the plots: black line corresponds to 5 iterations, red line to 15 iterations and blue line to 30 iterations.
5 it. 15 it. 30 it.
MC cluster 3 24.21 33.34 38.00
MC cluster 5 10.03 19.00 28.00
MC cluster 6 7.27 11.02 17.00
MS 2 0.82 1.07 1.66
MS 4 0.87 1.00 1.33
Table 2: Values of the CNR index computed after 5, 15 and 30 SGP iterations. The CNR value computed on microcalcifications is defined as in (23), whereas the CNR value computed on the masses is defined as in (22).
cluster 5 it. 15 it. 30 it. 5 it. 15 it. 30 it.
3 4.77 3.32 2.70 430 299 243
5 3.52 2.65 2.32 317 238 209
6 - 2.05 1.52 - 185 137
Table 3: FWHM index (24) and measures (25) computed on the reconstructed MCs of the BR3D phantom, after 5, 15 and 30 SGP iterations.

6.3 Experiments on a human data set

We now illustrate the results obtained by reconstructing a real breast volume with the SGP solver at different iterative stages. In Figure 8 (a)-(c) we report a crop of a reconstructed slice, where we can distinguish objects of interest, i.e. a spherical mass and a small microcalcification. The plots in Figure 8 (d)-(e) represent the PP calculated on the mass and the microcalcification, respectively. The mass is well distinguishable since the earliest reconstruction and its shape and gray level intensity do not change remarkably; however the regular blue and thin profile in Figure 8 (d) points out the denoising effects of the TV function in the last iterations. Also the microcalcification is detected in few iterations, even if a more time-consuming SGP execution enhances the contrast of the object with respect to the background. Table 4, reporting CNR and FWHM values computed on the objects in Figure 8, gives more insight on the quality of the reconstructions. In particular, it confirms that the noise progressively decreases and the microcalcification gets more and more defined, from 5 to 30 iterations.

In Figure 9, we report the reconstruction of two spiculated masses, which can occur in clinical cases. For such breast objects the previous measures of merits are not applicable. However we can observe that they both are well recognizable in the earliest reconstruction and the edges become sharper with increasing iterations.

(a) 5 iterations
(b) 15 iterations
(c) 30 iterations
(d) Plane Profile on the mass
(e) Plane Profile on the MC
Figure 8: Results obtained after 5, 15 and 30 SGP iterations on a human breast data set. (a)-(c) Reconstructions of a 440 400 pixels region presenting both a spherical mass (pointed by the arrow) and a microcalficication (identified by the circle). (d)-(e) Plane profiles on the mass and on the microcalcification. In the plots: black line corresponds to 5 iterations and blue line to 30 iterations.
5 it. 15 it. 30 it. 5 it. 15 it. 30 it.
MS 0.239 0.381 0.558 - - -
MC 8.78 16.59 16.49 8.57 7.81 7.29
Table 4: Values of CNR and FWHM measures on the mass and the microcalcification observable in Figure 8 (a)-(c).
(a) 5 iterations
(b) 15 iterations
(c) 30 iterations
Figure 9: Results obtained after 5, 15 and 30 SGP iterations on a human breast data set. The reported 558 480 pixels crops present two spiculated masses.

6.4 Experiments with a variable regularization parameter

In all the above experiments, we have set a constant value of the regularization parameter along the iterations, to let the solvers perform at their best on the prefixed model derived by the settled . In this paragraph we show the results achieved with the automatic rule (21) for the choice of a decreasing sequence applied to the SGP solver, to reconstruct the BR3D phantom.

Figure 10 plots the sequence with the blue line, while the red line represents the constant

value used in the SGP implementation in the previous experiments. We observe that the proposed strategy computes values greater than the heuristically fixed one

until the fifth iteration. In Figure 11 we compare the PP of one microcalcification from cluster 3, reconstructed at 5 and 15 iterations using both a fixed value and the proposed strategy for the regularization parameter. We can infer from Figure 11 (a) that the resulting larger TV weights in the first iterations produce more accurate results. However, on advanced reconstructions the differences are negligible. We can conclude that the proposed automatic strategy results very efficient.

Figure 10: Sequence of decreasing values versus the number of iterations (blue line) in the SGP execution on the phantom test. The red straight line represents the constant value used in SGP for the experiments presented in the previous sections.
(a) 5 iterations
(b) 15 iterations
Figure 11: Plane Profiles on one microcalcification of cluster number 3 of the phantom, obtained with SGP with different regularization parameters, in 5 and 15 iterations. In all the plots: red line corresponds to fix parameter, blue line the adaptive choice of .

7 Conclusions

In this paper, we have presented a general optimization framework including a TV regularized formulation for DBT image reconstruction. We have also proposed a user-independent rule for selecting suitable values of the regularization parameter.

The results obtained with three solvers are encouraging. In early reconstructions, objects of interest of size greater than 150 are visible and correctly located in the volume, whereas the object detection quality improves and the noise drastically reduces if more iterations are allowed. When extending the computation from 5 to 30 algorithms iterations, the increasing rate of the CNR value lies in a range to . At last, we have shown that varying the regularization parameter along the iterations produces better results, especially in the early stage of the algorithm execution, when compared to the use of a fixed value heuristically chosen.
Since the three considered solvers produce comparable high quality reconstructions, we can conclude that the proposed optimization problem statement can be successfully used to detect the most interesting objects in an early diagnosis of breast tumor.

8 Acknowledgments

The research has been funded by the Indam GNCS grant 2020 Ottimizzazione per l’apprendimento automatico e apprendimento automatico per l’ottimizzazione.


  • I. Andersson, D. M. Ikeda, S. Zackrisson, M. Ruschin, T. Svahn, P. Timberg, and A. Tingberg (2008) Breast tomosynthesis and digital mammography: a comparison of breast cancer visibility and birads classification in a population of cancers with subtle mammographic findings. European radiology 18 (12), pp. 2817–2825. Cited by: §1.
  • M. Beister, D. Kolditz, and W. A. Kalender (2012) Iterative reconstruction methods in x-ray ct. Physica Medica 28, pp. 94–108. Cited by: §1, §2.
  • K. Bliznakova, Z. Bliznakov, and I. Buliev (2012) Comparison of algorithms for out-of-plane artifacts removal in digital tomosynthesis reconstructions. Computer Methods and Programs in Biomedicine 107 (1), pp. 75 – 83. Note: Advances in Biomedical Engineering and Computing: the MEDICON conference case External Links: ISSN 0169-2607, Document, Link Cited by: §1, §2.
  • S. Bonettini, F. Porta, and V. Ruggiero (2016) A variable metric inertial method for convex optimization. SIAM J. Sci. Comput 31 (4), pp. A2558–A2584. Cited by: §4.1.
  • S. Bonettini and M. Prato (2015) New convergence results for the scaled gradient projection method. Inv. Probl. 31 (9), pp. 1196–1211. Cited by: §4.1.
  • R. Cavicchioli, J. Hu, E. Loli Piccolomini, E. Morotti, and L. Zanni (2020) A first-order primal-dual algorithm for convex problems with applications to imaging.gpu acceleration of a model-based iterative method for digital breast tomosynthesis. Scientific Reports 10 (1), pp. 120–145. External Links: Document Cited by: §2.
  • A. Chambolle and T. Pock (2011) A first-order primal-dual algorithm for convex problems with applications to imaging.. Journal of Mathematical Imaging and Vision 40 (1), pp. 120–145. Cited by: §4.3, §4.3.
  • T. F. Chan and P. Mulet (1999) On the convergence of the lagged diffusivity fixed point method in total variation image restoration. SIAM journal on numerical analysis 36 (2), pp. 354–367. Cited by: §4.2.
  • G. Chen, J. Tang, and S. Leng (2008) Prior image constrained compressed sensing (piccs). In Photons Plus Ultrasound: Imaging and Sensing 2008: The Ninth Conference on Biomedical Thermoacoustics, Optoacoustics, and Acousto-optics, Vol. 6856, pp. 685618. Cited by: §2.
  • Z. Chen, X. Jin, L. Li, and G. Wang (2013) A limited-angle CT reconstruction method based on anisotropic TV minimization.. Physics in medicine and biology 58 (7), pp. 2119–41. External Links: Document, ISSN 1361-6560 Cited by: §2.
  • K. Choi, J. Wang, L. Zhu, T. Suh, S. P. Boyd, and L. Xing (2010) Compressed sensing based cone-beam computed tomography reconstruction with a first-order method.. Medical physics 37 9, pp. 5113–25. Cited by: §1.
  • [12] Computerized Imaging Reference Systems. Note: https://www.cirsinc.com/products/a11/51/br3d-breast-imaging-phantom/BR3D Breast Imaging Phantom, Model 020 Cited by: §5.2, Table 1.
  • B. De Man and S. Basu (2004) Distance-driven projection and backprojection in three dimensions. Phys. Med. Biol.. Cited by: §4.5, §4.5.
  • D. Donoho (2006) Compressed sensing. IEEE Trans. inf. theory 52 (4), pp. 1289–1306. Cited by: §2.
  • C. L. Epstein (2007) Introduction to the mathematics of medical imaging. SIAM. Cited by: §3.
  • M. Ertas, I. Yildirim, M. Kamasak, and A. Akan (2013) Digital breast tomosynthesis image reconstruction using 2d and 3d total variation minimization. BioMedical Engineering OnLine 12 (1), pp. 112. External Links: Link Cited by: §2.
  • L. A. Feldkamp, L. C. Davis, and J. W. Kress (1984) Practical cone-beam algorithm. J. Opt. Soc. Am. A 1 (6), pp. 612–619. External Links: Document Cited by: §1.
  • G. H. Golub and C. F. Van Loan (2012) Matrix computations. Vol. 3, JHU press. Cited by: §4.3.
  • C. Graff and E. Sidky (2015) Compressive sensing in medical imaging. Appl. Opt. 54 (8), pp. C23–C44. Cited by: §1, §2, §2, §3.
  • M. R. Hestenes, E. Stiefel, et al. (1952) Methods of conjugate gradients for solving linear systems. Journal of research of the National Bureau of Standards 49 (6), pp. 409–436. Cited by: §4.2.
  • [21] IMS Giotto Class. Note: http://www.imsgiotto.com/Accessed: 2020-05-25 Cited by: §5.1.
  • T. L. Jensen, J. H. Jørgensen, P. C. Hansen, and S. H. Jensen (2012) Implementation of an optimal first-order method for strongly convex total variation regularization. BIT Numer. Math. 52 (), pp. 329–356. Cited by: §2.
  • A. Kak and M. Slaney (1988) Principles of computerized tomographic imaging ieee press. New York. Cited by: §2.
  • D. Lazzaro, E. L. Piccolomini, and F. Zama (2019) A nonconvex penalization algorithm with automatic choice of the regularization parameter in sparse imaging. Inverse Problems. External Links: Document Cited by: §4.4.
  • L. Liu (2014) Model-based iterative reconstruction: a promising algorithm for today’s computed tomography imaging. Journal of Medical Imaging and Radiation Sciences 45 (2), pp. 131 – 136. External Links: ISSN 1939-8654, Document Cited by: §1, §2.
  • E. Loli Piccolomini and E. Morotti (2016) A fast TV-based iterative algorithm for digital breast tomosynthesis image reconstruction. J. Algor. and Computat. Techn. 10 (4), pp. 277–289. Cited by: §2, §4.2.
  • E. Loli Piccolomini, V.L. Coli, E. Morotti, and L. Zanni (2018) Reconstruction of 3d x-ray ct images from reduced sampling by a scaled gradient projection algorithm. Comp. Opt. Appl. 71, pp. 171–191. External Links: Document Cited by: §2, §4.1.
  • Y. Lu, H. Chan, J. Wei, and L. M. Hadjiiski (2010) Selective-diffusion regularization for enhancement of microcalcifications in digital breast tomosynthesis reconstruction. Medical Physics 37 (11), pp. 6003. External Links: Document, ISSN 00942405, Link Cited by: §2.
  • X. Luo, W. Yu, and C. Wang (2017) An image reconstruction method based on total variation and wavelet tight frame for limited-angle ct. IEEE Access PP, pp. 1–1. External Links: Document Cited by: §1.
  • K. Niinimäki, M. Lassas, K. Hämäläinen, A. Kallonen, V. Kolehmainen, E. Niemi, and S. Siltanen (2016) Multiresolution parameter choice method for total variation regularized tomography. SIAM journal on imaging sciences 9 (3), pp. 938–974. External Links: Document, ISSN 1936-4954 Cited by: §2.
  • J. C. Park, B. Song, J. S. Kim, S. H. Park, H. K. Kim, Z. Liu, T. S. Suh, and W. Y. Song (2012) Fast compressed sensing-based cbct reconstruction using barzilai-borwein formulation for application to on-line igrt. Medical physics 39 (3), pp. 1207–1217. Cited by: §2.
  • J. Ramirez-Giraldo, J. Trzasko, S. Leng, L. Yu, A. Manduca, and C. H. McCollough (2011) Nonconvex prior image constrained compressed sensing (ncpiccs): theory and simulations on perfusion ct. Medical physics 38 (4), pp. 2157–2167. Cited by: §2.
  • I. Reiser, J. Bian, R. M. Nishikawa, E. Y. Sidky, and X. Pan (2009) Comparison of reconstruction algorithms for digital breast tomosynthesis. arXiv:0908.2610. Cited by: §1.
  • L. Ritschl, F. Bergner, C. Fleischmann, and M. Kachelrieß (2011) Improved total variation-based CT image reconstruction applied to clinical data. Physics in Medicine and Biology 56 (6), pp. 1545–1561. External Links: Document, Link Cited by: §1.
  • S. Rose, M. S. Andersen, E. Y. Sidky, and X. Pan (2015) Noise properties of CT images reconstructed by use of constrained total-variation, data-discrepancy minimization. Medical Physics 42 (5), pp. 2690–2698. External Links: Document, ISSN 0094-2405, Link Cited by: §2.
  • L. I. Rudin, S. Osher, and E. Fatemi (1992) Nonlinear total variation based noise removal algorithms. Physica D: Nonlinear Phenomena 60 (1), pp. 259 – 268. External Links: ISSN 0167-2789, Document, Link Cited by: §4.2.
  • E. Sidky, R. Chartrand, J. Boone, and X. Pan (2014) Constrained TpV-minimization for enhanced exploitation of gradient sparsity: application to CT image reconstruction. IEEE Journal of Translational Engineering in Health and Medicine 2 (1800418). External Links: Document, Link Cited by: §2, §4.3.
  • E. Y. Sidky, C. M. Kao, and X. Pan (2009a) Accurate image reconstruction from few-views and limited-angle data in divergent-beam CT. J. Xray Sci. Technol. 14 (2), pp. 119–139. Cited by: §1.
  • E. Y. Sidky, C. Kao, and X. Pan (2006) Effect of the data constraint on few-view, fan-beam CT image reconstruction by TV minimization. 2006 IEEE Nuclear Science Symposium Conference Record, pp. 2296–2298. External Links: Document Cited by: §2.
  • E. Y. Sidky, X. Pan, I. S. Reiser, R. M. Nishikawa, R. H. Moore, and D. B. Kopans (2009b) Enhanced imaging of microcalcifications in digital breast tomosynthesis through improved image-reconstruction algorithms. Medical Physics 36 (11), pp. 4920. External Links: Document, arXiv:0904.4495v1, ISSN 00942405, Link Cited by: §2.
  • E. Y. Sidky and X. Pan (2008) Image reconstruction in circular cone-beam computed tomography by constrained, total-variation minimization.. Physics in Medicine and Biology 53 (17), pp. 4777–807. External Links: Document, ISSN 0031-9155, Link Cited by: §2, §2.
  • E. Y. Sidky, I. S. Reiser, R. Nishikawa, and X. Pan (2008) Practical iterative image reconstruction in digital breast tomosynthesis by non-convex TpV optimization. Proceedings of SPIE Medical Imaging, pp. 2–7. External Links: Link Cited by: §2.
  • C. R. Vogel (2002) Computational methods for inverse problems. SIAM, Philadelphia, PA, USA. Cited by: §3, §4.2.
  • T. Wu, R. H. Moore, E. a. Rafferty, and D. B. Kopans (2004) A comparison of reconstruction algorithms for breast tomosynthesis. Medical Physics 31 (9), pp. 2636. External Links: Document, ISSN 00942405 Cited by: §1, §2.
  • T. Wu, A. Stewart, M. Stanton, T. McCauley, W. Phillips, D. B. Kopans, R. H. Moore, J. W. Eberhard, B. Opsahl-Ong, L. Niklason, et al. (2003) Tomographic mammography using a limited number of low-dose cone-beam projection images. Medical physics 30 (3), pp. 365–380. Cited by: §1.