Three-compartment-breast (3CB) prior-guided diffuse optical tomography based on dual-energy digital breast tomosynthesis (DBT).

Diffuse optical tomography (DOT) is a non-invasive functional imaging modality that uses near-infrared (NIR) light to measure the oxygenation state and the concentration of hemoglobin. By complementarily using DOT with other anatomical imaging modalities, physicians can diagnose more accurately through additional functional image information. In breast imaging, diagnosis of dense breasts is often challenging because the bulky fibrous tissues may hinder the correct tumor characterization. In this work, we proposed a three-compartment-breast (3CB) decomposition-based prior-guided optical tomography for enhancing DOT image quality. We conjectured that the 3CB prior would lead to improvement of the spatial resolution and also of the contrast of the reconstructed tumor image, particularly for the dense breasts. We conducted a Monte-Carlo simulation to acquire dual-energy X-ray projections of a realistic 3D numerical breast phantom and performed digital breast tomosynthesis (DBT) for setting up a 3CB model. The 3CB prior was then used as a structural guide in DOT image reconstruction. The proposed method resulted in the higher spatial resolution of the recovered tumor even when the tumor is surrounded by the fibroglandular tissues compared with the typical two-composition-prior method or the standard Tikhonov regularization method.


Introduction
Diffuse optical tomography (DOT) is a non-invasive functional imaging modality that uses near-infrared light in the spectral range of 600-1200 nm and reconstructs the oxygenation state and the concentration of hemoglobin in tissues for functional imaging [1,2]. Certain chromophores (oxy,deoxy-hemoglobin) determine the optical window of the tissue by dominantly absorbing the near-infrared at specific wavelengths [3]. DOT is useful for distinguishing malignant tumors from benign ones since the malignant tumor typically contains relatively high hemoglobin concentration due to the development of the blood vessels [4,5]. A breadth of DOT studies have been performed particularly in the brain and breast imaging for diagnosis and also for monitoring of chemotherapy [6][7][8][9]. There are clinical evidence that the functional image information acquired by DOT significantly enhances diagnostic performance [10][11][12][13][14]. Since an iterative scheme was first introduced in DOT image reconstruction using diffusion based forward model [15], a number of advanced reconstruction approaches have been developed such as non-iterative exact reconstruction method based on joint-sparsity [17,18] and time-resolved diffuse optical tomography [19] to name a few. However, the ill-posedness of DOT image reconstruction originated mostly from the intrinsic scattering nature causes rather poor spatial resolution and low contrast in the reconstructed images, and developing advanced image reconstruction methods is still an active area of research.
To ameliorate the ill-posedness, the structural prior guided methods using other imaging modalities have been investigated. The anatomical structure information can help the optical tomographic image reconstruction to enhance the spatial resolution of the region-of-interest (ROI) [20][21][22]. There have been developed a host of guiding approaches for DOT image reconstruction such as MRI-guided, CT-guided, DBT-guided, mammography-guided, and ultrasound-guided methods [20][21][22][23][24][25]. However, these co-registered imaging approaches have the common issue that the pixel value of the prior images is not directly related to the optical properties of the tissue. Therefore, using the structure prior derived from the co-registered image as a "hard" prior can introduce unwanted bias although the total number of unknowns in the inverse problem may decrease [20]. As an alternative, the "soft" prior method can be used. It uses a regularization scheme for incorporating the structural prior thus reducing the likelihood of spatial bias [21,22].
In breast imaging, there have been attempts using DBT images or mammographic images to guide the DOT reconstruction using the tissue composition [24,25]. In such x-ray imaging modalities, one can derive compositional tissue density maps (e.g., adipose and fibroglandular density fraction) using dual-energy scans [26]. Subsequently, it is possible to convert the obtained tissue density maps into the tissue optical attenuation maps through the calibrated or predefined conversion relationship [5]. Nevertheless, the above tissue decomposition has limitations particularly for the patients who have bulky fibrous breast tissue typically common in dense breasts. Even for a relatively large size of tumor, it is still hard to characterize the suspicious lesion after the tissue decomposition since the tumor has a similar x-ray attenuation property with the fibrous tissue [27][28][29].
Laidevant et al. proposed a three-compartment breast (3CB) decomposition technique using the dual-energy mammogram with the breast thickness information to obtain three selective material images (i.e., water, lipid, and protein) [30]. Based on the material specific images, biologically relevant image analysis in the context of tumor malignancy can be conducted. As mentioned earlier, malignant tumors usually contain relatively higher water content than benign ones due to the vascularization [4,5]. From an optical imaging point of view, water image is likely to reflect the concentration of chromophores (oxy,deoxy-hemoglobin) thereby providing useful information on the tumor malignancy. However, due to the 2D projection nature of mammography, this method may not be directly applicable for DOT image reconstruction in the form of an image prior. To fully exploit 3D optical property information of the breast, X-ray tomographic imaging such as DBT that can offer volumetric structural information as a prior would be a better alternative.
In the previous work, many researchers including us have successfully developed a DBT/DOT multi-imaging modality system for clinical implementation [12][13][14]25]. However, the acquired DBT image was either used as a selection guide to incorporate effective source-detector channels for DOT image reconstruction or used as a structural prior that does not necessarily correlate with the tissue optical properties. Thus, there were difficulties in diagnosis particularly for the case of dense breasts. To overcome this limitation, we propose a dual-energy DBT based 3CB decomposition technique in this work. We will not only use the structural information of DBT as a prior but also use the material information of the breast for DOT image reconstruction.
To show the feasibility of the proposed method, we have conducted a simulation study in this work since there is no clinically available dual-energy DBT system yet. We used a simple material box phantom and also realistic 3D numerical breast phantom, VICTRE. The dual-energy x-ray projections for the DBT scan system were acquired through the GATE Monte-Carlo simulation. By utilizing the obtained material selective volumes, the optical attenuation coefficient of the breast was derived from the linear sum of the contributions of the chromophores (oxy,deoxyhemoglobin, water, and lipid) [5]. In the "soft" prior framework, a Laplacian-type matrix L was constructed for guiding the suspicious lesion [21,22]. The NIRFAST Matlab toolkit [31] was used to generate numerically simulated optical measurement data in the specific source-detector geometry of the prototype system which has been used in our earlier work [13,14]. After the reconstruction, the results were compared with the two types of tissue composition prior methods [24], and the standard Tikhonov regularization. The two-composition-prior method decomposes the breast image into fibrous and adipose tissue fraction images using dual-energy measurements and the acquired fibrous tissue fraction image is utilized as a prior image. The standard Tikhonov regularization method uses L1-norm of the image for regularization to achieve noise suppression during iterations. Furthermore, the false-positive prior guided reconstruction was also conducted to check the robustness of the method.
In the Method section, the concepts of the 3CB decomposition technique and its associated calibration procedures are briefly explained. The numerical simulation steps of optical data acquisition and the proposed DOT image reconstruction algorithm are then described. Results of the DBT simulations are all included in the Method section. Since the focus of the paper is on DOT image reconstruction using the DBT image as a prior, we feel that it is more straightforward to explain DBT image part in the Method section. In the Results section, the maximum intensity projections of the reconstructed DOT images are compared. A quantitative analysis is also performed to validate the performance of the proposed method in reconstructing the tumor in the bulky fibrous region. Finally, the key findings related to the proposed method are summarized in the Discussion section including the interpretation of the results and the optimization parameter analysis.

