Model-based deep CNN-regularized reconstruction for digital breast tomosynthesis with a task-based CNN image assessment approach

Assume the post-log PV at the $i$th scan angle is ${y}_{i}\in {{\mathbb{R}}}^{M}$ where $M$ is the number of detector pixels, $i=1,\,\ldots ,{\,N}_{p},$ ${N}_{p}$ is the number of scan angles, and the unknown DBT volume is $x\in {{\mathbb{R}}}^{N}$ where $N$ is the number of voxels. In MBIR, image reconstruction is formulated as an optimization problem with the following cost function:

$$\begin{eqnarray}&&\hat{x}=\mathop{{\rm{argmin}}}\limits_{x}L\left(x\right)+\beta \cdot R\left(x\right),\end{eqnarray} \tag{ 1 }$$

where $L\left(x\right)$ is the negative log-likelihood or data-fit term that measures the fidelity between the measurement and the reconstructed image, $R\left(x\right)$ is the regularization term, and $\beta$ is a regularization parameter. If the measurement noise is additive Gaussian, ${y}_{i}={A}_{i}x+{\varepsilon }_{i},$ ${\varepsilon }_{i}{\sim }{\mathscr{N}}\left(0,{K}_{i}\right),$ where ${A}_{i}\in {{\mathbb{R}}}^{M\times N}$ is the linear system matrix at the $i$th scan angle, ${K}_{i}\in {{\mathbb{R}}}^{M\times M}$ is the noise covariance matrix, then the data-fit term becomes an inverse-covariance weighted Euclidean norm:

$$\begin{eqnarray}&&L\left(x\right)=\frac{1}{2}\displaystyle \sum _{i=1}^{{N}_{p}}{\unicode{x02016}{y}_{i}-{A}_{i}x\unicode{x02016}}_{{K}_{i}^{-1}}^{2},\end{eqnarray} \tag{ 2 }$$

where ${\unicode{x02016}v\unicode{x02016}}_{W}^{2}={v}^{{\prime} }{Wv}$ and $^{\prime}$ denotes matrix transpose. Usually, ${A}_{i}$ does not model the detector blur caused by finite detector pixel size, crosstalk, or the light spread in the scintillator of an indirect detector. Zheng et al (2018) introduced a detector blur matrix $B\in {{\mathbb{R}}}^{M\times M}$ in front of ${A}_{i}$ and modeled the resulting noise correlation in the covariance matrix ${K}_{i}:$

$$\begin{eqnarray}&&{L}_{\text{DBCN}}\left(x\right)=\frac{1}{2}\displaystyle \sum _{i=1}^{{N}_{p}}{\unicode{x02016}{y}_{i}-B{A}_{i}x\unicode{x02016}}_{{\left(B{K}_{q,i}{B}^{{\prime} }+{K}_{r}\right)}^{-1}}^{2},\end{eqnarray} \tag{ 3 }$$

where ${K}_{q,i}\in {{\mathbb{R}}}^{M\times M}$ is the diagonal quantum noise matrix of the $i$th scan angle, ${K}_{r}\in {{\mathbb{R}}}^{M\times M}$ is the diagonal detector readout noise matrix. In the actual implementation, Zheng et al (2018) assumed the quantum noise variances and readout noise variances to be constant across all detector pixels, i.e. ${K}_{q,i}={\sigma }_{q,i}^{2}\cdot {\boldsymbol{I}}$ and ${K}_{r}={\sigma }_{r}^{2}\cdot {\boldsymbol{I}}$ where ${\boldsymbol{I}}$ denotes the identity matrix. Let ${{S}_{i}=\left(B{K}_{q,i}{B}^{{\prime} }+{K}_{r}\right)}^{-1/2}.$ Then ${L}_{\text{DBCN}}\left(x\right)$ has the following equivalent form using the prewhitened PV, ${\widetilde{y}}_{i},$ and the DBCN system matrix ${\widetilde{A}}_{i}:$

$$\begin{eqnarray}&&{L}_{\text{DBCN}}\left(x\right)=\frac{1}{2}\displaystyle \sum _{i=1}^{{N}_{p}}{\unicode{x02016}{\widetilde{y}}_{i}-{\widetilde{A}}_{i}x\unicode{x02016}}_{2}^{2},\,{\text{where}\,\widetilde{y}}_{i}={S}_{i}{y}_{i},{\,\widetilde{A}}_{i}={S}_{i}B{A}_{i}.\end{eqnarray} \tag{ 4 }$$

