Target depth-regularized reconstruction in diffuse optical tomography using ultrasound segmentation as prior information

: Ultrasound (US)-guided diﬀuse optical tomography (DOT) is a promising non-invasive functional imaging technique for diagnosing breast cancer and monitoring breast cancer treatment response. However, because larger lesions are highly absorbing, reconstructions of these lesions using reﬂection geometry may exhibit light shadowing, which leads to inaccurate quantiﬁcation of their deeper portions. Here we propose a depth-regularized reconstruction algorithm combined with a semi-automated interactive neural network (CNN) for depth-dependent reconstruction of absorption distribution. CNN segments co-registered US to extract both spatial and depth priors, and the depth-regularized algorithm incorporates these parameters into the reconstruction. Through simulation and phantom data, the proposed algorithm is shown to signiﬁcantly improve the depth distribution of reconstructed absorption maps of large targets. Evaluated with 26 patients with larger breast lesions, the algorithm shows 2.4 to 3 times improvement in the top-to-bottom reconstructed homogeneity of the absorption maps for these lesions.


Introduction
Diffuse optical tomography (DOT) is a non-invasive functional imaging technique that reconstructs the optical properties of tissue [1]. Utilizing multiple wavelengths in the near infrared (NIR) spectrum enables quantitative estimations of oxygenated, deoxygenated, and total hemoglobin concentrations, which are directly related to tumor angiogenesis [2][3]. However, due to intense light scattering inside biological tissue, localizing a lesion area becomes challenging, and the DOT reconstruction problem becomes ill-posed and ill-conditioned. Other imaging modalities, such as ultrasound (US) [4][5][6], x-ray CT [7][8], and magnetic resonance imaging (MRI) [9][10], have been introduced to provide prior information to guide DOT reconstruction. As a result, the reconstruction is less ill-posed due to the significantly reduced number of imaging voxels. Several clinical studies have been reported for breast cancer diagnosis [4][5][6][7][9][10].
However, to reliably reconstruct the optical properties of MRI, CT, or US-identified breast lesions, a regularization scheme should be incorporated to constrain the solution space and obtain a quality image. The general approach is based on prior anatomical information from high resolution imaging modalities. Implementation of spatially encoded regularization using images obtained from MRI has been investigated, where a two-step image reconstruction procedure was required [11][12]. In the first step, a suspicious lesion was segmented from the rest of the tissue. Then this prior information was imposed to regularize the matrix inversion. Newer approaches include direct regularization imaging (DRI) in DOT [13] and depth compensation in fluorescence molecular tomography [14]. In DRI, the gray-scale values of MRI images are directly used to regularize the reconstruction, without the need for image segmentation [13]. For transmission-mode fluorescence molecular tomography, a depth compensation algorithm is introduced to reconstruct inclusions in deep tissues [14].
Our group has developed a US-guided DOT system and imaging algorithms that use coregistered US images to locate a breast lesion in reflection geometry and to segment the lesion into a fine-mesh lesion region of interest (ROI) and a coarse-mesh background tissue region [15]. This dual-mesh scheme improves the ill-posed reconstruction by significantly reducing the total number of imaging voxels. Regularization is empirically implemented as λ = p √ σ 1 to stabilize the solution, where p is proportional to the US measured tumor size and σ 1 is the largest eigenvalue of W † W, where W is the weight matrix [16]. However, images reconstructed using these algorithms without regularization in depth may show light shadowing due to the high absorption of the upper portion of the lesion, which may lead to inaccurate quantification of the lower portion [17]. Althobaiti et al. investigated DRI using US gray-scale images, and their method improved the malignant-to-benign lesion contrast as well as the lesion shape and distribution [18]. However, DRI using US gray-scale images is limited because the contrast of US images is poor due to speckle noise and posterior acoustic shadowing of larger tumors. In recent years, deep learning-based image segmentation has been widely used in many areas, and several segmentation neural networks have been introduced in medical imaging [19][20][21]. In this study, we have implemented a convolutional neural network (CNN) to automatically segment each US image into lesion and background regions [21]. The 2D target shape in one spatial dimension and in depth is automatically extracted, and the target symmetry in another spatial dimension is used to obtain an approximate 3D target shape. From the 3D target shape, coarse and fine meshes are defined and the regularization matrix in depth for image reconstruction is constructed. This approach allows near real-time image segmentation and speeds up image reconstruction, bringing US-guided DOT one step closer to clinical use.
Recently, sparse regularization has been shown to improve DOT image reconstruction because its sparsity is robust against noise and enables the preservation of edges in images. In this study, fast iterative shrinkage-thresholding (FISTA) [22] is used with the L1 regularization matrix to constrain the reconstruction [23][24]26]. We call our new approach with depthregularization "depth-regularized reconstruction". The resulting absorption maps have shown more homogeneous absorption distribution in depth. The performance of our depth-regularized reconstruction has been evaluated using simulated targets, phantom targets, and clinical data, all compared with the results of a regularized reconstruction without depth-dependent regularization [17].