Decomposition technique using dual-energy digital breast tomosynthesis
Breast tissue can be decomposed into three basis materials (water, lipid, and protein) in terms of X-ray imaging physics [32]. In polychromatic dual-energy x-ray imaging, the logarithm of the measurement at each energy is modeled by a nonlinear function [33,34]. It is often expanded by a power series of three variables: the high-energy measurements A HE , the ratio of low-and high-energy measurements R, and the total thickness of the object T [30]. Thus, the thickness of the basis materials can be presented in the calibration procedure as following.
, where a αβγ is the calibration coefficients to be determined during the calibration and t j is the thickness of each basis material. Using this technique, we obtained three compositional material volumes from the dual-energy digital breast tomosynthesis. For the calibration procedure, a total of 40mm-thick and also of 20mm-thick stack of digital calibration step phantoms (as shown in Fig. 1) were prepared. Each square unit has a thickness of 10 mm and a side width of 13 mm, and they are stacked in various combinations while maintaining the total thickness of the calibration phantoms. Each unit is composed of one of lipid (C 55 H 104 O 6 ), protein (C 100 H 159 N 26 O 32 S 0.7 ), or water (H 2 O).
To simulate a realistic DBT scanning system, the geometry of a prototype DBT system developed by Korea Electrotechnology Research Institute (KERI) was used and the scanning parameters are summarized in Table 1. The polychromatic projection data were acquired by use of the GATE Monte-Carlo simulation [35]. The X-ray spectrum was generated by the Spektr 2.2 software for high-energy(43 kV) and low-energy(28 kV) on the tungsten target with a 0.05mm-thick Rh filter. To increase the accuracy of material decomposition, an anti-scatter grid having the grid ratio of 5:1 was employed for reducing the scatter signals. For the detector materials, 0.2mm-thick selenium was used. Figure 1 shows the example projection images of the calibration phantoms. To validate the proposed material decomposition method, a 40mm-thick material box phantom was designed as shown in Fig. 2(a). This phantom is composed of a lipid rectangular box, protein ball, water ball, and blood balls having various radii of 10.0, 7.0, and 5.6 mm. The composition of the blood was determined from the reported value in [36]: 79.0 w/w of water, 0.6 w/w of lipid and 19.6 w/w of protein (w/w indicate mass fraction percent). Figure 2 , and (f) show the water, lipid, and protein-selective slices at the in-focus depth of each structure in the reconstructed images, respectively. The images were successfully decomposed into three basis material volumes including the blood balls as classified into the water volume map.  The number of projection 11