Here, ${S}_{i}$ serves as a prewhitening matrix that boosts high frequency signals but also amplifies noise. Zheng et al (2018) introduced an EP regularizer to control the noise level in DBCN reconstruction:

$$\begin{eqnarray}&&{R}_{\text{EP}}\left(x\right)=\frac{1}{1+\gamma }\displaystyle \sum _{j=1}^{N}\left(\eta \left({\left[{C}_{\to }x\right]}_{j}\right)+\eta \left({\left[{C}_{\downarrow }x\right]}_{j}\right)+\gamma \left(\eta \left({\left[{C}_{\nearrow }x\right]}_{j}\right)+\eta \left({\left[{C}_{\searrow }x\right]}_{j}\right)\right)\right),\end{eqnarray} \tag{ 5 }$$

where ${C}_{\to },{C}_{\downarrow },{C}_{\nearrow },{C}_{\searrow }{{\mathbb{\in }}{\mathbb{R}}}^{N\times N}$ are the finite differencing operators between neighboring pixels along the horizontal, vertical, and two diagonal directions in the DBT slices, $\gamma$ is an adjustable weight in the diagonal directions, and $\eta \left(t\right)={\delta }^{2}\left(\sqrt{1+{\left(t/\delta \right)}^{2}}-1\right)$ is the hyperbola potential function.

The DNGAN denoiser is a DCNN with trainable weights designed to denoise DBT images (Gao et al 2021a). DNGAN was trained with a supervised approach using pairs of noisy DBT image and target low-noise image of the same objects. In the training stage, the denoiser took pairs of image patches extracted from the two sets of DBT slices as input and target output and learned to produce the denoised image patch. The denoiser training loss was a weighted combination of the MSE loss and adversarial loss. The adversarial loss was derived from the Wasserstein GAN with gradient penalty (Gulrajani et al 2017) where the denoiser acted as the generator. We implemented the discriminator in the GAN as a trainable VGG-Net (Simonyan and Zisserman 2015). After training, the DNGAN denoiser can be deployed to full DBT slices of any size since it is fully convolutional. The performance of DNGAN denoiser was validated with both phantom and human subject DBT (Gao et al 2021a).

We integrated the DNGAN denoiser into the iterative reconstruction loop in the general reconstruction optimization problem (1). RED defines a regularizer based on a general image filter or denoiser (Romano et al 2017), which in this work is the trained DNGAN denoiser denoted as $G:$

$$\begin{eqnarray}&&{R}_{\text{RED}}\left(x\right)=\frac{1}{2}{x}^{{\prime} }\left(x-G\left(x\right)\right).\end{eqnarray} \tag{ 6 }$$

This regularizer promotes the cross-correlation between the residual after denoising and the image, or the residual itself, to be small.

By combining the DBCN data-fit term, the EP regularizer, and the RED regularizer with DNGAN denoiser, we formulated the overall optimization problem for the proposed model-based DCNN-regularized reconstruction (MDR):

$$\begin{eqnarray}&&\hat{x}=\mathop{{\rm{argmin}}}\limits_{x}{L}_{\text{DBCN}}\left(x\right)+{\beta }_{\text{EP}}\cdot {R}_{\text{EP}}\left(x\right)+{\beta }_{\text{RED}}\cdot {R}_{\text{RED}}\left(x\right),\end{eqnarray} \tag{ 7 }$$

where ${\beta }_{\text{EP}}$ and ${\beta }_{\text{RED}}$ are the regularization parameters.

A variety of optimization algorithms exist for inverse problems with ${R}_{\text{RED}}\left(x\right)$ (Romano et al 2017, Reehorst and Schniter 2019). We used the RED proximal gradient method (Reehorst and Schniter 2019) to solve (7). It introduces a secondary image variable $z\in {{\mathbb{R}}}^{N}$ and updates the primary image variable $x$ and the secondary image variable $z$ alternately:

$$\begin{eqnarray}&&{x}_{n}={\mathop{{\rm{argmin}}}\limits_{x}{L}_{\text{DBCN}}\left(x\right)+{\beta }_{\text{EP}}\cdot {R}_{\text{EP}}\left(x\right)+\frac{{\beta }_{\text{RED}}}{2}\unicode{x02016}x-{z}_{n-1}\unicode{x02016}}^{2}\end{eqnarray} \tag{ 8 }$$

