{'fraction_non_alphanumeric': 0.029030329807417418, 'fraction_numerical': 0.01114951789380669, 'mean_word_length': 4.744699910421021, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 3, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Optimization methods for very accurate Digital Breast Tomosynthesis image reconstruction

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.


Introduction
Digital Breast Tomosynthesis (DBT) is a 3D Xray cone-beam Computed Tomography (CT) technique for the early detection of breast tumors [1,2].
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 [3], produce poor noisy images in limited-angle tomography, hence they have been left in favour of Iterative Reconstruction (IR) algorithms [4,5,6]. 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 [7].
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 [8,9,10,11]. 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 [12,8,9,10,13,11]. 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 Z 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.

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 [7]. As a consequence, a wide amount of IR methods has been proposed to reconstruct tomographic images and an exhaustive analysis can be found in [13].
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 [14].
On the other hand, the worldwide increasing interest in Compressive Sensing (CS) [15] 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 [12,13] 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 [4,5,16,17,18,19,20,21,22,23]. 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 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 k-th iteration of the algorithm computing the sequence {x (k) } k 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. 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 [17,24,25]. 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 ap-plications [26]. In our previous works we have investigated a Fixed Point (FP) algorithm and a Scaled Gradient Projection (SGP) method in [27] and [28] respectively, and we have applied them to small simulated data sets. A parallel SGP implementation on GPUs has been presented in [29], where we were interested in showing the computational efficiency of the algorithm in terms of execution time. We also remark that in [30] 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 [31] 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 [32] 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 [32].

The optimization framework in modelbased 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 Z-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 Z vertical direction. In order to define the numerical model of tomographic image formation, we discretize the 3D object into N v voxels, whereas the 2D detector panel is made of N p recording units. For each fixed projection angle θ and i-th detector recording unit, the Lambert-Beer law relates the projections P θ i , along a ray R θ , to the attenuation coefficient function µ(w) of the voxel w crossed by R θ [33] as: where P 0 represents the intensity of the energy emitted by the X-ray source. The discretization of the integral in (1) for all the N θ scanning angles arises the following linear system: In equation (2) we denote with x the N v dimensional vector stacking the attenuation coefficients of all the voxels, while b is the vector of size N d = N p × N θ storing all the projections (i.e. the right hand sides of (1)) and M is the matrix of size N v × N d , 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 N v > N d ) 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 [13]. We consider here the former problem and express it as: where J(x) is a fit-to-data function, R(x) 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 x ≥ 0 reflecting the nonnegativity property of the linear attenuation coefficient µ in (1).
In particular, in this work we settle J(x) as the Least Squares (LS) function and R(x) as the Total Variation (TV) operator defined as [34]: 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 [34]. Exploiting the linearity of (3), the objective function gradient ∇f (x) = ∇LS(x) + λ∇T V (x) can be evaluated by separately computing ∇LS(x) as and ∇T V (x) through finite forward differences.

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 modelbased 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 matrixvector products), to efficiently run the solver in short time; • a limited request of memory, to solve realsize problems on commercially affordable hardware.
A challenging issue (common to the implementation of the three algorithms) is the computation of the projection matrix M : 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 M .

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 k-th iteration, the new solution is computed by moving along a descent direction d (k) of a quantity η k > 0, as: The direction d (k) is obtained through a projection P + onto the non-negative orthant: where α k is the step length and S k is the scaling matrix (step 7 in Algorithm 1). Essentially, the method follows a Gradient Projection approach accelerated by choosing the α k step length with Barzilai-Borwein techniques and by introducing a suitable scaling matrix improving the matrix conditioning [28]. In particular, the scaling matrix S k is a diagonal matrix with entries in a limited interval. To update S k (line 5 of Algorithm 1), we compute a splitting of the objective function gradient into its positive and negative parts, as: where where {ρ k } k is a decreasing positive sequence.
Regarding the convergence, it is proved in [35] that the SGP algorithm converges without any further restriction on the step length α k and on the scaling matrix S k to the unique minimum of (8).
In [36], the authors proved that the theoretical convergence rate of the SGP method is O(1/k).

The Fixed Point algorithm
The FP algorithm for the solution of the minimization problem: arg min has been firstly proposed for image denoising by Rudin, Osher and Fatemi in [37]. Starting from this approach we derived the lagged diffusivity FP Algorithm 2 for 3D tomographic image reconstruction.
Algorithm 2 Lagged diffusivity Fixed Point algorithm (FP) , with the Conjugate Gradient method. 5: At each k-th iteration, the FP algorithm updates the solution with the following rule: where the descent direction d (k) is computed by solving a linear system ) in line 4 approximates the Hessian matrix. It contains in fact the seven diagonals banded matrix L(x (k) ) which is the discretization matrix of the diffusion operator L(x) so that L(x)x = ∇T V (x) [34]. We solve the linear system with very few iterations of a Conjugate Gradient (CG) algorithm [38]: 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 H k and that, to save memory space, we perform it without storing the matrix H k : we only store L(x (k) ) and re-compute M and M T 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 [27] and [39] respectively.