Generation of the prior from the VICTRE digital breast phantom
A realistic 3D compressed digital breast phantom with its compressed thickness of 40 mm containing a 10mm-thick tumor inside was generated using the VICTRE software toolkit [37,38].
We simulated a dense breast where the tumor is embedded in the fibrous tissue and also a fatty breast with tumor for comparison. The composition of the breast tissue was determined according to the tissue information reported in [36] and the used compositional fractions are summarized in Table 2. Projections were then acquired under the same scanning conditions with the ones used in the calibration procedure. The reconstructed DE-DBT images and the resulting material decomposition maps from the reconstructed DBT images are shown in Fig. 3. As shown in Fig. 3(f), it appears that the protein fractional map is overestimated near the breast boundary, whereas other material maps are reasonably estimated. This is thought to be due to the relatively small thickness in the breast boundary part compared to the thicknesses of the calibration points. Those material volumes were then used to estimate the optical attenuation of the breast, based on the known optical properties of major chromophores at a specific wavelength (830 nm) using the generic Eq. (2) as following [5].
, where µ a is the tissue optical absorption coefficient that we want to estimate. The optical absorption coefficients of chromophores (µ a,oxy , µ a,deoxy , µ a,water , µ a,fat ) and the HGb oxygen saturation S were determined from the reported values in [5,[39][40][41]. W and F indicate the water and lipid volume fractions respectively. Meanwhile, the miner chromophores of the breast such as melanin monomers or bilirubin were neglected for their contributions to the optical imaging are known to be substantially lower than the major ones [5]. Additionally, the protein chromophores including amino acids, tryptophan, tyrosine and cysteine were excluded from the consideration due to the negligible optical absorption contribution in the range of wavelength of 830 nm. Rather, they have significant absorption properties near the wavelength of 280 nm, or ultraviolet region [42]. To determine the blood volume fraction k, the reported clinical substances fraction data of breast tissues [43][44][45] were used with an assumed proportionality with the water volume. Due to the characteristic water contents between tissues, the suspicious lesion can be segmented using a threshold from the recovered breast optical attenuation map as shown in Fig. 4(b). The segmented image was then converted to mesh domain to use as prior information for further "soft-prior" scheme of regularization.

Numerical simulation for optical tomography using NIRFAST software
The meshes were constructed with tetrahedral elements for the finite element method (FEM) using Near-Infrared Fluorescence and Spectral Tomography (NIRFAST), an open-source MATLABbased software package [31]. The forward mesh was constructed with 13,512 nodes and 2.5 mm of node distance. For the reconstruction mesh, 8,029 nodes and 3 mm of node distance were      used. As illustrated in Fig. 4(d), the 55 pairs of source-and-detector were used in our numerical simulation and its spatial distributions followed the top-bottom projection types of a prototype DOT system [13,14]. The optical properties of breast tissues were predefined using clinical breast measurements in [12] and reported in Table 3. A diffusion-based forward model [15] was adopted to generate the optical measurements by numerically solving it through the FEM on the forward mesh with additional random noise at its level of 1%.

3CB-prior guided reconstruction algorithm
According to the previous procedures, we derived a prior mesh that is relevant to the tumor malignancy based on the estimated optical attenuation map of the breast. To suitably exploit the prior mesh for structural prior guided DOT reconstruction, a Laplacian-type matrix L was constructed for regularization by the following equation.
, where N is the total number of the nodes, i and j represents the node indices. The matrix L effectively smooths the region-of-interest while allowing the discontinuity of the boundaries of different regions [21,22]. The objective function of the guided DOT reconstruction algorithm is given by Eq. (4).
, where Ω(µ) is an objective function with optical properties µ; Ψ d is a measured data; f is a forward operator;λ is regularization parameter;L is the matrix for the structural prior. The first term represents the model error or a data fidelity, and the second term represents the deviation of the solution from the prior knowledge working as a regularizer. For the optimization, the iterative Gauss-Newton approach was used. In this under-determined inverse problem, the general update equation for the k-th iteration can be expressed as following.
, where ∆µ k is the update for the optical properties at the k-th iteration, and J k is the Jacobian matrix. The absorption coefficients µ a and scattering coefficients µ s were recovered through the reconstruction, and µ a was used as a target unknown for evaluating the reconstruction performance quantitatively. Although there exist noniterative image reconstruction approaches such as truncated singular value decomposition method [16] or sparsity-exploiting algorithms [17,18], we adopted the popularly used iterative approach in this work and focused on highlighting the use of image prior acquired by the proposed dual-energy DBT. We do not, however, exclude chances of such noniterative image reconstruction algorithms utilizing the proposed image prior. In this work, the iteration stops when change of objective function between successive iterations become smaller than 2% of the previous change. In addition, the slope of each log intensity times distance and the slope of phase with respect to distance were calculated, in Fig. 5, to estimate initial bulk optical properties using the homogeneous fitting algorithm in [46].

Regularization parameter
For testing the effects of the regularization strength, the value of parameter λ was chosen to be one of 0.1, 1, and 10. We then compared the proposed method with the standard Tikhonov regularization method and the two-composition-prior method i.e., fibroglandular-adipose tissue priors. The ground-truth tissue composition was used to generate the two-composition prior mesh. To efficiently analyze the accuracy of the recovered absorption coefficient, the maximum intensity projection (MIP) of the reconstructed volume images were compared. For a quantitative comparison, the RMSE of the recovered absorption coefficients of the fibrous region were compared between the ground-truth. The line-profiles of the recovered tumor were then displayed for each λ.
where µ i is the recovered optical coefficient at pixel i, µ o i is the ground-truth optical coefficient, and n is the number of pixels.

False-positive analysis
To address the robustness of the proposed algorithm, the inaccurate structural prior guided reconstruction was intentionally conducted. The new prior mesh was created by shifting the original suspicious lesion to the bottom-right of the breast. This inaccurate prior was used to reconstruct the optical measurement obtained under the same conditions as above. The influence of the false-positive prior was analyzed through the comparison with the ground-truth for each λ.

Comparisons between reconstructed results
The pseudo-colored MIP images of the reconstructed absorption coefficient of fatty and dense breasts are displayed in Fig. 6 and 7. The first row shows the results of the proposed-prior method and the second row shows the results of the two-composition prior method. The last row displays the results of the Tikhonov regularization. For both cases, the Tikhonov regularization method turns out to fail in characterizing the tumor for all λ used in this work. Other two prior-guided methods perform better in tumor characterization. Yet, the two-composition prior method shows similar tumor contrast and resolution in Fig. 6 while it shows weaker tumor contrast and poorer resolution than the proposed method at λ = 0.1 and 10 in Fig. 7. The proposed method shows the highest tumor contrast and resolution for both cases. But, in Fig. 7, it visually slightly better than other methods at λ = 1. The data fidelity term and the regularization term in the objective function compete in terms of noise characteristics. At λ = 0.1, it seems that the regularization term rarely contributes. By increasing the value of λ, image noise can be suppressed with the target contrast being possibly degraded particularly at a suboptimal point. Thus, the result of λ=1 in Fig. 7, the proposed method may be regarded as a suboptimal point in our study. At λ = 10, the noise is substantially reduced while respecting the structural prior. The later quantitative studies were done only for dense breast results what we want to show better performance than others as a research goal. To competitively quantify the contrast and resolution of tumor between compared methods, the structural similarity index of normalized MIP image was also calculated with a ground-truth image along the λ in Table 4.  In Fig. 8, the magnitudes of the RMSE are plotted in bar graphs along the values of λ. The bar graphs were compared in the order of Tikhonov regularization, the two-composition prior, and the proposed methods in orange, green and purple colors with increasing λ. In Fig. 9 the pixel absorption coefficient values (mm-1) across the tumor were plotted versus the position x (mm) as a line-profile. The line-profiles were overlayed for comparison with the ground-truth. The red, blue, green profiles represent the results of the proposed, the two-composition prior and Tikhonov regularization methods, respectively.

Influence of the false-positive prior
We also conducted a case study where a false-positive structural prior exists to test whether the regularization term overwhelms or not. In Fig. 9, the results are shown along increasing λ. The falsely guided region is highlighted by a red dashed circle, whereas the true region is marked by a     black dashed circle. There was no significant increase in the absorption coefficient in the falsely guided region even though the true target signal was not enhanced as much as in the correctly guided case as shown in Fig. 6.

Discussion and conclusions
In Fig. 2, the material box phantom images were decomposed into the water, lipid, and protein maps using the proposed dual-energy DBT based 3CB decomposition technique. It is noted that the blood balls were mostly decomposed into the water volume fraction, and indicates the blood in the breast would be represented in the water image in practical cases. Although the fibrous tissue also contains substantial fraction of water, malignant tumors usually have vascularization that would lead to concentrated water fractions. Therefore, the contrast of the water map is considered useful information in the context of tumor malignancy. In Fig. 3, different types of breast tissues were decomposed into the water, lipid, and protein maps. In dual-energy DBT reconstructed images of the breast, the tumor is hard to locate due to its low contrast relative to the fibrous background tissue. However, the tumor contrast is much enhanced in the water fraction maps thereby increasing the detectability of the tumor. In the meantime, the material decomposition was rather poorly estimated in the breast boundary regions due to the relatively small thickness of the breast boundary. In the future work, we are going to investigate more relevant calibration method in conjunction with the breast thickness estimation particularly near the breast boundary. A deep-neural network-based breast thickness estimation algorithm in DBT projections is also in the progress [47].
In Fig. 6 and 7, the Tikhonov regularization method failed to appropriately characterize the tumor. This is because the Tikhonov regularization does not incorporate the structural prior but simply denoises the images acquired by the data-fidelity-driven iterative reconstruction. Structural prior is believed to guide the solution to physically nontrivial point in the solution space. The DOT imaging problem in this work employs a top-to-bottom projection geometry resulting in a limited-angle inverse problem which is typically considered a highly ill-posed problem. An appropriate structural prior information can help the image reconstruction converge into a desirable solution. Although the Tikhonov regularization is generally known to provide a smoothed solution with its structure edges conserved in angularly fully-sampled tomographic image reconstruction, it has a limitation in handling the ill-posedness of the DOT imaging problem under our interest.
The two-composition prior method is also limited in that the structural prior contains not only the tumor but also the fibrous background tissue in the dense breast. In the fatty breast, however, tumor may have a substantially high contrast in the fibrous background and the two-composition prior method can sufficiently guide the tumor accordingly. Figure 6 indeed supports that the two-composition prior method also works well in the fatty breast case. On the other hand, the proposed method provides more precise structural prior thereby increasing the tumor contrast and resolution in both fatty and dense breast cases. It is the key advantage of the proposed method over the existing ones that the tumor can better be localized in the water map through the 3CB decomposition method even in the dense breast case. Underlying physics is that the tumor and the fibrous background usually have different relative water compositions that may not be readily discernable in the two-composition maps. Quantitative analysis consistently supports the performance of the proposed method over the other two algorithms as shown in Figs. 8, 9 and Table 4.
It was empirically found that the value of λ between 1 to 10 would be working reasonably well although we only show the results corresponding to 0.1, 1, and 10. Our aim was on illustrating the trend in the image quality change rather than tightly optimizing the value in this work. A judicious way of choosing the value of λ in the real experimental environment would be left as our future work. As shown in Fig. 10, there was no significant increase in the absorption coefficients in the red dashed circle region defined by a false-positive prior. This agrees well with the existing false-positive analysis of the soft prior methods [21,22]. Therefore, we believe the proposed method has a practical utility for enhancing image quality of the breast DOT imaging applications.   In this work, the scattering coefficient is not regarded as an active biological target as the absorption coefficient in the context of tumor malignancy since the scattering coefficients are not much different among the tissue types. On the other hand, the absorption coefficient difference between the tumor and the fibrous tissue is primarily due to the distinct hemoglobin concentration along with the tumor malignancy. Figure 11 shows the reconstructed images of scattering coefficients of the dense breast phantom in our study. The image contrast is much weaker than the absorption contrast, and the difference of the reconstructed images between the proposed method and the two-composition method is rather negligible.  Fig.11. The maximum intensity projections of the reconstructed scattering coefficients of the dense breast using different methods with different regularization parameters: the first row represents the results of the proposed prior method, the second row the results of the two-composition-prior method, and the third row represents the standard Tikhonov reconstruction results. The window level is [0.68 0.90]. In summary, we proposed an effective approach to constructing more accurate prior information for the dense breast patients using dual-energy DBT. It was successfully demonstrated that the dual-energy DBT based 3CB technique can provide material-selective volume images which can be utilized for optical property conversion. The utility of the proposed method was shown in-silico in the structural prior-guided DOT image reconstruction framework. Particularly, it is encouraging that the contrast and spatial resolution of the tumor can be increased in a dense breast, which is typically considered a diagnostically challenging case.