$$\begin{eqnarray}&&{z}_{n}=G\left({x}_{n}\right),\end{eqnarray} \tag{ 9 }$$

where $n=1,\,\ldots ,{\,N}_{\text{iter}}$ is the iteration index. Both ${L}_{\text{DBCN}}\left(x\right)$ and ${R}_{\text{EP}}\left(x\right)$ are convex and differentiable in $x.$ We used the diagonally preconditioned gradient descent with ordered subsets for the inner minimization problem (8):

$$\begin{eqnarray}&&{x}_{n}^{i+1}={x}_{n}^{i}-\alpha \cdot P\left({N}_{p}\cdot {\widetilde{A}}_{i}^{\prime} \left({\widetilde{A}}_{i}{x}_{n}^{i}-{\widetilde{y}}_{i}\right)+{\beta }_{\text{EP}}\cdot {\rm{\nabla }}{R}_{\text{EP}}\left({x}_{n}^{i}\right)+{\beta }_{\text{RED}}\cdot \left({x}_{n}^{i}-{z}_{n-1}\right)\right)\end{eqnarray} \tag{ 10 }$$

where $\alpha$ is the step size, $n$ is incremented after $i$ goes through ${N}_{\text{inner}}$ cycles of $1,\,\ldots ,\,{N}_{p},$ ${x}_{n}^{0}={x}_{n-1}^{{N}_{\text{inner}}{N}_{p}}.$ We used the following preconditioning matrix $P$ whose inverse majorizes the Hessian of the cost function (8):

$$\begin{eqnarray}&&P={\left(\left(\displaystyle \sum _{i=1}^{{N}_{p}}\frac{1}{{\sigma }_{q,i}^{2}+{\sigma }_{r}^{2}}\text{diag}\left\{{A}_{i}^{{\prime} }{A}_{i}{{\bf{1}}}_{N}\right\}\right)+8{\beta }_{\text{EP}}\cdot {\boldsymbol{I}}+{\beta }_{\text{RED}}\cdot {\boldsymbol{I}}\right)}^{-1},\end{eqnarray} \tag{ 11 }$$

where ${{\bf{1}}}_{N}$ denotes the vector of ones of length $N.$