Diffuse optical tomography
Photon migration inside a highly scattering medium, such as breast tissue, can be modeled by a diffusion equation of the photon-density wave [25]. Assuming the optical properties inside a small voxel are constant and the optical properties change smoothly across all the voxels, the diffusion equation can be reformulated as a Helmholtz wave equation. Then, by linearizing the problem based on Born approximation and discretizing the imaging space into N disjoint voxels, the equation becomes where U SC represents the perturbation measurements and M is the number of measurements, which is the number of sources multiplied by the number of detectors. W is the weight matrix computed from the analytic solution of the diffusion equation for a semi-infinite medium, using the optical properties of the background breast tissue, and δµ a is the absorption distribution to be reconstructed. We solve the inverse problem of Eq. (1) with a regularization term R(δµ a ) = ||diag(λi)δµ a || 1 1 : Ultrasound has been introduced to provide prior information of the lesion's location and size. With co-registered US images, a dual-mesh method segments the medium into a fine-mesh ellipsoid lesion region and a coarse mesh background tissue region [15]. Mathematically, we separate W and δµ a into two categories, W = [W f , W C ] and δµ a = [δµ af , δµ aC ], where W f and δµ af are the weight matrix and the changes of absorption distribution inside the lesion region, and W C and δµ aC are the weight matrix and changes of absorption distribution in the background region. This simple ellipsoid dual-mesh schedule is used for both σ 1 -regularized reconstruction (defined below) and depth-regularized reconstruction.

σ 1 -regularized reconstruction using FISTA
The inverse problem in solving absorption distributions is described in Eq. (2), where ||U SC − Wδµ a || 2 is used to measure how well the data fit the model, R(δµ a ) = ||diag(λ)δµ a || 1 1 is used to regularize the reconstruction, where the l 1 norm is used in R(δµ a ) to obtain a sparse solution to improve the accuracy of reconstruction. Regularization is implemented as λ = p √ σ 1 , where p is proportional to the US measured tumor size and σ 1 is the largest eigenvalue of W H W [16].
We refer to this algorithm as σ 1 -regularized reconstruction in the rest of the manuscript. FISTA with a constant step size was used to solve the inverse problem [22,26]. The proximal gradient method solves the optimization problem with a gradient step followed by a proximal step. The reconstruction method of the inverse problem is encapsulated in Algorithm 1 given below. We use a zero initialization for the absorption coefficient distribution δµ a = [0] N×1 . The intermediate variables s and q are also initialized accordingly, as described in step 2 [22]. We then go through the iteration. The gradient of the ||U SC − Wδµ a || 2 used in step 4 can be computed as where W H is the Hermitian adjoint of the weight matrix W. A step size of τ is chosen as 2/norm(W H W) in this paper. In step 5 of the algorithm, we compute δµ a t using the proximal operator associated with the regularization term R(δµ a ), where the proximal operator is defined as This proximal operator can be efficiently calculated with the soft threshold function S γ (z), defined as [22,26] where is the element-wise multiplication operator and γ = τλ. Here, the sign(z) is a sign function that extracts the sign of each element in z, and |z| computes the absolute value of each element in z. We then update the intermediate variables s t and q t , following the procedures listed as step 6 and step 7, which help accelerate the convergence rate [22].

