Deep learning reconstruction of digital breast tomosynthesis images for accurate breast density and patient-specific radiation dose estimation
Jonas Teuwen, Nikita Moriakov, Christian Fedon, Marco Caballo, Ingrid Reiser, Pedrag Bakic, Eloy García, Oliver Diaz, Koen Michielsen, Ioannis Sechopoulos
DDeep learning reconstruction of digital breast tomosynthesis images for accuratebreast density and patient-specific radiation dose estimation
Jonas Teuwen a,b,1 , Nikita Moriakov a,b,1 , Christian Fedon a , Marco Caballo a , Ingrid Reiser e , Pedrag Bakic f,g , EloyGarc´ıa c , Oliver Diaz d , Koen Michielsen a , Ioannis Sechopoulos a,h a Department of Medical Imaging, Radboud University Medical Center, the Netherlands b Department of Radiation Oncology, Netherlands Cancer Institute, the Netherlands c Vall d’Hebron Institute of Oncology, VHIO, Spain d Department of Mathematics and Computer Science, University of Barcelona, Spain e Department of Radiology, The University of Chicago, USA f Department of Radiology, University of Pennsylvania, USA g Department of Translational Medicine, Lund University, Sweden h Dutch Expert Centre for Screening (LRCB), the Netherlands
Abstract
The two-dimensional nature of mammography makes estimation of the overall breast density challenging, and es-timation of the true patient-specific radiation dose impossible. Digital breast tomosynthesis (DBT), a pseudo-3Dtechnique, is now commonly used in breast cancer screening and diagnostics. Still, the severely limited 3rd dimensioninformation in DBT has not been used, until now, to estimate the true breast density or the patient-specific dose. Inthis study, we propose a reconstruction algorithm for DBT based on deep learning specifically optimized for thesetasks. The algorithm, which we name DBToR, is based on unrolling a proximal primal-dual optimization method,where the proximal operators are replaced with convolutional neural networks and prior knowledge is included in themodel. This extends previous work on a deep learning based reconstruction model by providing both the primal andthe dual blocks with breast thickness information, which is available in DBT. Training and testing of the model wereperformed using virtual patient phantoms from two di ff erent sources. Reconstruction performance, as well as accu-racy in estimation of breast density and radiation dose, was estimated, showing high accuracy (density < ± < ± Keywords: digital breast tomosynthesis, reconstruction, deep learning
1. Introduction
Breast cancer screening with mammography hasproven e ff ective in reducing breast cancer related mortal-ity (Plevritis et al., 2018). However, since mammography ∗ Corresponding author: Ioannis Sechopoulos
Email address: [email protected] (Ioannis Sechopoulos) Equal contribution. is a 2D imaging modality, it results in the projection of theinternal tissue of the breast onto a single plane, yieldingtissue superposition. This results in a lower sensitivity andspecificity, especially with dense breasts. In addition, thislack of information on the true tissue distribution in thevertical (x-ray source to detector) direction limits the abil-ity to accurately estimate the breast density (i.e., the pro-portion of the breast that consists of fibroglandular tissue).Up to now, estimates of breast density from mammogra-phy necessitated the use of models of the image acquisi-
Preprint June 12, 2020 a r X i v : . [ phy s i c s . m e d - ph ] J un ion process and assumptions and simplifications regard-ing tissue distribution. Although these models have beencharacterized for consistency and precision, there accu-racy is unknown. Estimation of breast density is of intenseinterest due to it being an important factor for both mask-ing e ff ects and breast cancer risk (McCormack, 2006). Asa result of its significance, breast density reporting is nowmandated by law in several states in the USA. This makesit especially important that we use methods that providereliable, robust, and consistent breast density estimates.In addition to limiting the accuracy of breast densityestimation, the 2D nature of mammography makes it im-possible to estimate the radiation dose received by a spe-cific patient during image acquisition. This is because,in breast imaging, the dose of interest is only that to thefibroglandular tissue, since this is the tissue at risk of de-veloping breast cancer (Dance and Sechopoulos, 2016),and to determine this dose it is necessary to know its ver-tical location within the breast. Therefore, currently allbreast dosimetry is based on dose estimates using a modelbreast, which does not reflect the dose deposited in theactual patient breast. It has been shown that the use of amodel breast results in an average over-estimation of thetrue patient breast of 30% and that this error can be ashigh as 120%, if not more (Dance et al., 2005; Sechopou-los et al., 2012; Hernandez et al., 2015). Even if a moreaccurate, unbiased model of the average breast is devel-oped, an e ff ort which is currently ongoing (Arana Pe˜naet al., 2020), the over 100% error in patient-specific doseestimates by use of a population-wide model will not beameliorated. To obtain accurate radiation dose estimates,the actual amount and position of the fibroglandular tis-sue in the actual patient’s breast needs to be considered.In this work we address both the estimation of amount andlocation of the fibroglandular tissue within each patient’sbreast, resulting in accurate estimates of both the breastdensity and fibroglandular dose.Digital breast tomosynthesis (DBT) has been intro-duced over the last two decades to decrease the impact oftissue superposition in mammography, by providing lim-ited depth information, resulting in improved detectionand diagnosis performance (Zackrisson et al., 2018; Zu-ley et al., 2013). Hence, DBT is rapidly replacing digitalmammography as the primary x-ray technique for breastimaging (Niklason et al., 1997). DBT imaging consists ofacquiring several low-dose planar x-ray projections over a limited angle. These projections are then used to re-construct a pseudo-3D volume, albeit with very limitedvertical spatial resolution. Even so, the introduction of animaging modality that provides some vertical informationinto widespread use for screening provides the opportu-nity, for the first time, to obtain accurate estimates of thebreast density and fibroglandular dose. However, giventhe very poor spatial resolution in the vertical directionobtained with all current DBT reconstruction algorithms,the localization of the fibroglandular tissue in DBT im-ages has proven extremely challenging. Previous e ff ortsto achieve this have not resulted in improvements overthe model-based dose estimates obtained from 2D imag-ing (Geeraert, 2014), or have involved algorithms that re-quire a lot of manual input, making them challenging toimplement clinically (F¨ornvik et al., 2018).In general, reconstruction of a 3D volume from com-puted tomography (CT) or DBT projections is an inverseproblem, and there exist multiple algorithms for solvingsuch problems. For DBT, in particular, the choice ofthe reconstruction algorithm and regularization parame-ters can greatly influence the reconstruction quality, aswas shown in previous work (Rodriguez-Ruiz et al., 2017;Michielsen et al., 2016). In this work, we propose anew approach to DBT reconstruction, based on our earlierwork (Moriakov et al., 2019), using deep learning meth-ods, that results in the 3D representation of the imagedbreast optimized for estimation of the true distribution ofthe fibroglandular tissue. This in turn allows for accu-rate estimation of both the breast density and the radia-tion dose imparted on it, improving significantly upon thestate-of-the-art.A recent and very promising development in medicalimaging reconstruction is using methods that rely on deeplearning (Arridge et al., 2019). The goal of this paperis to show the potential of such methods for the problemof DBT reconstruction, using the quantitative estimationof breast density and radiation dose as the target applica-tion. For this, we extended our previous results on DBTreconstruction (Moriakov et al., 2019) and investigate adata-driven reconstruction algorithm called Deep BreastTomographic Reconstruction (DBToR), where the recon-struction is computed from projection data with a deepneural network. To train the model, and to test the per-formance of DBToR for the tasks of breast density andradiation dose estimation, we used dedicated breast CT2mages, where the tissue distribution is known, and use afinite-element model to simulate the change in the tissuedistribution when under compression in DBT. In additionto these images, we also evaluate the model on simulatedbreast phantoms. In this work, we show: (i) the feasibility of recon-struction of DBT using deep learning, having developeda novel deep learning-based architecture to reconstructDBT; (ii) that the resulting reconstructions have greatlyimproved vertical resolution compared to state-of-the-artanalytical and iterative reconstruction methods; (iii) thatthe proposed reconstruction method is able to provide ac-curate breast density estimates; and (iv) that the densetissue location information results in accurate patient-specific dosimetric estimates.
2. Deep learning-based reconstruction
Before we describe our model in more detail, we give abrief overview of reconstruction problems formulated asinverse problems.
Image reconstruction can be stated as an inverse prob-lem. Formally this means that given an image x ∈ X andmeasured projection data y ∈ Y , we have y = Ax + η (1)where A : X → Y is the forward operator, or projec-tion operator, that models how the image x gives rise tothe projection Ax in the absence of noise, and η is a Y -valued random variable modeling the noise componentof the measurements. Typically, the spaces X and Y areBanach spaces, and here X and Y denote the spaces offunctions describing true breast anatomy and measure-ments, respectively. The noise vector η can typically beconsidered to be Gaussian when the photon count is largeenough, as is typically the case in transmission imaging.This problem can also be considered from a Bayesianperspective, where X and Y are probability spaces and the probability distribution x ∼ P ( x ) on X is called the prior .Bayes theorem states that P ( x | y ) = P ( y | x ) P ( x ) P ( y ) , (2)where the conditional probability P ( y | x ) is called the like-lihood of data y , which expresses the probability of takingmeasurements y from given data x . The likelihood is de-rived from the forward model, using the knowledge of theforward operator A and the noise distribution η .The goal of reconstruction is to retrieve the image x from the noise-corrupted measurements y . Inversion ofthe operator A is generally an ill-posed problem, and thereexist di ff erent approaches towards estimating x . From aprobabilistic perspective, two of the most important es-timators are Maximum Likelihood and
Bayes estimators.In the
Maximum Likelihood estimator, the estimated x isgiven by finding x ML that maximizes the likelihood of data y . For Bayes estimator, the goal is minimizing the ex-pected loss over all estimators ˆ x : Y → X that give anestimate of x given measurements y . That is,ˆ x Bayes : = argmin ˆ x : Y → X E x ∼ X ( ˆ x ( y ) − x ) . (3)It can be shown that the Bayes estimate of the unknownparameter ˆ x Bayes from measurements y is simply the meanof the posterior distribution P ( x | y ), that is,ˆ x Bayes ( y ) = (cid:90) X xP ( x | y ) dx . (4)In our approach, we will obtain a neural network ap-proximation to the Bayes estimator ˆ x Bayes , where argminin Equation (3) is taken over the estimators given by afamily of neural networks and where optimization is per-formed using minibatch stochastic gradient optimization.
The DBToR algorithm, a data-driven algorithm de-signed for DBT reconstruction, extends the LearnedPrimal-Dual (LPD) reconstruction algorithm (Adler and¨Oktem, 2018) to incorporate additional prior informationabout the geometry in the form of the thickness measure-ment of the breast under compression in DBT (as mea-sured by the DBT device and stored in the DICOM image3 h gf m + f m Outcopyconv 3x3 + PReLU conv 3x3Projection Backprojectionm - mask, g - sinograms, h - initial dual vector (zeros), f - initial primal vector (zeros), Out - final reconstruction Figure 1: Network architecture of DBToR. Dual blocks (green) are on the upper row and primal blocks (blue) are in the bottom row. The blockshave the same architecture, elaborated in the first blocks. m : the breast thickness mask, h : initial dual vector, f : initial primal vector, g : sinogramdata, Out: final reconstruction header). The algorithm provides a Bayes estimator givenby a deep neural network, which is trained to reconstructthe images directly from projection data. The choice ofa particular neural network architecture can greatly a ff ectthe results and is non-trivial. The DBToR neural networkconsists of several ‘reconstruction blocks’, which take inprojection data, together with information on the thick-ness of the breast under compression as the initial input,perform a forward and a backward pass by taking projec-tions and back-projections, and use a convolutional neuralnetwork to produce an intermediate reconstruction result,which is then improved further by each successive recon-struction block. The architecture of the network is illus-trated in Figure 1.We denote the training set of images D train . For animage x , we let sinogram x be the corresponding projec-tion data. We assume that the training data D train is arepresentative sample from the domain of images thatwe want to reconstruct (DBT images). The neural net-work is trained in a supervised fashion as follows. Werepeatedly sample image x ∼ D train and the correspond-ing input projection data y = sinogram x from the trainingdataset. The corresponding thickness mask is denoted by m = thickness mask x and is represented by a rectangu-lar mask with the same width as the detector and wherethe height is given by the measured breast thickness aftercompression. These measurements are provided by theDBT system, and are available both at training and at testtime.Using this information, the reconstruction z : = m , y (cid:55)→ compute reconstruction( m , y )is given by a deep neural network with parameters θ . Weselect the L -loss l θ : = (cid:107) x − z (cid:107) to train the network ina supervised fashion. The input projection data is log-transformed, and scaled such that the standard deviationand mean over the training set is 1. The algorithm is de-scribed in Figure 2 and the network architecture in Fig-ure 1. In all experiments we set the number of primalblocks N prim and dual blocks N dual to 10. The parame-ters θ are updated using the Adam optimizer with a cosineannealing learning rate schedule (Loshchilov and Hutter,2019), i.e. the learning rate in step t was η t = η (cid:16) + cos (cid:16) π tN iter (cid:17)(cid:17) : procedure compute reconstruction ( m , g ) f ← ∈ X N p rim (cid:46) Initialize primal vector h ← ∈ U N d ual (cid:46) Initialize dual vector for i ← , N do h i ← Γ θ di ( h i − , P ( f (2) i − ) , g , P ( m ))) f i ← Λ θ pi ( f i − , P ∗ ( h (1) i ) , m ) end for return f (1) I end procedure for j ← , N iter − do if j is divisible by freq then x ∼ D train (cid:46) Sample train data y ← sinograms x (cid:46) Sample sinograms m ← thickness mask x (cid:46) Create masks end if z ← compute reconstruction ( m , y ) loss ← (cid:107) z − x (cid:107) change parameters θ pi , θ di , i = , . . . , N to reduce loss end for Figure 2: Pseudocode of the DBToR reconstruction algorithm. starting at a learning rate η of 10 − . The other Adam pa-rameters were the default choices. N iter is the total num-ber of iterations (as in Algorithm 2). Γ θ di and Λ θ pi for i = , . . . , I are ResNet-type blocks consisting of threeconvolutional layers with kernel size 3 × P and P ∗ are the forward and back-ward operators respectively. At test time, the algorithmtakes the input projection data sinogram x and the breastthickness information to compute the reconstruction us-ing the function compute reconstruction.In total, we trained three versions of the DBToR al-gorithm: two versions were trained on virtual phantomprojections, both without any noise and with varying lev-els of noise, and one version was pretrained on noisyvirtual phantom projections and subsequently finetunedon noisy deformed Breast Computed Tomography (BCT)phantoms. This data is described in full detail in Section3.1 and Section 3.2 respectively. For comparison, we alsotrained the basic LPD algorithm in addition to DBToR onnoise-free virtual phantom projections. This was achievedby removing the height mask from the algorithm input. The final step of the reconstruction for estimation ofbreast density and radiation dose involved in image ac-quisition, is the classification of the reconstructed imageinto skin, adipose, and fibroglandular tissue voxels. Thisis required because rather than relying on voxel attenua-tion values, which can be quite similar for di ff erent typesof tissue, we need correct tissue labels to compute thebreast density and the dose absorbed by the fibroglandu-lar tissue. Given its high contrast, the skin layer was seg-mented through a fast seeded region-growing algorithm(Adams and Bischof, 1994), which grows the segmentedregion starting from a subset of seeds (corresponding tothe voxels located on the outer edge of the skin layer) bysubsequently including voxels whose intensity was higherthan or equal to the mean seed intensity value. For fi-broglandular tissue classification, we extend our previouswork on breast CT classification (Caballo et al., 2018),and in the first step remove the skin temporarily fromthe image, with the resulting representation undergoinga well-established, automatic thresholding method basedon fuzzy c-means clustering (Bezdek, 1981), an algorithmgenerally adopted in the case of images with low noisecontent, accompanied by a non-negligible degree of blur-ring, as is the case for our images. Briefly, voxels are iter-atively assigned to a given class (adipose or fibroglandulartissue) in an unsupervised fashion, with the iteration stop-ping criterion aiming at maximizing the distance betweenthe average voxel values of the two classes. As opposedto traditional cluster analysis, this method allows for a de-gree of fuzzy overlap between the classes over each iter-ation, which helps classify the boundary voxels in eachsubsequent iteration. The fuzzy partition term (Bezdek,1981) was experimentally tuned to a value of 2.0.
3. Materials and Methods
To train and evaluate the algorithm, we created twodatasets of 3D breast phantoms from which we extractedthe coronal slices and their corresponding DBT projec-tions. The first dataset consists of virtual 3D breast phan-toms generated using a stochastic model, while the sec-ond is based on patient dedicated BCT images. DBT pro-jections of these phantoms were simulated using deter-ministic simulation methods, with the posterior addition5f Poisson noise. The use of virtual phantoms not onlyprovided training data, but also allowed for assessing theaccuracy of the density and dose estimates, since groundtruth is known.Each voxel of these phantoms was indexed with a la-bel denoting the corresponding tissue type: skin, adiposetissue, fibroglandular tissue, and Cooper’s ligaments. Theelemental compositions of these materials were obtainedfrom the work of Hammerstein et al. (1979), except forthe composition of Cooper’s ligaments, which was as-sumed to be identical to that of fibroglandular tissue. Lin-ear attenuation coe ffi cients at 20 keV, a typical averageenergy of the spectra in clinical DBT imaging, were cal-culated for each material using the software from Booneand Chavez (Boone and Chavez, 1996). The resulting lin-ear attenuation coe ffi cients were 0 .
512 cm − for adiposetissue, 0 .
798 cm − for fibroglandular tissue (and Cooper’sligaments), and 0 .
854 cm − for skin.This process was done for the virtual phantom datasetand the patient BCT dataset. We extracted 41499 coronal slices from 50 breast phan-toms generated using the method of Lau et al. (2012).This method generates breast phantoms in two steps: first,the breast structure is simulated on a coarse scale by gen-erating large compartments of adipose tissue Zhang et al.(2008); Bakic et al. (2011). Second, finer detail for fi-broglandular tissue is added subsequently in the form ofpower-law noise (Reiser and Nishikawa, 2010). An exam-ple of the simulated phantoms is shown in Figure 3. These
Figure 3: 2D coronal breast phantom containing skin (darkest gray), adi-pose tissue (dark gray), fibroglandular tissue (light gray), and Cooper’sligaments (black).
Patient dedicated breast CT images were acquired foran unrelated, ethical-board approved patient study evalu-ating this imaging technology. The images were releasedfor other research purposes after anonymization. In orderto compute the density and the accumulated dose to thefibroglandular tissue, the patient breast CT images wereautomatically classified into four categories (air, skin, adi-pose and fibroglandular tissue) using a previously de-veloped algorithm for BCT image classification (Caballoet al., 2018) (see Section 2.3). The classified breasts thenunderwent simulated mechanical deformation as previ-ously described (Fedon et al., 2019; Garc´ıa et al., 2020).Briefly, the breasts were converted into a finite element(FE) biomechanical model using the package iso2mesh(v.1.8; Matlab v.13a). A large number of 4-node tetra-hedral elements, between 100k and 500k, were used inorder to minimize the numerical error during the FE anal-ysis (del Palomar et al., 2008). Nearly incompressible(Poisson ratio equal to 0.495), homogeneous and isotropicNeo-Hookean material models for each tissue were usedto describe their mechanical behaviour. The Youngs mod-ulus for fibroglandular, adipose and skin tissue were setto 4 .
46 kPa, 15 .
10 kPa, and 60 .
00 kPa, respectively (Well-man, 1999). The mechanical compression was then sim-ulated using the open-source package NiftySim (v.2.3.1;University College London) (Johnsen et al., 2015), whichuses a Total Explicity Dynamic Lagrangian approach tosolve the mechanical FE problem (Miller et al., 2007).Each breast model was compressed to the thicknessrecorded in the corresponding DICOM header of thecranio-caudal DBT view of that breast, which was ac-quired, for clinical purposes, during the same visit as theacquisition of the BCT image.The total phantom population includes compressedbreast thicknesses from 3 . . . . . × . × . ffi cient for dosimetric applications (Fedonet al., 2019).Using this method, of which an example is given inFigure 4, we obtained a total of 28891 deformed BCT6 a) (b) (c)Figure 4: (a) Coronal slice of a breast CT image, (b) the same image classified into skin (white), adipose (dark gray) and fibroglandular (light gray)tissue voxels, and (c) the classified deformed image with the technique described in Section 3.2. slices extracted from 91 di ff erent patient breasts. Giventhat the number of deformed BCT slices was substantiallylower than that of virtual phantom slices, we pre-trainedthe model with the latter, and then fine-tuned the modelusing the BCT slices from 46 patient BCT images. Theother 45 patient BCT image phantoms were used for test-ing the reconstruction performance of the model and theaccuracy of density and dosimetry estimates. Each patientbreast was either completely included or excluded whenselecting slices for fine-tuning and testing the model inorder to prevent data contamination and bias. Limited angle fan-beam projections were simulated forall coronal phantom slices using a geometry with the cen-ter of rotation placed at the center of the phantom as seenin Figure 5. The x-ray source was placed 65 cm above thecenter of rotation, and the source-detector distance was70 cm. A total of 25 equally spaced projections between − ◦ and 24 ◦ were generated, with the detector rotatingwith the x-ray source. The detector was a perfect pho-ton counting system (100% e ffi ciency) consisting of 1280elements with a resolution of 0 . y i ( x ) = b i exp (cid:16) − (cid:88) j l i j x j (cid:17) was used for the simulations, with ˆ y being the simulatedprojection data, b i the number of x-ray photons emittedtowards detector pixel i , l i j the length of the intersectionbetween voxel j and the line between the source and de-tector pixel i . The linear attenuation in voxel j is denotedby x j . For the virtual phantom data, we generated a se-ries of data sets at 3 noise levels from the noiseless simu-lated projections. This was accomplished by setting pho-ton count b i = · √ N with N = { , , } . For eachnoise level, a single Poisson noise realization was gen-erated. For the deformed BCT phantoms, only photoncounts corresponding to N = ff ects were not includedin the computation and thus were not needed in the recon-struction. The smoothing prior was not included becauseit does not influence the spatial distribution of fibroglan-dular tissue in the reconstruction.7 .4. Density computation Breast density by mass, also called glandularity ( G ),was computed as follows: G = N g ρ g N g ρ g + N a ρ a (5)where N g is the number of voxels classified as fibrog-landular tissue, N a is the number of voxels classified asadipose tissue, and ρ g and ρ a are the corresponding den-sity for fibroglandular and adipose tissue, respectively, ac-cording to Hammerstein et al. Hammerstein et al. (1979).The true breast density of each BCT phantom was ob-tained by applying this equation to the phantom volumesthemselves. The estimated breast density resulting fromthe proposed method was obtained by applying the equa-tion to the classified reconstructed images. The mean glandular dose (MGD) estimations were per-formed using a previously described and validated MonteCarlo code (Fedon et al., 2018a,b), based on the Geant4toolkit (release 10.05, December 2018). The Monte Carlogeometry replicates the one used to generate the projec-tions and is shown in Figure 5.As during simulation of the DBT projections above,each voxel was labelled with an index related to its com-position: air, adipose tissue, fibroglandular tissue, andskin, using the chemical compositions reported by Ham-merstein et. al. Hammerstein et al. (1979).The energy deposited in the fibroglandular voxels wasrecorded and then converted into dose according to theformula MGD = (cid:80) i E i M g (6)where E i is the energy deposited at the interaction event i ,and M g is the total fibroglandular breast mass.10 primary x-rays were emitted by an isotropic pointsource placed at 70 cm from the detector and collimatedto irradiate the entire detector. In order to replicate thetomosynthesis acquisition mode, a total of 25 projectionswere simulated from − ◦ to 24 ◦ , every 2 ◦ . The projec-tion at 0 ◦ replicates the mammographic acquisition. Thenumber of primary particles ensured a statistical uncer-tainty on the total dose of less than 0 .
22 23 27 28 29 x y z cm X ‐ ray source Breast support table
Detector
Breast
Model
Compression paddle S o u r c e ‐ D e t e c t o r d i s t a n c e c m Pa � ent Body ‐ S o u r c e ‐ c e n t e r o f r o t a t i o n d i s t a n c e c m Figure 5: Imaging geometry implemented in the Monte Carlo simula-tion: the x-ray source is placed at 70 cm from the detector, a 3 mmthick polyethylene terephthalate (PET) compression paddle was simu-lated and a large water cuboid was included to take into account thepatient-body backscatter. The x-ray field irradiated the breast model atdi ff erent angles (from -24 ◦ to 24 ◦ ). The center of rotation is placed at 65cm from the x-ray source. Drawing is not to scale and rotation of thedetector is not shown.Table 1: X-ray spectra used in the Monte Carlo simulation. HVL: 1sthalf value layer Breast Spectrum HVL / Rh - 27 kV 0 .
519 240-49 W / Rh - 28 kV 0 .
530 750-59 W / Rh - 29 kV 0 .
538 660-69 W / Rh - 30 kV 0 .
547 1870-79 W / Rh - 31 kV 0 .
557 12the method proposed by Sempau et al. (2001). Photoelec-tric interactions, and coherent and incoherent scatter wereincluded in the simulations without modifying the defaultcut range for photons (1 mm, corresponding to an energyof 2 .
45 keV and 2 .
88 keV for adipose and fibroglandulartissue, respectively).The x-ray spectra were modeled using the TASMICS8odel (Hernandez and Boone, 2014) by adjusting thethickness of the modeled rhodium filter to match thefirst half-value layer measured with a solid state detector(RaySafe X2-MAM sensor, Billdal, Sweden) in the mod-elled system as shown in Table 1. The Monte Carlo sim-ulations for estimating the MGD were performed twicefor each breast; once for each of the 45 BCT phantoms,and once each for the corresponding labelled DBToR re-constructions. In this way, the accuracy of the resultingpatient-specific dosimetry estimates could be assessed.
The results of the DBToR reconstruction, prior to thevoxel classification for estimation of breast density anddose, were compared to the baseline iterative MLTR re-construction algorithm and the LPD algorithm. The LPDalgorithm was trained using the noiseless version of thesimulated projections.
4. Results
All models were trained on a single NVidia Titan XGPU with 12 GB of memory and batch size 1, the numberof iterations N iter was set to 4 · for all models. Train-ing each model took approximately a week. At inferencestage, the reconstruction takes around 1 second for a sin-gle coronal slice. For DBToR trained on noise-free vir-tual phantom data we report the corresponding L loss,Structural Similarity Index (SSIM) (Wang et al., 2004)and Peak Signal-to-Noise Ratio (PSNR) on noise-free vir-tual phantom test data in Table 2. Lower L loss andhigher SSIM and PSNR values indicate better reconstruc-tion performance. For DBToR trained on noisy virtualphantom projections, where we used noise level N = N = , ,
12 in Table 3. In both tables, we report themean and standard deviation of the metrics obtained us-ing 3 cross-validation folds with 50% of the data used fortraining and 50% used for testing in each fold.The LPD algorithm is significantly outperformed in thenoise-free case. Visual inspection of the slices producedby LPD revealed that the LPD often reconstructs breast re-gions adjacent to the compression paddle very poorly forboth test and training data. In particular, it frequently failsto reproduce the flatness of the part of the skin surface which is in contact with the compression paddle. Sincethe same type of reconstruction artifacts appear for bothtraining and test data, we ruled out overfitting of LPDas the cause of the artifacts. While it is possible thata much larger version of LPD with more reconstructionblocks would learn to correct these artifacts, we will seethat DBToR resolves them without any further increase inthe number of parameters. Since LPD performed poorlyon noise-free projections (see Table 2), we excluded itfrom further training on noisy projections. The proposedDBToR algorithm outperforms the MLTR reconstructionalgorithm at all noise levels and for all metrics being con-sidered, while yielding visually more accurate reconstruc-tions as well. It is also interesting to note from Table 3 thatthe performance of DBToR at noise level N = N =
8, which corresponds to a 4 times higher photoncount. At noise level N =
12, the performance of MLTRis slightly below the performance of MLTR for the noise-free case. At the same time, DBToR at N =
12 reachescomparable level of performance to that of the DBToR onnoise-free data, which can be seen from slightly higherPSNR and slightly lower SSIM as compared to the noise-free version.For DBToR trained on noisy virtual phantom projec-tions and subsequently finetuned on deformed breast CTslices, where we used noise level N = Figure 6 shows a box-whisker plot of the absolute per-centage di ff erence in glandularity between the DBToR es-timate and the ground truth (GT). It can be seen that, onaverage, no bias is observed ( p -value 0 . able 2: Results on noise-free phantom projections, mean ± standard deviation (in bold best result) Model L -loss SSIM PSNRMLTR 0 . ± . · − . ± . · − . ± . · − LPD 0 . ± . · − . ± . · − . ± . . ± . · − . ± . · − . ± . · − Table 3: Results on noisy phantom projections for di ff erent noise levels N , mean ± standard deviation (in bold best result) Model L -loss SSIM PSNRMLTR ( N =
4) 0 . ± . · − . ± · − . ± . · − DBToR ( N = . ± . · − . ± . · − . ± . MLTR ( N =
8) 0 . ± . · − . ± . · − . ± . · − DBToR ( N = . ± . · − . ± . · − . ± . · − MLTR ( N =
12) 0 . ± . · − . ± . · − . ± . · − DBToR ( N = . ± . · − . ± . · − . ± . · − Table 4: Results on noisy BCT phantom projections (in bold best result)
Model L -loss SSIM PSNRMLTR 0 .
006 0 .
82 21 . . .
93 27 . The comparison between the MGD evaluated from theDBToR reconstructions and the GT is shown in the Bland-Altman plots of Figure 7, for both mammography andDBT geometries.As can be seen, no bias is observerd in the proposeddose estimation method ( p -value 0 . ff erences in the reconstructed fi-broglandular distribution obtained with the DBToR modelcompared to the GT, as shown in an example in Figure 8.For this case, the absolute di ff erence on the glandular-ity is 1 .
5% (namely 12 .
0% for the DBToR and 13 .
5% forGT). Thus, due to a higher amount of fibroglandular tis-sue closest to the x-ray source in the case of the GT, theradiation dose estimated by DBToR is lower by 13 .
5. Discussion and conclusion
We presented a deep learning-based method for the re-construction of DBT, which we call DBToR. The model isboth data driven and model-based, since the forward andbackprojection operators for a given DBT geometry area part of our neural network architecture and at the sametime the model is trained to reduce the tomographic arti-facts of the reconstruction. As training data we used twosources of data, one based on random samples with sta-tistical properties similar to real breast volumes, and onedataset of patient breast CT images that have been com-pressed with a finite element model to simulate the samebreast under compression in a DBT system.Compared to LPD, in DBToR we added the com-pressed breast thickness as prior information. Since thelimited angle causes a severely ill-posed problem, it wasexpected that this information would be definitely re-quired, and the experiments (Table 2) confirmed thatthe result dramatically degrades, compared to the ’full’DBToR, when this prior knowledge is not provided to thealgorithm. Requiring this additional information does notlimit the generalizability of the method, since it is readilyavailable in all DBT systems.The results indicate that the proposed algorithm outper-forms the MLTR iterative reconstruction in terms of re-10 ● − − − A b s o l u t e d i ff e r en c e i n g l andu l a r i t y pe r c en t age ( % ) GTDL th PercentileMedian25 th PercentileQ1 – (1.5 x IQR)Q3 + (1.5 x IQR) A b s o l u t e d i ff e r e n c e i n G l a ndu l a r i t y ( % ) -5-4-3-2-1012345 0 10 20 30 40 50 60 70 D i ff i n G l a ndu l a r i t y ( D B T o R ‐ G T ) ( % ) Glandularity (%)
Mean + 2SDMean (0.19)Mean ‐ Figure 6: (top) Box-whisker plot and (bottom) Bland-Altman plot of theabsolute di ff erence (DBToR - GT) of glandularity (in percentage). Thesymbols in the box-whisker plot denote outliers. construction quality for this application. Furthermore, thealgorithm generalizes well even when trained on a smalldataset and is robust to noise.The simulated acquisitions in this work used a mono-energetic beam and and did not include x-ray scatteredradiation, so it remains to be seen how our new recon-struction method will handle these factors. In practice thee ff ect of not modeling the spectrum will likely be minimalas regular filtered backprojection reconstructions also donot account for these physical e ff ects and apply a seriesof precorrection steps to the projection data instead, suchas the beam hardening correction described by Herman(1979). We foresee using the same approach to extendour method to work in clinical data.We have shown that the method achieves robust and ac-curate predictions of breast density, which is an importantmetric relating to masking and cancer risk. As opposed tocurrent density estimation methods based on mammogra-phy and DBT projections, which require assumptions andmodeling of the image acquisition process, the use of theimages produced by DBToR allows for a direct estima-tion of the amount of dense tissue present in the volume,resulting in estimates to within 3%. Accurate determina-tion of breast density opens up further opportunities forpersonalized risk-based screening. As is crucial for anaccurate dosimetric estimate, the location in the verticaldirection of the dense tissue is also estimated with ac-curacy, resulting in a state-of-the-art dosimetric estimate.It is known that current model dose estimates introducean average bias of 30%, and can misrepresent the actualpatient-specific dose by up to 120% (Dance et al., 2005;Sechopoulos et al., 2012; Hernandez et al., 2015). Incomparison, the results obtained here achieve errors be-low 20% with no systematic bias. True patient-specificdosimetry could be used, for the first time, to gather doseregistries, especially for screening, ensuring the optimaluse of this imaging technology, and allowing for continu-ous monitoring of dose trends and providing valuable datafor additional optimization and development of existingand new imaging technologies.The main limitation of the current work is that it workson a slice by slice basis rather than on a full 3D volume.With this simplification we were able to concentrate onthe network structure rather than on the logistics of han-dling the enormous datasets needed to train a 3D model.The logical next step for our algorithm is an extension11 % D i ff e r e n c e i n D o s e [ ( D L ‐ G T ) / G T ] Dose(nGy per million of photons)
Mammography
Mean + 2SDMeanMean ‐ (a) -20%-15%-10%-5%0%5%10%15%20%0,00 2,00 4,00 6,00 8,00 10,00 12,00 % D i ff e r e n c e i n D o s e [ ( D L ‐ G T ) / G T ] Dose(nGy per million of photons)
Tomosynthesis
Mean + 2SDMeanMean ‐ (b)Figure 7: Bland-Altman plot of the di ff erence in dose estimates resulting from the DBToR reconstruction and the ground truth, in percentage, forboth mammography (a) and DBT (b). The red line represents the mean, while the two blue-dashed lines represent the 95% limits of agreement. GT DL
Coronal View Side View
Figure 8: Example of fibroglandular tissue distribution for an outliercase: the ground truth (in white) depicts a higher amount of fibroglan-dular tissue in the top breast layer (i.e., facing the x-ray tube), while theDBToR model (in green) predicts a fibroglandular distribution spreadtowards the anterior, and more inferior, part of the breast. to fully 3D data instead of 2D slices. From there on, wecould extend the model-based parts of the deep learningnetwork instead of starting from precorrected projectiondata, as in current filtered backprojection methods, byincluding a polychromatic x-ray spectrum, x-ray scatter,and other relevant factors in the forward model. It couldalso be valuable to optimize the network output specif-ically for artificial intelligence reading by training thecomposite network in an end-to-end fashion. Finally, thereconstruction of a diagnostic-quality volume, for inter-pretation by radiologists, with the potential for the highervertical resolution obtained here, would be a valuable im-provement for clinical performance.To conclude, we created a deep learning based recon- (a) (b)(c) (d)Figure 9: (a) Coronal slice of a breast CT phantom, (b) MLTR recon-struction, (c) DBToR reconstruction, and (d) classification of DBToRreconstruction. struction for DBT that was able to achieve accurate pre-dictions of breast density and from there an accurate cal-culation of patient specific MGD.
Acknowledgments
The research was partially supported by NationalCancer Institute, National Institutes of Health:12 a) (b)(c) (d)Figure 10: (a) Axial slice of a breast CT phantom, (b) MLTR recon-struction, (c) DBToR reconstruction, and (d) classification of DBToRreconstruction.
References
Adams, R., Bischof, L., 1994. Seeded Region Growing.IEEE Transactions on Pattern Analysis and MachineIntelligence 16.Adler, J., ¨Oktem, O., 2018. Learned Primal-Dual Recon-struction. IEEE Transactions on Medical Imaging 37,1322–1332. doi: .Arana Pe˜na, L.M., Fedon, C., Garcia, E., Diaz, O., Longo,R., Dance, D.R., Sechopoulos, I., 2020. Monte Carlodose evaluation of di ff erent fibroglandular tissue distri-bution in breast imaging, in: Van Ongeval, C., Mar-shall, N., Bosmans, H. (Eds.), 15th International Work-shop on Breast Imaging (IWBI2020), SPIE. p. 76.doi: .Arridge, S., Maass, P., ¨Oktem, O., Sch¨onlieb, C.B.,2019. Solving inverse problems using data-drivenmodels. Acta Numerica 28, 1174. doi: .Bakic, P.R., Zhang, C., Maidment, A.D.A., 2011. De-velopment and characterization of an anthropomor-phic breast software phantom based upon region-growing algorithm. Medical Physics 38, 3165– 3176. URL: http://doi.wiley.com/10.1118/1.3590357 , doi: .Bezdek, J.C., 1981. Pattern recognition with fuzzy objec-tive function algorithms. Plenum, New York, NY.Boone, J.M., Chavez, A.E., 1996. Comparison ofx-ray cross sections for diagnostic and therapeu-tic medical physics. Medical Physics 23, 1997–2005. URL: http://doi.wiley.com/10.1118/1.597899 , doi: .Caballo, M., Boone, J., Mann, R., Sechopoulos, I., 2018.An unsupervised automatic segmentation algorithm forbreast tissue classification of dedicated breast com-puted tomography images. Medical Physics 45, 2542–2559.Dance, D.R., Sechopoulos, I., 2016. Dosimetry in x-ray-based breast imaging. Physics in Medicine and Biol-ogy 61, R271–R304. doi: .Dance, D.R., et al., 2005. Breast dosimetry usinghigh-resolution voxel phantoms. Radiation ProtectionDosimetry 114, 359–363. doi: .Fedon, C., Caballo, M., Longo, R., Trianni, A., Se-chopoulos, I., 2018a. Internal breast dosimetry in mam-mography: Experimental methods and Monte Carlovalidation with a monoenergetic x-ray beam. MedicalPhysics 45, 1724–1737. doi: .Fedon, C., Caballo, M., Sechopoulos, I., 2018b. Internalbreast dosimetry in mammography:Monte Carlo val-idation in homogeneous and anthropomorphic breastphantom with a clinical mammography system. Medi-cal Physics 45, 3950–3961. doi: .Fedon, C., et al., 2019. Monte Carlo study on optimalbreast voxel resolution for dosimetry estimates in digi-tal breast tomosynthesis. Physics in Medicine & Biol-ogy 64, 015003. doi: .F¨ornvik, D., Timberg, P., Zackrisson, S., Tingberg, A.,Petersson, H., Dustler, M., 2018. Towards determina-tion of individual glandular dose, in: Chen, G.H., Lo,J.Y., Gilat Schmidt, T. (Eds.), Medical Imaging 2018:Physics of Medical Imaging, SPIE. p. 3. doi: .13 a) (b) (c) (d)Figure 11: (a) Sagittal slice of a breast CT phantom, (b) MLTR reconstruction, (c) DBToR reconstruction, and (d) classification of DBToRreconstruction. Garc´ıa, E., Fedon, C., Caballo, M., Mart´ı, R., Sechopou-los, I., Diaz, O., 2020. Realistic compressed breastphantoms for medical physics applications, in: VanOngeval, C., Marshall, N., Bosmans, H. (Eds.), 15th In-ternational Workshop on Breast Imaging (IWBI2020),SPIE. p. 73. doi: .Geeraert, N., 2014. Quantitative evaluation of fibroglan-dular tissue for estimation of tissue-di ff erentiated ab-sorbed energy in breast tomosynthesis. Ph.D. thesis.KU Leuven. Faculty of medicine.Hammerstein, R., Miller, D.W., White, D.R., Masterson,E., Woodard, H.Q., Laughlin, J.S., 1979. Absorbed Ra-diation Dose in Mammography. Radiology 130, 485–491. URL: http://pubs.rsna.org/doi/10.1148/130.2.485 , doi: .Herman, G.T., 1979. Correction for beam hardening incomputed tomography. Physics in Medicine and Biol-ogy 24, 008. URL: https://iopscience.iop.org/article/10.1088/0031-9155/24/1/008 , doi: .Hernandez, A., Boone, J., 2014. Tungsten anode spectralmodel using interpolating cubic splines: unfiltered x-ray spectra from 20 kV to 640 kV. Medical Physics 41,042101.Hernandez, A.M., Seibert, J.A., Boone, J.M., 2015.Breast dose in mammography is about 30% lower whenrealistic heterogeneous glandular distributions are con-sidered. Medical Physics 42, 6337–6348. doi: . Johnsen, S.F., et al., 2015. NiftySim: A GPU-based non-linear finite element package for simulation of soft tis-sue biomechanics. Int J Comput Assist Radiol Surg 10,1077–1095. doi: .Lau, B.A., Reiser, I., Nishikawa, R.M., Bakic, P.R.,2012. A statistically defined anthropomorphic soft-ware breast phantom. Medical Physics 39, 3375–3385. URL: http://doi.wiley.com/10.1118/1.4718576 , doi: .Loshchilov, I., Hutter, F., 2019. SGDR: Stochastic gra-dient descent with warm restarts, in: 5th InternationalConference on Learning Representations, ICLR 2017 -Conference Track Proceedings, pp. 1–16.McCormack, V.A., 2006. Breast Density andParenchymal Patterns as Markers of Breast Can-cer Risk: A Meta-analysis. Cancer Epidemi-ology Biomarkers & Prevention 15, 1159–1169.URL: http://cebp.aacrjournals.org/cgi/doi/10.1158/1055-9965.EPI-06-0034 , doi: .Michielsen, K., 2016. Maximum a posteriori reconstruc-tion for digital breast tomosynthesis. Ph.D. thesis. KULeuven. Faculty of medicine. Leuven.Michielsen, K., Nuyts, J., Cockmartin, L., Marshall, N.,Bosmans, H., 2016. Design of a model observer to eval-uate calcification detectability in breast tomosynthesisand application to smoothing prior optimization. Medi-cal Physics 43, 6577–6587. doi: .Miller, K., Joldes, G., Lance, D., Wittek, A., 2007. TotalLagrangian explicit dynamics finite element algorithm14or computing soft tissue deformation. Communica-tions in Numerical Methods in Engineering 23, 121–134. doi: .Moriakov, N., Michielsen, K., Adler, J., Mann, R., Se-chopoulos, I., Teuwen, J., 2019. Deep learning frame-work for digital breast tomosynthesis reconstruction,in: Bosmans, H., Chen, G.H., Gilat Schmidt, T. (Eds.),Medical Imaging 2019: Physics of Medical Imaging,SPIE. p. 220. doi: .Niklason, L.T., et al., 1997. Digital tomosynthesis inbreast imaging. Radiology 205, 399–406. doi: .Nuyts, J., Man, B.D., Dupont, P., Defrise, M., Suetens,P., Mortelmans, L., 1998. Iterative reconstruction forhelical CT: a simulation study. Physics in Medicineand Biology 43, 729–737. doi: .del Palomar, A.P., Calvo, B., Herrero, J., Lopez, J.,Doblare, M., 2008. A finite element model to accu-rately predict real deformations of the breast. Med EngPhys 30, 1089–1097. doi: .Plevritis, S.K., et al., 2018. Association of screeningand treatment with breast cancer mortality by molecu-lar subtype in US women, 2000-2012. JAMA - Journalof the American Medical Association 319, 154–164.doi: .Reiser, I., Nishikawa, R.M., 2010. Task-based assessmentof breast tomosynthesis: E ff ect of acquisition parame-ters and quantum noisea). Medical Physics 37, 1591–1600. URL: http://doi.wiley.com/10.1118/1.3357288 , doi: .Rodriguez-Ruiz, A., et al., 2017. New reconstruction al-gorithm for digital breast tomosynthesis: better imagequality for humans and computers. Acta Radiologica .Sechopoulos, I., Bliznakova, K., Qin, X., Fei, B., Feng,S.S.J., 2012. Characterization of the homogeneous tis-sue mixture approximation in breast imaging dosime-try. Medical Physics 39, 5050–5059. doi: . Sempau, J., Sanchez-Reyes, A., Salvat, F., Tahar, H.,Jiang, S., Fernandez-Varea, J., 2001. Monte Carlo sim-ulation of electron beams from an accelerator head us-ing PENELOPE. Physics in Medicine and Biology 46,1163–1186.Wang, Z., Bovik, A., Sheikh, H., Simoncelli, E.,2004. Image Quality Assessment: From Er-ror Visibility to Structural Similarity. IEEETransactions on Image Processing 13, 600–612.URL: http://ieeexplore.ieee.org/document/1284395/ , doi: .Wellman, P.S., 1999. Tactile imaging. Ph.D. thesis. Har-vard University.Zackrisson, S., et al., 2018. One-view breast tomosyn-thesis versus two-view mammography in the malmbreast tomosynthesis screening trial (mbtst): a prospec-tive, population-based, diagnostic accuracy study. TheLancet Oncology 19, 1493–1503.Zhang, C., Bakic, P.R., Maidment, A.D.A., 2008. Devel-opment of an anthropomorphic breast software phan-tom based on region growing algorithm, in: Miga, M.I.,Cleary, K.R. (Eds.), Medical Imaging 2008: Visual-ization, Image-Guided Procedures, and Modeling, p.69180V. doi: .Zuley, M.L., et al., 2013. Digital breast tomosynthesisversus supplemental diagnostic mammographic viewsfor evaluation of noncalcified breast lesions. Radiology266, 89–95. doi:10.1148/radiol.12120552