Improving DOT reconstruction with a Born iterative method and US-guided sparse regularization

: Ultrasound (US)-guided diﬀuse optical tomography (DOT) is a promising low-cost imaging technique for diagnosis and assessment of breast cancer. US-guided DOT is best implemented in reﬂection geometry, which can be co-registered with US pulse-echo imaging and also minimizes the tissue depth for adequate light penetration. However, due to intense light scattering, the DOT reconstruction problem is ill-posed. In this communication, we describe a new non-linear Born iterative reconstruction method with US-guided depth-dependent (cid:96) 1 sparse regularization for improving DOT reconstruction by incorporating a priori lesion depth and shape information from the co-registered US image. Our method iteratively solves the inverse problem by updating the photon-density wave using the ﬁnite diﬀerence method, computing the weight matrix based on Born approximation, and reconstructing the absorption map using the fast iterative shrinkage-thresholding optimization algorithm (FISTA). We validate our method using both phantom and patient data and compare the results with those using the ﬁrst order linear Born method. Phantom experiments demonstrate that the non-linear Born method provides more accurate target absorption reconstruction and better resolution than the linear Born method. Clinical studies on 20 patients show that non-linear Born reconstructs more realistic tumor shapes than linear Born, and improves the malignant-to-benign lesion contrast ratio from 2 . 73 to 3 . 07, which is a 12 . 5% improvement. For lesions approximately more than 2 . 0 cm in diameter, the average malignant-to-benign lesion contrast ratio is increased from 2 . 68 to 3 . 31,