Ultrasound segmentation using CNN and depth-regularized reconstruction
Co-registered US images can help to localize the lesion region and reduce the number of voxels with unknown optical properties, easing the ill-posed problem. However, due to the low contrast, speckle noise, and posterior shadow in most breast US images, it remains challenging to accurately perform automated lesion segmentation [20][21]. Typically, an ellipsoid is used to approximate the lesion, and more than twice the lesion size in the spatial dimension is used to define the dual-mesh region, which accounts for larger target region in low resolution DOT imaging. The depth is manually measured from co-registered US images. Here we take a further step and extract lesion shape information from US images as a prior for depth-regularized DOT reconstruction.
, where the index i is the i-th power of the absolute depth at i-th target layer. For example, index i is 1 for the top target layer in the dual-mesh, 2 for the second layer, … and i for the i-th target layer. Here, ℎ( ) is the width of the i-th target layer automatically measured by the CNN. c is a constant empirically chosen as 4 based on simulations. The units of width and depth are in cm.
outside the dual-mesh region is given as zeros. Thus, the regularization is tighter for shallow target layers and looser for deeper layers. Algorithm 1. Fast iterative shrinkage-threshold algorithm with constant step size .

t = t + 1 9. End For
The automated convolutional neural network (CNN) for breast cancer detection was developed by Yap et al. [20]. Their CNN has 33 convolution layers. The network architecture is one convolution layer with a 7×7 kernel, one pooling layer, 6 convolution layers with 3×3 kernels, one pooling layer, 8 convolution layers with 3×3 kernels, one pooling layer, 12 convolution layers with 3×3 kernels, one pooling layer, and 6 convolution layers with 3×3 kernels. Instead of a regular convolution, the last two stages use dilated convolution (also known as atrous convolution) with a dilation rate of 2. A pyramid scene parsing network [21] extracts the feature map from the output of the final convolution layer. The CNN model is pre-trained on ImageNet, which is an image database with 1.2 million images. The input RGB image has three channels, for red, green, and blue in general. We fine-tuned the adapted CNN by freezing the weights of the first two convolutional layers, which captures universal features like boundaries and edges, then we trained the model using PASCAL, a well-known standard dataset with US images [21], and a small learning rate to modify the parameters to fit our task.
The segmentation of US images using the pre-trained CNN was performed automatically by selecting four markers, two each in the lateral and depth dimensions. Those four marks perform as boundary points. We introduced an extra channel in the image with four 2D Gaussian functions centered at each marker. After selecting the four markers, this CNN automatically generated the 2D shape of the lesion. The target symmetry in another spatial dimension was used to obtain an approximate 3D target shape. By performing this semi-automatic interactive segmentation, we could efficiently localize the lesion and extract its approximate shape, including its width at each depth layer and depth as a prior, and could incorporate these parameters into the depth-regularized DOT reconstruction with R(δµ a ) = ||diag(λi)δµ a || 1 1 , as given below. The depth-dependent regularization matrix [diag(λi)] is computed as follows. At target layer i within the dual-mesh, which is generally 2 to 3 times the spatial dimension measured by co-registered US using CNN, the regularization matrix is computed as λi = C width(i) × depth i , where the index i is the i-th power of the absolute depth at i-th target layer. For example, index i is 1 for the top target layer in the dual-mesh, 2 for the second layer, . . . and i for the i-th target layer. Here, width(i) is the width of the i-th target layer automatically measured by the CNN. c is a constant empirically chosen as 4 based on simulations. The units of width and depth are in cm. λi outside the dual-mesh region is given as zeros. Thus, the regularization is tighter for shallow target layers and looser for deeper layers.
To quantify the light shadowing effect, the ratio of the summation of the reconstructed δµ a top from the top depth layer to the summation of δµ a bottom from the bottom layers is calculated as R = δµ a top layer δµ a bottom layers .