The Chambolle-Pock algorithm
The CP algorithm has been firstly proposed in [40] for the solution of the general minimization problem: arg min where G is a convex, lower-semicontinuous and proper function, F is convex and lowersemicontinuous and K is a continuous linear operator. To fit the problem statement (15) and fully exploit the linearity of K, we assign: which results in defining the operator K as a matrix composed by the four following blocks: where ∇ x , ∇ y and ∇ z are the forward differences operators acting along the X, Y and Z axes respectively. In order to include the non-negative constraints, we fixed G(x) = δ Ω (x) as the indicator function of the convex set Ω = {x : x ≥ 0}, i.e. δ Considering the convex conjugate F * of F , defined as F * (y) = max x {x T y − F (x)}, and the proximal mappings of G and F * , i.e.
the k-th CP iteration can be described with the following three steps: with an extrapolation step: ) and θ > 0.
In particular, the proximal mapping prox σ [F * ] 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 [31]. The proximal mapping of G is defined as: hence it is exactly the projection P + (x) of x onto the feasible set Ω (lines 9-10 of Algorithm 3). The updated iterate x (k+1) is computed with a FISTA strategy as in line 11 of Algorithm 3. The algorithm convergence is demonstrated in [40]. We finally observe that the algorithm needs to compute the value Γ (line 1 of Algorithm 3): to estimate the matrix 2-norm as Γ ≈ K 2 = ρ(K T K) (where ρ is the spectral radius of a matrix), we perform two iterations of the power method for the maximum eigenvalue computation [41].

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 {λ k } k , k = 1, . . . of regularization parameters defining a sequence of minimization problems stated as (13), whose solutions converge to a good reconstructed image. See [42] for more details and the convergence proof.
We propose the following fully-automatic strategy to compute a decreasing sequence {λ k } k . At the beginning of our algorithm, we leave out the regularization by setting the first parameter λ 0 = 0: we are in fact interested in a very good data fitting, to recover as many image features as possible. Next, the starting value λ 1 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 1/k at each k-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 λ 0 = 0 to initialize the algorithm and run the first iteration (labelled with k=0) to compute x (1) ; (1) ) and use it to compute x (2) ; • For each k ≥ 2, set and use it to compute x (k+1) .

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 M 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 M T . Different algorithms have been proposed in literature for the computation of the matrix M . We have adopted the Distance Driven (DD), which accurately models the discretization of the Lambert-Beer's law (1) for cone-beam projections [43]. In DD, M of size N v × (N p × N θ ) is constituted by N θ submatrices M θ of size N v × N p . Each element M θ i,j represents the contribution of the j-th voxel (for j = 1, . . . N v ) to the projection onto the i-th detector pixel (for i = 1, . . . N p ), 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 i-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 j-th voxel and the whole magenta extension is proportional to the value M i,j of the matrix. For all the voxels j not contributing to the i-th projection the corresponding matrix element M i,j = 0; hence M is extremely sparse. However, despite the huge number of nonzero elements, for its very large size, in real applications M 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 [43] 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.

DBT system configuration
Our tests are performed on the digital system Giotto Class of the Italian I.M.S. Giotto Spa company in Bologna [44]. The source executes N θ = 11 scans from equally spaced angles in approximately 30 degrees range; in the highest vertical position, the source is about 70cm over the detector. The stationary digital detector has a sensitive area of 24 cm × 30 cm and squared pixel pitch of 0.085 mm; the reconstructed voxel dimensions along the three cartesian axes are ∆ x = ∆ y = 0.090 mm and ∆ z = 1 mm 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.

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 [45]. 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 10 cm × 18 cm. Inside one of them, we find acrylic spheres simulating breast masses (MSs), 1 cm 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 11.4 cm × 21 cm 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 µm diameter, respectively) and the second and fourth mass (with diameter 4.7 and 3.1 mm, 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 10.5 cm × 20.7 cm.  Table 1: Diameters of a microcalcification (MC) in a cluster and of a mass (MS) in the BR3D phantom as reported in [45]. Measures are in micrometers (µm).

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 M M C 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 XY plane) where the microcalcification lies and we extract the Plane Profile (PP) along the Y axes. The FWMH index is computed as: where d 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 Z 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: is the mean of the reconstructed values inside a circular region of three pixels diameter inside the considered MC and in the background,z corresponds to the slice where the object is on focus and N z is the total number of discrete slices. Similarly, we compute the ASF for the masses.

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.

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 k in Algorithm 2), respectively.
The value of β in (6) has been fixed as β = 0.001. 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 λ = 0.005 for both SGP and CP methods and λ = 0.001 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 125 × 125 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.

SGP algorithm insights
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 Z 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 M 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 w 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 iter- ations (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.

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     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.

Experiments with a variable regularization pa-
rameter 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  Table 4: Values of CNR and FWHM measures on the mass and the microcalcification observable in Figure 8 (a)-(c).
λ. In this paragraph we show the results achieved with the automatic rule (21) for the choice of a decreasing sequence {λ k } k applied to the SGP solver, to reconstruct the BR3D phantom. Figure 10 plots the sequence {λ k } k 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 λ = 0.005 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.

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 µm 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 +150% to +280%. 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.

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