breast cancers and assessment of breast cancer neoadjuvant chemotherapy response, where a priori lesion location information provided by coregistered US images is used to assist DOT reconstruction [4,[11][12][13]. US-guided DOT is best implemented in reflection geometry, which is compatible with conventional US pulse-echo imaging. Additionally, the light propagation depth is reduced to several centimeters when the patient is in a supine position, favoring reflection mode imaging.
In breast DOT reconstruction, linear Born approximation is widely used to compute the weight matrix for modeling the forward problem. With background optical properties estimated from the contralateral breast or a homogeneous medium, the weight matrix can be computed analytically and optimization methods can be used to search for the lesion absorption distribution iteratively [14,15]. Unfortunately, linear Born approximation tends to underestimate the absorption coefficients when the lesion is large and highly absorbing. To solve this problem, Yao et al. proposed an Born iterative method for more accurate photon-density wave estimation [16]. They showed, through simulation, that the non-linear Born iterative method could achieve more accurate absorption reconstruction results in cases where first order linear Born failed.
Sparse regularization is a newer approach for improving DOT reconstruction accuracy and robustness. Because of the intense light scattering, the DOT image reconstruction problem is usually ill-posed [14]. In the past, Correia et al. developed a DOT algorithm based on edge-preserving total variation (TV) regularization [17]. Lee et al. proposed forming the imaging problem as a joint sparsity recovery problem [18]. In recent years, 1 regularization has gained popularity in solving image reconstruction problems. It has been shown, both analytically and experimentally, that incorporating a priori structure information using 1 regularization can reconstruct high quality images, even when insufficient measurements are provided [19,20].
In this article, we present a new depth-dependent 1 regularized non-linear Born iterative reconstruction method for US-guided DOT. Our approach iteratively updates the photondensity wave using a finite difference method and computes the weight matrix based on Born approximation. We validate our method using phantom and patient data and compare the results with those obtained using the first order linear Born method. Phantom experiments demonstrate that the proposed method reconstructs the absorption distribution more accurately than linear Born. Clinical studies of 20 patients have shown that our method provides more realistic tumor shapes than linear Born, and improves the average malignant-to-benign lesion contrast ratio from 2.73 to 3.07, which is a 12.5% improvement. For lesions approximately more than 2.0 cm in diameter, the average malignant-to-benign lesion contrast ratio is increased from 2.68 to 3.31, which is a 23.5% improvement.

Theory and methods
In the following subsections, we discuss how to solve each individual sub-problem and combine them together. In Section 2.1, we review how to model light propagation inside breast tissue. In Section 2.2, we explain how to formulate and solve the inverse problem. In Section 2.3, we discuss how to use the estimated optical properties to correct the photon-density wave estimation. In Section 2.4 and 2.5, we describe how the experimental phantom and patient data were acquired using our DOT system. Finally, Section 2.6 summarizes the steps of the proposed algorithm for DOT reconstruction.

Dual-grid Born approximation
NIR photon migration in breast tissue can be modeled by the frequency-domain diffusion equation of the photon-density wave [21]. Assuming optical properties change smoothly inside the breast tissue, we can approximate the frequency-domain diffusion equation as a Holmholtz equation [22], where U(r) is the photon-density wave and S(r) is the source distribution normalized by the diffusion coefficient. The wavenumber is given as where µ a (r) and D(r) = 1/ 3 µ a (r) + µ s (r) are the absorption and diffusion coefficient distributions, respectively. Further, µ s (r) is the reduced scattering coefficient distribution, v is the speed of light in tissue, and ω is the modulation frequency of the photon-density wave. Let U(r, r s ) = U 0 (r, r s ) + U SC (r, r s ), where U 0 (r, r s ) represents the photon-density wave in a homogeneous medium with constant background optical properties µ a0 and µ s0 , and U SC (r, r s ) represents the perturbed photon-density wave due to heterogeneities. U SC (r, r s ) can be written in an integral form as [22] U SC (r, r s ) = in which G(·) is the Green's function satisfying the extrapolated boundary condition [23].
Since there are no analytical solutions for the photon-density wave distribution U(r, r s ) in an inhomogeneous medium, we choose to approximate it with a numerical method, which will be discussed in Section 2.3. We further assume the reduced scattering coefficient varies slowly inside the human breast, and is significantly higher than the absorption coefficient; hence we can approximate the diffusion coefficient distribution with the diffusion constant D(r) ≈ D 0 = 1/ 3 µ a0 (r) + µ s0 (r) [24]. Moreover, k 0 and O(r) can be written as In US-guided DOT reconstruction, to improving the ill-posed inverse problem, we have used a dual grid schedule, discretizing the lesion volume with a fine grid and the background bolume with a coarse grid. Thus, the integration in Eq. (3) can be formulated in matrix form as where W ∈ C M×N = W L , W B , x ∈ R N = x L ; x B , y ∈ C M is the measurement, and is the measurement noise. [· , ·] and [· ; ·] are the horizontal and vertical matrix concatenations, respectively. Also, and x L ∈ R N L and x B ∈ R N B are vectors representing the absorption coefficient distributions of the lesion and the background volume, respectively. Further, and are weight matrices for the lesion and the background volume, respectively.

Gradient-based reconstruction
With the definition of W, x, and y above, we formulate the inverse problem aŝ Algorithm 1 Sparsely regularized DOT reconstruction 1: input: initial guessx 0 = 0, weight matrix W, step size τ, regularization parameters λ 2: set: s 0 ←x 0 , q 0 ← 1 3: for t = 1, 2, . . . do 4: 6: t ← t + 1 9: end for 10: return:x t where y ∈ C M is the measurement, x ∈ R N is the absorption distribution, and λ ∈ R M is a vector of depth-dependent non-negative regularization parameters determined from the lesion height information given by the co-registered ultrasound image. The first term, D(x) = 1 2 y − Wx 2 2 , measures the data fidelity of the forward model, while the 1 -term, R(x) = diag(λ)x 1 , promotes the sparsity level of the reconstructed images. Each element of the vector λ is empirically chosen as 0.01 , where d i is the width of the tumor at the depth of the i th reconstruction slice, measured in centimeters from the co-registered US image. The optimization problem described in Eq. (8) is solved with the fast iterative shrinkage-thresholding algorithm (FISTA) [25]. FISTA is a widely used proximal gradient method with Nesterov's acceleration [26]. The proximal gradient method solves the optimization problem with a gradient step followed by a proximal step. We choose to use FISTA because it is relatively economical to compute and gives the fastest convergence rate among first order optimization methods [26].
The reconstruction method of the inverse problem is encapsulated in Algorithm 1. We use a zero initialization for the absorption coefficient distribution x. The intermediate variables s and q are also initialized accordingly, as described in step 2. We then go through the iteration. The gradient of D(x) in step 4 can be calculated as and τ is the step size for the proximal gradient method, where W H is the Hermitian adjoint of W.
For the experiments reported in this article, we use τ = 1/norm(W H W). In step 5, after updating the intermediate variable z t , we constrainx t using the proximal operator associated with the convex regularization term R(x), defined as This proximal operator can be efficiently calculated with the soft threshold function S γ (z), defined as [25] S where is the element-wise multiplication operator and γ = τλ. Here, sgn(·) is the sign function that extracts the sign of a real number, and |·| calculates 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 [26].

Finite difference photon-density wave estimation
The photon-density wave required for formulating the Born weighting matrices in Eq. (3) can be estimated with the finite difference method [16,27]. If we use Cartesian coordinates to discretize the volume and approximate the differential operations with Frechet derivatives, Eq. (1) can be numerically written as where i, j, and k are indices along the x, y, and z directions; U i, j,k is the discrete sample of the photon-density wave U(r) at the < i, j, k > position; and S i, j,k is the source distribution. By defining u and s as vector representations of the distributions of photon-density wave U and the source s, respectively, Eq. (12) can be implemented as a linear operation on the photon-density wave u, where L is the system matrix, which depends on the optical property distribution of the medium. We can calculate the initial photon-density wave distribution with estimated background optical properties from the contralateral breast, using a fitting method described by Zhou et. al. [28] in detail. After constructing the matrix L based on the estimated absorption distribution, we update the photon-density wave distribution using the conjugate gradient method. Because the absorption coefficients in the coarse-grid region are very close to the background, we update the photon-density wave only inside the fine-grid area mentioned earlier in Section 2.1.

Data acquisition
To perform imaging experiments, we use the US-guided DOT system described in [29], which consists of a NIR imaging system and a commercial US system. Four laser diodes at wavelengths of 740, 780, 808, and 830 nm are used to generate NIR light modulated at a 140 MHz carrier frequency. The light is then multiplexed to nine positions and delivered to the medium sequentially. Fourteen parallel photo-multiplier detector (PMT) tubes detect the reflected light from the medium. A custom-made A/D board digitizes the fourteen-channel data.

Experimental procedures
Phantom experiments and patient studies were conducted to evaluate the performance of the proposed algorithm. We used phantom experiments to validate both the accuracy and resolution of the proposed method. We first submerged phantom targets having different optical contrasts in an intralipid solution. The high contrast target was made of material with µ a = 0.23 cm −1 and µ s = 7 cm −1 , while for the low contrast target material, µ a = 0.11 cm −1 and µ s = 7.5 cm −1 . We submerged three spherical targets with diameters of 1, 2, and 3 cm at depths of 0.5, 1.0, 1.5, and 2.0 cm. Note that these depths were measured at the surface of the target using co-registered US images. The intralipid solution had an absorption coefficient µ a 0 = 0.02 − 0.03 cm −1 and a reduced scattering coefficient µ s 0 = 7 − 8 cm −1 , which were acquired by fitting [28]. Then we explored the resolution of the proposed algorithm by submerging two 1.0 cm diameter high contrast (µ a = 0.23 cm −1 ) spherical targets inside the intralipid solution at 1.5 cm depth. The two balls were both placed in the center region along the US B-scan direction. Finally, we tested the performance of the proposed algorithm using data from 20 patients, of which 10 patients had malignant lesions, and 10 patients had benign lesions, based on biopsy results. The study was approved by the local Institutional Review Board (IRB) and was compliant with the Health Insurance Portability and Accountability Act (HIPPA). Informed consent was given by every patient, and the data used in this study have been de-identified. We imaged both the lesion and the normal contralateral breast with our US-guided DOT system. Measurements from the contralateral normal breast were used to estimate the average background optical properties of the breast [30].

Summary of the proposed algorithm
In summary, the proposed algorithm combines depth-dependent sparse regularization with a non-linear Born iterative method. First, we reconstruct the lesion absorption by solving the inverse problem using FISTA. The strength of 1 regularization for each depth is determined by the height of the lesion, measured from the co-registered US images. Then we re-calculate the photon-density wave using the finite-difference method and obtain updated estimations of the weight matrix and the target absorption distribution. The overall process is pictured in the flowchart shown in Fig. 1. Given the perturbed photon-density wave measurement u sc (r), in step one, we initialize the wavenumber distribution k(r) and the photon-density wave distribution u(r) with those of the homogeneous background media. In step two, we use the photon-density wave distribution to form the weight matrix W. In step three, we reconstruct the absorption coefficient distribution O(r), and recalculate the photon-density wave u(r) and wave number k(r) distributions. We then repeat steps two and three until adequate iterations are reached. For the experiments reported in this article, the iteration stops after 10 iterations.

Accuracy experiments
Phantom experiments were performed with the methods described in Section 2.5. We compared our reconstruction results with the first order linear Born method developed by us earlier [11]. Figure 2 shows reconstructed images of a high optical contrast ball phantom located at 1.5 cm (top surface) depth inside the intralipid solution. The 3D absorption distribution is displayed as slices at different depths, labeled above each column: (a-c) and (d-f) are reconstructions of one 2 cm diameter ball using non-linear Born (µ a ma x = 0.22cm −1 ) and linear Born (µ a ma x = 0.18cm −1 ), respectively. A more comprehensive analysis of the accuracy of absorption coefficients is shown in Figure 3. HC and LC stand for high contrast (µ a = 0.23cm −1 ) and low contrast (µ a = 0.11cm −1 ), respectively. S, M, and L stand for small (1 cm diameter), medium (2 cm), and large (3 cm), respectively. The color bars in the right upper legend indicate the depth of the top layer of the phantom target. The average absorption coefficients were estimated with 89.6% accuracy for high contrast phantoms and 86.1% for low contrast phantoms. The accuracy is  calculated as µ a ma x /µ a t r u t h × 100%. Further, we present the convergence analysis of our iterative image reconstruction method using phantom data shown in Fig. 3. To compare our method with the conjugate gradient optimization method for linear Born discussed in [31], we normalized the least squares error (LSE) for each method to the power of the scattered field, y 2 . The mean and standard deviation of least square errors (LSE) for each method are plotted as a function of iterations in Fig. 4, where zero initialization is used for both methods. We observe that FISTA converges faster than the conjugate gradient method. Note that, on average, the objective function converges to a lower value for non-linear modeling, because more accurate estimation of the photon-density wave better fits the perturbed photon-density wave measurement u sc (r) in Eq. (3), reducing the LSE.

Two target resolution test
Additionally, we compared the resolution of reconstruction from non-linear Born with that from linear Born by submerging two 1.0 cm diameter high contrast (µ a = 0.23 cm −1 ) ball shaped targets separated by 2 cm along the US B-scan direction inside the intralipid solution at 1.5 cm depth. The reconstruction results using non-linear Born with regularization and linear Born without regularization are illustrated in Fig. 5. The non-linear Born algorithm with sparse regularization gives a smaller full width at half maximum (FWHM) value [32], which resolves the two targets much better than linear Born. Notice that our method does not require employing of two fine-grid regions, as reported in [33].

Patient study
We compared non-linear Born with linear Born across 20 patients, 10 with benign lesions and 10 with malignant ones. Patient data were acquired from the lesion side of the breast and the contralateral mirror position of the healthy breast. The perturbed photo-density wave was calculated as U lesion − U r e f er ence U r e f er ence ,   where U lesion and U r e f er ence are measurements from the lesion and reference breast, respectively. In the past, we have compared the use of a contralateral mirror position of a lesion breast and a symmetric area of the same lesion breast as a healtht breast reference; however, the contralateral reference is more robust because the tissue curvature and the chest wall depth can be made symmetrical under the real-time assessment of co-registered ultrasound [28]. Figure 6 shows a reconstructed tHb,oxyHb, and deoxyHb map of a medium size malignant lesion. The tHb is calculated from absorption coefficients of four wavelengths, with the extinction coefficients for deoxygenated and oxygenated hemoglobin given in the literature [34]. The co-registered US image indicates that the lesion is centered at 2 cm depth from the surface of the breast. The functional maximum tHb concentration reconstructed with non-linear Born and linear Born are 95.0µM and 84.4µM, respectively. The oxyHb distribution closely follows the tHb distribution, but is more heterogeneous, with slightly periphery enhancement. The deoxyHb distribution is more centered in the tumor core. This type of peripheral oxyHb distribution and core deoxyHb distribution is often seen in larger cancers due to the necrotic tissue in the center and rapid tumor growth at the peripery [13]. Figure 7 shows reconstruction results on a benign lesion, and the co-regesitered US image suggests the lesion is located at 2 cm depth. The maximum tHb concentrations reconstructed with non-linear Born and linear Born are 29.8 µM and 28.7 µM, respectively. It is interesting to note that this benign lesion had higher deoxyHb than oxyHb, but both are low. This benign lesion is diagnosed as a proliferate lesion, which may account for the relatively higher deoxyHb component. Finally, we calculate the tHb,oxyHb, and deoxyHb values across all 20 cases. Figure 8 presents the statistics of the reconstructed functional maximum tHb, oxyHb, and deoxyHb values in box plots. Again, we compare non-linear Born with linear Born. The non-linear Born algorithm improves the average malignant-to-benign lesion contrast ratio from 2.73 to 3.07, which is a 12.5% improvement. 12.4% improvement. For oxyHb and deoxyHb, the non-linear Born algorithm does not improve the average malignant-to-benign lesion ratio than that of linear Born. However, the mean oxyHb of non-linear Born of malignant group is higher than that of the linear Born (p=0.04), where p is the p-value from the t-test. The mean oxyHb of non-linear Born of benign group is statistically the same as the linear Born (p=0. 16). This suggests that non-linear Born statistically improves the linear Born on oxyHb estimate for malignant group. For deoxyHb, non-linear Born improves deoxyHb than linear Born for both malignant and bening groups. These interesting premilnary data will be validated with a large patient pool.

Discussion and summary
To summarize, we have experimentally demonstrated and validated that the proposed method can successfully reconstruct functional images of phantom targets and breast lesions. Phantom experiments confirm that the non-linear Born method yields better resolution and more accurate absorption coefficient distributions than the linear Born method. However, non-linear Born also underestimates the target absorption for a high contrast phantom target located shallower than 1 cm deep, due to the lack of a center source in the probe [35].
In clinical cases, we see that non-linear Born reconstructs higher absorption coefficient value for large malignant cases than the linear Born method. Based on the results from 20 patients' data, the average malignant-to-benign lesion contrast is increased from 2.73, using the linear Born method, to 3.07, which is a 12.5% improvement. For lesions approximately more than 2.0 cm in diameter, the average malignant-to-benign contrast is increased from 2.68 to 3.31, which is a 23.5% improvement. Our method can achieve more faithful results than the linear Born method because the photon-density wave attenuation is calculated more accurately with the iterative update, and the US a priori structure information is incorporated adequately through sparsity-promoting regularization. Moreover, our method also presents more realistic tumor absorption distributions.
In the past, we have attempted to simultaneously reconstruct both absorption and scattering distributions of breast lesions [36]. We had success in phantom data, but these algorithms are not robust for patient data when applied to a large patient database. In the simultaneous reconstruction, the distribution of the lesion diffusion coefficient, D(r) = 1/(3 × µ s (r)), is reconstructed with the distribution of the lesion absorption, µ a (r). However,D(r) is one order of magnitude smaller than µ a (r) and cannot be reconstructed reliably for all patient data. Additionally, simultaneous reconstruction of both absorption and scattering distributions doubles the unknowns in the image reconstruction. Since µ a (r) is directly related to tumor angiogenesis, we have focused on this important parameter in our algorithm development.
Although the reconstruction results are improved by the new approach, several issues remain. First, the Born-type modeling requires a reference medium, and contralateral breast data is still needed. Second, the numerical updating of the photon-density wave distribution slows the reconstruction speed. Designing and implementing a GPU-based fast computation method is needed for clinical translation of this algorithm.
To conclude, the proposed non-linear Born method with US-guided depth regularization significantly improves the reconstructed target shape, accuracy, and resolution. Our method uses a non-linear forward model for better photon-density distribution estimation and a fast converging algorithm for solving the inverse problem, incorporating lesion structure information provided by the US image. Moreover, with selective modifications, the method is also applicable to MRI-or X-ray-guided DOT.