Phantom and patient experiments
To evaluate the performance of the proposed algorithm, phantom and patient data were acquired from the US-guided DOT system. The DOT system uses four laser diodes with wavelengths of 740, 780, 808, and 830 nm, and it incorporates 10 or 14 parallel photomultiplier (PMT) detectors and a commercial US system [6,27]. Each laser diode is modulated at 140 MHz and delivers light sequentially to nine source locations on a hand-held probe. The sources are placed on one side of the ultrasound transducer and detectors are placed on the other side, with source-to-detector separations varying from 3 cm to 7.6 cm. The probe diameter is 9 cm, and the volume underneath the probe is represented by voxels. The patient study protocol for this study was approved by the local Institutional Review Boards and was compliant with the Health Insurance Portability and Accountability Act. All patients signed an informed consent form, and all patient data used in this study were de-identified.
A total of 26 patients with large lesions were studied to evaluate the algorithms. Based on biopsy results, 12 patients had malignant lesions (mean age 61 years; range 47-73 years) and 14 had benign lesions (mean age, 44 years; range 22-71 years). For malignant lesions, the average size was 1.8 ± 0.63 cm and average depth was 2.2 ± 0.39 cm; while for benign lesions, the average size was 2.2 ± 0.92 cm and average depth was 1.6 ± 0.43 cm.

Ultrasound segmentation
Ultrasound segmentation, the first step in the proposed algorithm, provides the lesion location and shape information. With accurate segmentation, we can then achieve more accurate reconstruction results. The pre-trained CNN model is applied to our co-registered US images to perform US segmentation.
To evaluate the performance of the US segmentation, an experienced US imager delineated the boundaries of all breast lesions studied. One example is presented in Fig. 1. Based on the evaluations by US experts, the automated interactive neural network gives very similar results to the manual measurements.

Monte Carlo simulation
To evaluate the performance of the proposed algorithm, MC simulations were performed using the known optical properties of the target and background tissue. In the MC simulations, photons propagated in a 10 cm × 10 cm × 6 cm volume whose background optical properties were µ a 0 = 0.02 cm −1 and µ sp 0 = 7 cm −1 . Three target shapes were simulated using these same target optical properties of µ a = 0.20 cm −1 and µ sp = 7 cm −1 . For the first shape, the target was set as two cylinders of different diameters, stacked symmetrically with their centers aligned, to mimic a shaped lesion. The top cylinder was 1.5 cm in diameter and 0.75 cm high, and the bottom one was 3 cm in diameter and 0.75 cm high. The center of this stacked target was located at (0, 0, 2 cm) inside the background volume. Together, the stacked cylinders measured 1.5 cm in height. Within the image background volume, the interfacial plane between the cylinders was set at a depth of 2 cm. Imaging two planes at 2.5 mm above and below that depth, as is typical practice in our DOT, would yield reconstructions at depths of 1.75 and 2.25 cm. Two sets of reconstructed images were generated by using the σ 1 -regularized reconstruction and the proposed depth-regularized reconstruction. For the second target shape, both the top and bottom cylinders were 1.5 cm in diameter and 0.75 cm high, while for the third shape, the larger cylinder was placed on top and the smaller one below.
Longitudinal sections of simulated target and reconstructed absorption maps from both algorithms are presented in Fig. 2. For the first shape, the σ 1 -regularized reconstruction has a maximum reconstructed value of µ a = 0.199 cm −1 in (d), but only the top portion of the lesion is resolved. The reconstructed FWHM width is 3 cm, which is two times larger than the true target dimension of 1.5 cm. The depth-regularized reconstruction provides a similar maximum reconstructed value of µ a = 0.197 cm −1 , along with the second reconstructed depth layer with 2.25-cm FWHM, which is closer to the real target shape shown in (e). Similar results, (f) and (g), are shown for the second shape, where the σ 1 -regularized reconstruction has a maximum reconstructed value of µ a = 0.170 cm −1 , but fails to resolve the target in the second depth layer. The depth-regularized reconstruction provides a maximum reconstructed value of µ a = 0.171 cm −1 along with a second reconstructed target layer of R = 1.65. The second reconstructed target layer has a similar reconstructed FWHM width to the first layer's target width, which is closer to the real target size. However, for the third target shape, shown in (h) and (i), due to the large absorbing top cylinder, both algorithms failed to resolve the target in the second depth layer.
The importance of accurate US measurement in resolving the target depth distribution was also evaluated by simulations. A simple ellipsoid fitting was applied to all three simulated target shapes to acquire the approximate width of each target at different depths. Two sets of images were reconstructed for each shape, based on the depth-regulated algorithm with ellipsoidal fitting and the same algorithm with the actual target width. (f) and (g) in Fig. 3 show that the simple ellipsoid fitting works fine when the target width is closer to the measurement from the ellipsoid. However, (d) and (e) show the ellipsoid fitting leads to a larger ratio, R = 2.45, as compared with using a more accurate measurement of the width, which yields R = 1.32. For target shape (c), due to the highly absorbing target on the top, both methods failed to resolve the bottom portion.