Clustered MCs are one of the important signs of early breast cancer and image noise can negatively impact its diagnosis. The detectability of MC clusters and the image noise are therefore important indicators of image quality for DBT reconstruction. We developed an MC classifier as a model observer for the detection task of differentiating breast structured background with and without clustered MCs in human subject DBTs for image quality evaluation. Importantly, reconstruction or other image processing processes can inadvertently enhance artifacts mimicking calcifications in the normal tissue background. Image quality assessment by MC classifiers takes into account the false positives (FPs) while other measures such as the detectability index (d') focus only on the visibility of the target objects. The MC classifier was implemented as a DCNN with trainable weights, and therefore called CNN-MC. It took an image region of interest (ROI) or patch as input. The training set included MC patches as positives with label of 1 and MC-free background patches as negatives with label of 0. The training loss was the binary cross entropy (BCE) loss, which has been shown to give the maximum likelihood estimation of the DCNN weights from the training data (Kupinski et al 2001). The CNN-MC output a score between 0 and 1 indicating the likelihood that an input patch contained clustered MCs. After training, the CNN-MC model with frozen weights was applied to the test patches.

We trained the CNN-MC for each image condition and tested it accordingly. The area under the receiver operating characteristic (ROC) curve (AUC) of the test set scores was used as an MC detectability measure. If the underlying image condition enhanced the MCs effectively, the CNN-MC should learn the MC features more accurately during training and have better classification ability between image patches with and without MCs during testing, resulting in a higher AUC metric. Any MC-like artifacts in the background patches would degrade the classification performance of CNN-MC and reduce the AUC. We studied the use of CNN-MC to rank the relative performances of different DBT reconstruction and denoising methods.

There are variations in the CNN-MC observer modeling. DCNN classifiers with different network structures learn image features differently. Even for the same network structure, the randomness in DCNN training such as kernel initialization or data batching can lead to different local minima. These are akin to the interobserver and intraobserver variabilities of human observers. To account for these variations, we investigated the VGG-Net (Simonyan and Zisserman 2015), the ResNet (He et al 2016) and the ConvNeXt (Liu et al 2022) of similar sizes, as shown in table 1, as the backbone structures of CNN-MC. We also repeated the training multiple times with different random initialization for each structure. The image condition rankings from the individual models and the final combined rankings were analyzed.

Table 1. Network structures of CNN-MC and CNN-NE in this study.

* The input to the models is a 128 × 128 pixel image patch. BN: batch normalization. LN: layer normalization. Convolutional layers are specified by the kernel size and number of filters. 'd' denotes the depth-wise convolution. 's2' denotes the stride-2 operation. The brackets denote ResNet block or ConvNeXt block.

It is important to assess the image noise for an image processing technique because enhancing the signals is always associated with a change in noise. We developed a DCNN noise estimator (CNN-NE) to quantify the root-mean-square (RMS) noise level of its input image patch. We designed their network structures to be the same as those for CNN-MC (table 1), except for without the exit softmax function, to facilitate transfer learning (discussed next). The training set contained DBT image patches of breast structured background with the training labels calculated by an ROI-based method for each patch as follows: the patch was divided into 10 × 10 pixel ROIs with 25 pixel spacing on a grid, the ROI background trend was removed by a quadratic fitting, the RMS variations of the ROI pixel values were then calculated and averaged over all ROIs as the RMS noise. The CNN-NE was implemented as a regression model. The training loss was the MSE loss between the estimated RMS noise and the training labels. After training, the CNN-NE model was applied to the test set patches. The average of the estimated RMS noise over all patches was used as a noise measure of the entire set.

For training CNN-NE and CNN-MC, we adopted the training technique of transfer learning (Yosinski et al 2014, Samala et al 2019). The CNN-NE required only breast structured background images with a range of noise levels so that a large training set could be obtained relatively easily. The trained CNN-NE then served as the pre-trained model to be further fine-tuned as CNN-MC by transfer learning. Since the availability of MC-positive samples was limited, transfer learning reduced the training sample size required and improved the model robustness. It also reused the weights that the model learned from the source task of noise estimation to encourage the CNN-MC to focus on the background noise patterns in the downstream task of MC detection.

We generated a collection of virtual breast phantoms at a voxel resolution of 0.05 mm using the anthropomorphic breast model from the Virtual Imaging Clinical Trial for Regulatory Evaluation (VICTRE) package (Badano et al 2018). The VICTRE project considered four breast density categories and one compressed thickness for each category. They stated that this design mirrored the cohort demographics of their comparative clinical trial. In our study, we used the same breast density setting as VICTRE, and slightly varied the thicknesses to make it more varied like human data. In particular, a total of 70 phantoms were simulated at a range of glandular volume fractions (GVF), including 10 almost entirely fatty (5% GVF), 10 scattered fibroglandular dense (15% GVF), 25 heterogeneously dense (34% GVF), and 25 extremely dense (60% GVF). The thicknesses of the compressed phantoms for the four density categories were 52–70 mm at every 2 mm, 46–64 mm at every 2 mm, 36–60 mm at every 1 mm, and 31–55 mm at every 1 mm, respectively.

We generated the PVs using the Monte Carlo x-ray imaging simulator (Badal et al 2021). We configured the scan geometry as 9 PVs in 25° scan angle and 3.125° increments and the x-ray spectrum as 34 kVp Rh/Ag to model the Pristina DBT system (GE Healthcare, Waukesha, WI). The distances between source and isocenter, isocenter and breast support, breast support and detector were 617 mm, 20 mm, 23 mm, respectively. Scatter and electronic noise were not simulated. The exposure was adjusted for each phantom so that the estimated mean glandular dose matched the measured value for that breast thickness from the automatic exposure control (AEC) of Pristina. To simulate a wide range of noise levels, we repeated the simulations and varied the exposure by factors of 0.25, 0.3, 0.35, 0.4, 0.5, 0.6, 0.75, 1, 1.4, 2, 3. These exposure factors were selected such that the noise standard deviations of the post-log x-ray intensity on the PVs, and hence the reconstructed DBT images, were evenly spaced. Finally, we reconstructed the DBT images using 3 iterations of the simultaneous algebraic reconstruction technique (SART) (Zhang et al 2006).

For DNGAN training, we used the virtual phantom DBTs with an exposure factor of 3 as the training high dose (HD) targets and the DBTs with the remaining 10 exposure levels as the training low dose (LD) inputs. To ensure that DNGAN covered a wide range of DBT noise levels due to variations in patient sizes and other factors in clinical settings, we trained a suite of DNGAN denoisers, each of which can handle a certain range of image noise levels. To do so, we calculated the average RMS noise for each DBT volume using the method described in section 2.2.2 and obtained the range of the noise values for the entire set of virtual DBT volumes. We empirically grouped the volumes into 12 groups so that the noise interval of each group was reasonably small. We then extracted training patches within the breast regions from each group. This resulted in 12 groups of DNGAN training data with an average of 449 245 (min: 343 498; max: 547 595) LD/HD pairs of 32 × 32 pixel training patches for training a denoiser for each group.

We randomly extracted 128 × 128 pixel patches within the breast regions from all the 70 virtual phantom DBTs and 11 exposure levels to form the CNN-NE training set. This gave a total of 256 194 patches. The low-frequency background was removed from the patches to reduce the nonuniformity caused by heterogeneous breast structures (Chan et al 1995).

We had 238 human subject cases with biopsy-proven MCs collected with Institutional Review Board approval and written informed consent. Two-viewed DBTs of the breast with MCs were acquired for each case using the GEN2 prototype DBT system (GE Global Research, Niskayuna, NY). The system acquired 21 PVs in 60° scan angle and 3° increments with a 29 kVp Rh/Rh x-ray beam. The distances between source and isocenter and between isocenter and detector were 640 mm and 20 mm, respectively. The breast support was at the same height as the isocenter. For our reconstructions, we used the central 9 PVs in 24° angle so that the scan geometry was similar to that of the GE commercial Pristina DBT system. The radiation doses at 9 PVs were similar to those of the Pristina system. A Mammography Quality Standards Act (MQSA) approved radiologist marked the biopsied MCs with 3D bounding boxes in the reconstructed DBT images based on all available clinical information.

We split the human subject DBTs by case into three disjoint sets for training, validation and testing. The training set consisted of 64 cases, or 127 DBT views (one view was lost due to technical issue). The validation set consisted of 52 cases, or 104 DBT views. The remaining 122 cases with 246 views (one case had bilateral MCs) were sequestered as independent test set. For each data set, positive patches of size 128 × 128 × 3 pixels containing clustered MCs were extracted inside the radiologist's 3D boxes on a regular grid with centers separated by 128 pixels in the two directions along a slice and separated by 1 slice with offset centers along the depth direction. An MC cluster therefore might be extracted multiple times in part or in whole but spatially shifted in each patch, which served as augmented samples to reduce the imbalance between the two classes. Non-overlapping negative patches of the same size were randomly extracted outside the 3D boxes in the structured breast background. For the CNN-MC, we obtained 751 positive and 19 079 negative patches for the training set, 709 positive and 1955 negative patches for the validation set, and 2289 positive and 6176 negative patches for the test set. For each patch, we removed the low-frequency background and took maximum intensity projection (MIP) over the three slices to emphasize the MCs, if any. Figure 1 shows example MIP patches from the CNN-MC training set. For the validation and test sets of CNN-NE, the same sets of negative patches for CNN-MC were extracted except that only the central slices without MIP were used.

Zoom In Zoom Out Reset image size

The physical phantom DBTs were used in a previous observer study with radiologists to detect MC clusters under three image conditions (Chan et al 2023). Each of the six breast phantoms were 5 cm thick and made of 50% glandular/50% adipose heterogeneous breast-tissue-equivalent materials. Clusters of simulated MCs of four nominal speck diameters (0.125–0.150 mm, 0.150–0.180 mm, 0.180–0.212 mm, 0.212–0.250 mm) were embedded in the phantoms. There were 36 clusters for each diameter range, giving a total of 144 clusters. The phantoms were imaged twice using two automatic exposure modes of the Pristina DBT system: the standard (STD) mode, which was for routine patient imaging, and the STD+ mode, which used about 54% more dose than that of STD. The DBT images were reconstructed by the built-in commercial algorithm. A third image set called dnSTD was obtained by denoising the STD set using the DNGAN denoiser. The denoiser was trained and validated by Gao et al (2021a) using a separate set of physical phantoms that shared the same materials but had different designs from the test phantoms in the Chan et al (2023) study. For each image set and each MC speck size, we extracted a total of 36 128 × 128 pixel background-corrected 3-slice MIP patches centered at the MC clusters as positives and paired them with 75 MC-free background MIP patches as negatives for CNN-MC testing. The set of 75 negative patches was extracted from different random locations for each MC positive set.