Phantom experiments
Phantom experiments were conducted at a wavelength of 780 nm to evaluate the performance of the depth-regularized algorithm. Solid ball phantoms made of polyester resin were used to mimic a high contrast malignant lesion with calibrated optical properties of µ a = 0.23 cm −1 and µ sp = 7 cm −1 . These target phantoms were submerged 2 cm deep in background intralipid solution with µ a 0 in the range of 0.02 to 0.03 cm −1 and µ sp 0 in the range of 7 to 8 cm −1 , typical bulk optical properties for normal breast tissue. An example of a 2 cm diameter target is shown in Fig. 4. Figure 4(a) is a co-registered US image showing the spherical phantom in the center of the image. Absorption distribution maps (b) reconstructed from the σ 1 -regularized algorithm have a maximum µ a = 0.214 cm −1 , but resolve the target from only the 1.5 cm depth. (c) The depth-regularized algorithm has a similar maximum µ a = 0.204 cm −1 , and resolves the target in both depth layers. However, it cannot recover the third target layer at 2.5 cm. Reconstructed images from the σ 1 -regularized algorithm, with a maximum reconstructed µ a = 0.214 cm −1 . Only the target top-layer was resolved in depth, which leads to an infinite top-to-bottom ratio. (c) Reconstructed images from the proposed depth-regularized algorithm, with a maximum reconstructed µ a = 0.204cm −1 . The proposed algorithm resolves both target layers, which is closer to the phantom shape and has a R = 1.071 top-to-bottom ratio. Note that each image slice is 9 cm by 9 cm. The depth spacing between layers is 0.5 cm, and the absorption coefficient is measured in cm −1 . Note that the target is marked by a circle, and some sound reverberations from the hard resin ball are shown outside the circle.
Three high contrast phantoms, with µ a = 0.23 cm −1 , and three low contrast phantoms, with µ a = 0.11 cm −1 , were reconstructed using both algorithms. The σ 1 -regularized algorithm failed to resolve the second target layer in the high contrast cases and resolved two out of three low contrast cases of R = 1.948 ± 0.276. The low contrast case in which the σ 1 -regularized algorithm could not resolve the second target layer was located at 3 cm depth. The proposed depth-regularized algorithm provides an average ratio R = 1.022 ± 0.155 for the low contrast phantoms and an average ratio R = 1.269 ± 0.280 for the high contrast phantoms.

Clinical study
To evaluate the performance of the proposed algorithm, clinical data were acquired from both the lesion side of a breast, as target data, and the contralateral position of the normal breast, as reference data. Background breast tissue optical properties were computed by fitting the reference data. Four wavelength absorption distribution maps were reconstructed for each case, then used to compute hemoglobin concentrations. Figure 5 shows example reconstructed 780 nm absorption distribution maps of a stage 2 invasive cancer with mixed ductal and lobular features. A co-registered US image is shown in Fig. 5(a), and the resulting image after applying the CNN, with the lesion shape marked, is shown in Fig. 5(b). Images (c) reconstructed with the σ 1 -regularized algorithm show that the top portion of the lesion in slice 4 (2 cm depth) has a much higher µ a of 0.307 cm −1 and is much larger than the bottom portion in slice 5 (2.5 cm depth). The depth-regularized algorithm (d) improves the reconstruction for the second layer, with R = 1.577. One example of a benign fibroadenoma, reconstructed from imaging at 780 nm, is presented in Fig. 6 to demonstrate the performance of the proposed algorithm. Image (c), reconstructed with the σ 1 -regularized algorithm, has an R = 2.556. The depth-regularized algorithm (d) improves reconstruction for the second layer (slice 3), with R = 1.085. Finally, after applying both the σ 1 -regularized algorithm and the depth-regularized algorithm to all the patient data, the total hemoglobin concentrations were computed from the reconstructed absorption distributions of all four wavelengths. The boxplot in Fig. 7 shows the maximum total hemoglobin concentrations (HbT) for all 26 cases for both the σ 1 -regularized algorithm and depth-regularized algorithm. For the benign cases, the total hemoglobin from the σ 1 -regularized algorithm is 63.44 ± 26.16µM, and for the depth-regularized algorithm, it is 50.32 ± 24.19 µM. For the malignant cases, the value is 122.76 ± 26.29 µM for the σ 1 -regularized algorithm, and 117.18 ± 28.79 µM for the depth-regularized algorithm. From the reconstructed values, we can conclude that the depth-regularized reconstruction algorithm has similar performance for both absorption coefficients and total hemoglobin, as compared to the σ 1 -regularized algorithm. However, the depth-regularized algorithm provides improved depth distribution.

Quantitative analysis of the light shadowing effect
The top-to-bottom ratio of the light shadowing effect was quantitatively analyzed for 26 cases of patient data. Due light shadowing, some clinical cases do not have reconstructed values for their bottom portions when the σ 1 -regularized algorithm is used, and so have an infinite ratio for their quantitative measurements of light shadowing effect. After excluding those cases, we  reconstructed the absorption distribution maps of 12 benign and 11 malignant cases, and then computed the ratios, R, of light shadowing effect. From the boxplots in Fig. 8, we can see clear differences between the σ 1 -regularized and the depth-regularized algorithms. The proposed depth-regularized algorithm has a small ratio, which means more homogeneously reconstructed absorption distribution maps in depth. Statistically, the σ 1 -regularized algorithm has mean ratios of 3.53 in benign cases and 3.82 in malignant cases. These values are expected, because in benign cases the lesions are less absorbing than in malignant cases, which leads to reduced light shadowing. The depth-regularized algorithm reconstructs absorption distributions with reduced light shadowing effect, with a mean ratio of 1.19 for benign cases and 1.59 for malignant cases, corresponding to 3.0 to 2.4 times improvements in top-to-bottom reconstructed homogeneity for benign and malignant large lesions.

Evaluation on small targets
Although the proposed algorithm focused on large lesions, we want to show that the depthregularized algorithm will not change reconstruction results with small targets. We evaluated the proposed algorithm's performance on 5 benign and 5 malignant cases with small targets, where the lesion diameters were smaller than 1 cm. One example of a small 6 mm ductal carcinoma in situ is presented in Fig. 9, with the lesion located 1.3 cm beneath the skin. The σ 1 -regularization algorithm and depth-regularized algorithm provide both similar results for reconstructed values and reconstructed image shape.
For five benign lesions, the average size was 1.0 ± 0.33 cm and average depth was 1.33 ± 0.38 cm; while for five malignant lesions, the average size was 1.0 ± 0.44 cm and average depth was 1.6 ± 0.44 cm. For the benign cases, the total hemoglobin of the σ 1 -regularized algorithm is 47.16 ± 12.27 µM, and for the depth-regularized algorithm it is 46.27 ± 14.97 µM. For the malignant cases, the total hemoglobin of the σ 1 -regularized algorithm is 137.73 ± 21.93 µM, and for the depth-regularized algorithm it is 138.53 ± 21.37 µM.

Robustness test
The depth-regularized algorithm has shown promising results in reducing the light shadowing effect and improving target depth distribution. The key is the construction of the depth-dependent regularization matrix, which relies on the US segmentation. Due to the low contrast and speckle noise of US images and the user dependency of our semi-automated segmentation algorithm, different users may realize different segmentation results. Here we evaluated the influence of segmentation errors in reconstruction of a simulated target, where the target shape was known. The simulated target was a 2 cm diameter sphere with an absorption coefficient of µ a = 0.2 cm −1 , located 2 cm below the surface of the tissue. Three target layers of 1.5 cm, 2 cm, and 2.5 cm were reconstructed. One set of images, reconstructed using the ground truth target shape, has a maximum reconstructed absorption coefficient of 0.184 cm −1 . The segmentation error of each layer was defined as the percentage error of the width of the segmented region in the corresponding reconstructed layer. The error is positive when the segmented region is larger than the true size of the lesion and negative otherwise. We manually added 5%, 10%, 15%, and 20% positive and negative errors to each individual reconstructed layer and to all reconstructed layers simultaneously. Using the proposed depth-regularized algorithm, 24 sets of images were reconstructed. The reconstructed maximum absorption coefficients, versus segmentation errors, are presented in Fig. 10. The largest error differs by only 11% from the original reconstructed absorption coefficients, which means the proposed algorithm is fairly robust against small errors in segmentation. It is interesting to note that when the error of top target portion is positive, the width is larger, the λ1 is smaller, and the regularization is looser. Therefore the reconstructed µ a is larger than the use of the true width at the same depth layer. When the error is negative, the trend is reversed. However, for the second target layer, the positive error causes looser regularization and therefore a higher reconstructed µ a value in the second target layer. Since the maximum is always located in the top target layer, this reduces the maximum µ a . The negative error in the second target layer increases the value of the reconstructed maximum target µ a . The overall error follows the trend of the first target layer errors. The highest error is 11%, when the top layer is measured as 20% smaller than the true width, and the error is about 5% when the top layer is measured as 20% larger than the true width at the same depth. The second layer error is about 5% when the second layer is measured as either smaller or larger than the true width at the same depth. These findings suggest that when the uncertainty is higher in segmentation, it is the best to use a larger width for the top target layer. This could be done by providing guidelines to operators.

Discussion and summary
As presented here, deep learning-based US segmentation extracts the lesion location and the lesion shape from co-registered US images as a prior, which is incorporated into a regularization matrix in the depth-regularized reconstruction algorithm. The choice of regularization matrix is both important and difficult. With the help of an accurate and efficient US segmentation neural network, a lesion's width and its reconstruction depth can both be used for constructing a depth-dependent regularization matrix. This adaptive regularization matrix gives us a better estimate of the true target absorption distribution in depth.
The depth-regularized reconstruction does not change the maximum reconstructed target absorption coefficient, and therefore does not change the maximum total hemoglobin concentration when it is compared with the σ 1 -regularized algorithm. Thus, it does not change the quantification of the malignant-to-benign lesion contrast ratio. However, in this work, it improved the lesion depth distribution more in the benign group, with a ratio of 3, than in the malignant group, with a ratio of 2.4. Thus, the improved depth distribution could be used as another parameter to distinguish benign vs. malignant lesions, i.e, benign lesions have a more homogenous depth distribution than malignant lesions. However, this initial result will need to be further tested on a large patient population.
Depth-regularized reconstruction is limited in recovering the depth distribution if the top of the target mass is much larger than the bottom of the mass, as seen from target shape 3 in the simulations (Figs. 2 and 3). Additionally, the algorithm cannot resolve the deep portion of the target beyond the second layer, as seen in the Fig. 4 phantom data, which is more than one centimeter below the top target layer in our reconstruction mesh set-up. In general, the target depth sensitivity is about 3 cm when the center of the target is located at this depth. This depth is adequate to image almost all breast lesions. Note that we can always position a patient differently by using a wedge underneath the patient in a supine position to make sure the lesion is within the depth range.
The automated CNN segmentation algorithm provides the lesion shape to extract the spatial and depth priors and reduce the need for manually measuring these parameters. It will also minimize user dependence in the reconstruction of lesion absorption maps, which is an important step toward the clinical translation of US-guided DOT technology. However, the CNN algorithm still requires users to input the boundary points, and the user must have basic knowledge of US images. This task can be readily accomplished by ultrasound sonographers who are trained to assist radiologists in breast exams. Ideally, the 3D shape of the target gained from several co-registered 2D US images can refine the regularization matrix in depth and further improve the accuracy of the reconstructed absorption distribution. This approach will require scanning the lesion in another dimension using a linear array. In future work, we will evaluate the additional improvement and its clinical utility.
In conclusion, we demonstrate a depth-regularized reconstruction algorithm combined with a semi-automated interactive neural network to reconstruct the absorption distribution and total hemoglobin concentration. Co-registered US images are segmented as both spatial and depth priors and are incorporated into the reconstruction. The depth-regularized algorithm is validated using both simulations and phantom and clinical data. The light shadowing effect is reduced by 2-3 times in large lesions. The algorithm improved depth distribution more in a benign group than in a malignant group, a difference that could provide another parameter for breast cancer classification. We believe our work can be generalized for US-guided diffuse optical tomography/fluorescence tomography in reflection geometry.