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

Shiqi Xu; K. M. Shihab Uddin; Quing Zhu

doi:10.1364/BOE.10.002528

1. Introduction

Diffuse optical tomography (DOT) is a low-cost, non-invasive functional imaging modality that utilizes near-infrared (NIR) light to image biological tissue. In DOT imaging, tissue is illuminated with diffused light, and the reflected or transmitted light is measured at the tissue surface. The tissue’s optical properties are then estimated with image reconstruction algorithms. In the past decades, DOT has been actively studied for diagnosis of breast cancer and assessment of treatment response. Clinical studies have demonstrated that the hemoglobin (Hb) concentration (oxy-Hb, deoxy-Hb, and total Hb) calculated from DOT reconstructed absorption images is directly related to breast tumor angiogenesis, and it can be used to differentiate cancers from benign lesions as well as to assess breast cancer treatment response [1–5].

However, due to the intense light scattering in breast tissue, lesion localization is quite challenging; therefore, multimodality imaging systems have been developed to assist in reconstructing more informative DOT images [6–10]. Recently, ultrasound-guided DOT has emerged as a promising low-cost adjunct technique to ultrasound (US) for the diagnosis of 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 [11–14]. 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, $l_{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 $l_{1}$ regularization can reconstruct high quality images, even when insufficient measurements are provided [19, 20].

In this article, we present a new depth-dependent $l_{1}$ regularized non-linear Born iterative reconstruction method for US-guided DOT. Our approach iteratively updates the photon-density 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.

2. 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.

2.1. 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],

\nabla^{2} U (r) + k^{2} (r) U (r) = - S (r),

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

k (r) = \sqrt{(\frac{i ω}{v} - μ_{a} (r)) / D (r)},

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_{S C} (r, r_{s})

, where

U_{0} (r, r_{s})

represents the photon-density wave in a homogeneous medium with constant background optical properties

μ_{a 0}

and

μ_{s 0}^{'}

, and

U_{S C} (r, r_{s})

represents the perturbed photon-density wave due to heterogeneities.

U_{S C} (r, r_{s})

can be written in an integral form as [22]

U_{S C} (r, r_{s}) = \int U (r, r_{s}) O (r^{'}) G (r - r^{'}) d r^{'},

in which

G (\cdot)

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) \approx D_{0} = 1 / (3 (μ_{a 0} (r) + μ_{s 0}^{'} (r)))

[24]. Moreover, k₀ and

O (r)

can be written as

k_{0} = \sqrt{(\frac{i ω}{v} - μ_{a 0}) / D_{0}}, O (r) = - \frac{δ μ_{a} (r)}{D_{0}} .

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

y = Wx + ϵ,

where

W \in ℂ^{M \times N} = [W_{L}, W_{B}]

,

x \in ℝ^{N} = [x_{L}; x_{B}]

,

y \in ℂ^{M}

is the measurement, and

ϵ

is the measurement noise.

[\cdot, \cdot]

and

[\cdot; \cdot]

are the horizontal and vertical matrix concatenations, respectively. Also, and

x_{L} \in ℝ^{N_{L}}

and

x_{B} \in ℝ^{N_{B}}

are vectors representing the absorption coefficient distributions of the lesion and the background volume, respectively. Further,

W_{L} = {[- \frac{1}{D_{0}} G (r_{v j}, r_{d i}) U (r_{v j}, r_{s i})]}_{M \times N_{L}}

and

W_{B} = {[- \frac{1}{D_{0}} G (r_{v j}, r_{d i}) U (r_{v j}, r_{s i})]}_{M \times N_{B}}

are weight matrices for the lesion and the background volume, respectively.

2.2. Gradient-based reconstruction

With the definition of W, x, and y above, we formulate the inverse problem as

\hat{x} = \underset{x \in ℝ^{N}}{argmin} {\frac{1}{2} ‖ y - Wx ‖_{2}^{2} + ‖ diag (λ) x ‖_{1}},

where

y \in ℂ^{M}

is the measurement,

x \in ℝ^{N}

is the absorption distribution, and

λ \in ℝ^{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) = \frac{1}{2} ‖ y - Wx ‖_{2}^{2}

, measures the data fidelity of the forward model, while the

l_{1}

-term,

R (x) = ‖ diag (λ) x ‖_{1}

, promotes the sparsity level of the reconstructed images. Each element of the vector

λ

is empirically chosen as

\frac{0.01}{d_{i}^{2}}

, where d_i is the width of the tumor at the depth of the

i^{t h}

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].

Table 1. Sparsely regularized DOT reconstruction

View Table

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

\frac{\partial}{\partial x} D (x) = W^{H} (Wx - y),

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 constrain

{\hat{x}}^{t}

using the proximal operator associated with the convex regularization term

R (x)

, defined as

{prox}_{R} (x) = \underset{z \in ℝ^{N}}{argmin} {R (x) + \frac{1}{2} ‖ x - z ‖_{2}^{2}},

This proximal operator can be efficiently calculated with the soft threshold function $S_{γ} (z)$ , defined as [25]

S_{γ} (z) = sgn (z) ⊙ max (0, | z | - γ),

where

⊙

is the element-wise multiplication operator and

γ = τ λ

. Here, sgn

(\cdot)

is the sign function that extracts the sign of a real number, and

| \cdot |

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].

2.3. 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

\begin{array}{l} \frac{U_{i + 1, j, k} + U_{i - 1, j, k}}{Δ x^{2}} + \frac{U_{i, j + 1, k} + U_{i, j - 1, k}}{Δ y^{2}} + \frac{U_{i, j, k + 1} + U_{i, j, k - 1}}{Δ z^{2}} \\ - (\frac{2}{Δ x^{2}} + \frac{2}{Δ y^{2}} + \frac{2}{Δ z^{2}} - k_{i, j, k}^{2}) U_{i, j, k} = - S_{i, j, k}, \end{array}

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,

Lu = s,

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.

2.4. 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 NIRlight 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.

2.5. 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.

Fig. 1 Flowchart of the proposed algorithm as summarized in Section 2.6

Download Full Size | PDF

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].

2.6. 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 $l_{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_{s c} (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.

Fig. 2 (a) 3D contour plot of the phantom target. (b-g) Reconstructed absorption maps from a high contrast 2.0 cm diameter ball phantom placed at 1.5 cm (top surface) depth. The 3D absorption distribution is displayed as slices at different depthslabeled above each column. (b-d) and (e-g) are reconstructions using linear Born ( $μ_{a_{m a x}} = 0.18 {cm}^{- 1}$ ) and non-linear Born ( $μ_{a_{m a x}} = 0.22 {cm}^{- 1}$ ), respectively. The target depth distribution is reconstructed more accurately using the proposed non-linear Born than with the linear Born. The color bar and the scale bar are used used for all images in this figure.

Download Full Size | PDF

Fig. 3 Accuracy of the absorption coefficients reconstructed using non-linear Born for phantoms submerged at various depths. HC and LC stand for high contrast and low contrast, respectively. S,M, and L stand for small, medium, and large, respectively. Color bars in the legend indicate the depth of the top layer of the target.

Download Full Size | PDF

3. Experimental results

3.1. 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_{m a x}} = 0.22 {cm}^{- 1}$ ) and linear Born ( $μ_{a_{m a x}} = 0.18 {cm}^{- 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.23 {cm}^{- 1}$ ) and low contrast ( $μ_{a} = 0.11 {cm}^{- 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_{m a x}} / μ_{a_{t r u t h}} \times 100 %$ .

Fig. 4 (a) Schematic of the probe used for phantom experiments.(b) Normalized least squares error of the conjugate gradient method for linear Born (CG linear) and FISTA for non-linear Born (FISTA non-linear) using phantom data.

Download Full Size | PDF

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_{s c} (r)$ in Eq. (3), thus reducing the LSE.

Fig. 5 Reconstructed absorption maps of two 1 cm diameter high contrast ball phantoms placed at 1.5 cm depth, using linear Born with regularization ( $μ_{a_{m a x}} = 0.23 {cm}^{- 1}$ ) and non-linear Born ( $μ_{a_{m a x}} = 0.23 {cm}^{- 1}$ ), respectively. Both targets are resolved better using non-linear Born without regularization. The color bar and the scale bar are used for (b)-(c).

Download Full Size | PDF

Fig. 6 Reconstructed Hb map of a stage 3 malignant lesion. The 3D Hb distribution is displayed as slices at different depths, labeled above each column. (a) A co-registered US image. (b) A center slice of the reconstructed tHb distribution at the orthogonal plane. (c)-(e) Reconstructed tHb concentration distributions using linear Born without regularization; maximum tHb = 84.4 µM. (f)-(h) Reconstructed tHb concentration distributions using non-linear Born with regularization; maximum tHb = 95.0 µM. (i)-(k) Reconstructed oxyHb concentration distributions using non-linear Born with regularization; maximum oxyHb = 65.33 µM. (l)-(n) Reconstructed deoxyHb concentration distributions using non-linear Born with regularization; maximum deoxyHb = 47.88 µM. The color bar in (n) is used for (b)-(n) and the scale bar in (n) is used for (c)-(n).

Download Full Size | PDF

Fig. 7 Reconstructed Hb map of a benign but proliferative lesion. The 3D Hb distribution is displayed as slices at different depths, labeled above each column. (a) A co-registered US image. (b) A center slice of the reconstructed tHb distribution at the orthogonal plane. (c)-(e) Reconstructed tHb concentration distributions using linear Born without regularization; maximum tHb = 28.7 µM. (f)-(h) Reconstructed tHb concentration distributions using non-linear Born with regularization; maximum tHb = 29.8 µM. (i)-(k) Reconstructed oxyHb concentration distributions using non-linear Born with regularization; maximum oxyHb = 8.2 µM. (l)-(n) Reconstructed deoxyHb concentration distributions using non-linear Born with regularization; maximum deoxyHb = 25.3 µM. The color bar in (n) is used for (b)-(n) and the scale bar in (n) is used for (c)-(n).

Download Full Size | PDF

Fig. 8 Box plot of maximum tHb,oxyHb, and deoxyHb concentrations obtained from 10 malignant and 10 benign cases using linear Born with regularization and non-linear Born without regularization, respectively. The p-values from a t-test are labeled

Download Full Size | PDF

3.2. 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].

3.3. 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

\frac{U_{l e s i o n} - U_{r e f e r e n c e}}{U_{r e f e r e n c e}},

where

U_{l e s i o n}

and

U_{r e f e r e n c e}

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.

4. 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, . , 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.

Funding

National Institute of Health (R01EB002136, R01CA228047).

Acknowledgments

The authors would like to express our great appreciation to Prof. James Ballard for editing, Shuying Li for figure preparation, and Yifeng Zeng for inspirational discussion.

Disclosures

All authors declare that there are no conflicts of interest.

References

1. B. W. Pogue, S. P. Poplack, T. O. McBride, W. A. Wells, K. S. Osterman, U. L. Osterberg, and K. D. Paulsen, “Quantitative hemoglobin tomography with diffuse near-infrared spectroscopy: pilot results in the breast,” Radiology 218, 261–266 (2001). [CrossRef] [PubMed]

2. B. J. Tromberg, B. W. Pogue, K. D. Paulsen, A. G. Yodh, D. A. Boas, and A. E. Cerussi, “Assessing the future of diffuse optical imaging technologies for breast cancer management,” Med. Phys. 35, 2443–2451 (2008). [CrossRef] [PubMed]

3. J. Culver, R. Choe, M. Holboke, L. Zubkov, T. Durduran, A. Slemp, V. Ntziachristos, B. Chance, and A. Yodh, “Three-dimensional diffuse optical tomography in the parallel plane transmission geometry: Evaluation of a hybrid frequency domain/continuous wave clinical system for breast imaging,” Med. Phys. 30, 235–247 (2003). [CrossRef] [PubMed]

4. Q. Zhu, A. Ricci Jr, P. Hegde, M. Kane, E. Cronin, A. Merkulov, Y. Xu, B. Tavakoli, and S. Tannenbaum, “Assessment of functional differences in malignant and benign breast lesions and improvement of diagnostic accuracy by using us-guided diffuse optical tomography in conjunction with conventional us,” Radiology 280, 387–397 (2016). [CrossRef] [PubMed]

5. S. Jiang, B. W. Pogue, P. A. Kaufman, J. Gui, M. Jermyn, T. E. Frazee, S. P. Poplack, R. DiFlorio-Alexander, W. A. Wells, and K. D. Paulsen, “Predicting breast tumor response to neoadjuvant chemotherapy with diffuse optical spectroscopic tomography prior to treatment,” Clin. Cancer Res. 20, 6006–6015 (2014). [CrossRef] [PubMed]

6. Z. Deng, Y. Lin, J. Zimmermann, and G. Gulsen, “Fully automatic ultrasound guided diffuse optical tomography (us-dot) system for whole breast imaging,” in Biomedical Optics, (Optical Society of America, 2012), pp. BTu3A–6.

7. B. Brooksby, B. W. Pogue, S. Jiang, H. Dehghani, S. Srinivasan, C. Kogel, T. D. Tosteson, J. Weaver, S. P. Poplack, and K. D. Paulsen, “Imaging breast adipose and fibroglandular tissue molecular signatures by using hybrid mri-guided near-infrared spectral tomography,” Proc. Natl. Acad. Sci. 103, 8828–8833 (2006). [CrossRef] [PubMed]

8. V. Ntziachristos, A. Yodh, M. D. Schnall, and B. Chance, “Mri-guided diffuse optical spectroscopy of malignant and benign breast lesions,” Neoplasia 4, 347–354 (2002). [CrossRef] [PubMed]

9. Q. Fang, J. Selb, S. A. Carp, G. Boverman, E. L. Miller, D. H. Brooks, R. H. Moore, D. B. Kopans, and D. A. Boas, “Combined optical and x-ray tomosynthesis breast imaging,” Radiology 258, 89–97 (2011). [CrossRef]

10. V. Krishnaswamy, K. E. Michaelsen, B. W. Pogue, S. P. Poplack, I. Shaw, K. Defrietas, K. Brooks, and K. D. Paulsen, “A digital x-ray tomosynthesis coupled near infrared spectral tomography system for dual-modality breast imaging,” Opt. Express 20, 19125–19136 (2012). [CrossRef] [PubMed]

11. Q. Zhu, M. Huang, N. Chen, K. Zarfos, B. Jagjivan, M. Kane, P. Hedge, and S. H. Kurtzman, “Ultrasound-guided optical tomographic imaging of malignant and benign breast lesions: initial clinical results of 19 cases,” Neoplasia 5, 379–388 (2003). [CrossRef] [PubMed]

12. Q. Zhu, E. B. Cronin, A. A. Currier, H. S. Vine, M. Huang, N. Chen, and C. Xu, “Benign versus malignant breast masses: optical differentiation with us-guided optical imaging reconstruction,” Radiology 237, 57–66 (2005). [CrossRef] [PubMed]

13. Q. Zhu, P. U. Hegde, A. Ricci Jr, M. Kane, E. B. Cronin, Y. Ardeshirpour, C. Xu, A. Aguirre, S. H. Kurtzman, P. J. Deckers, et al., “Early-stage invasive breast cancers: potential role of optical tomography with us localization in assisting diagnosis,” Radiology 256, 367–378 (2010). [CrossRef] [PubMed]

14. S. R. Arridge and J. C. Schotland, “Optical tomography: forward and inverse problems,” Inverse Probl. 25, 123010 (2009). [CrossRef]

15. H. Dehghani, S. Srinivasan, B. W. Pogue, and A. Gibson, “Numerical modelling and image reconstruction in diffuse optical tomography,” Philos. Transactions Royal Soc. A: Math. Phys. Eng. Sci. 367, 3073–3093 (2009). [CrossRef]

16. Y. Yao, Y. Wang, Y. Pei, W. Zhu, and R. L. Barbour, “Frequency-domain optical imaging of absorption and scattering distributions by a born iterative method,” JOSA A 14, 325–342 (1997). [CrossRef] [PubMed]

17. T. Correia, J. Aguirre, A. Sisniega, J. Chamorro-Servent, J. Abascal, J. J. Vaquero, M. Desco, V. Kolehmainen, and S. Arridge, “Split operator method for fluorescence diffuse optical tomography using anisotropic diffusion regularisation with prior anatomical information,” Biomed. Opt. Express 2, 2632–2648 (2011). [CrossRef] [PubMed]

18. O. Lee, J. M. Kim, Y. Bresler, and J. C. Ye, “Compressive diffuse optical tomography: noniterative exact reconstruction using joint sparsity,” IEEE Transactions on Med. Imaging 30, 1129–1142 (2011). [CrossRef]

19. E. J. Candès, J. Romberg, and T. Tao, “Robust uncertainty principles: Exact signal reconstruction from highly incomplete frequency information,” IEEE Transactions on Inf. Theory 52, 489–509 (2006). [CrossRef]

20. M. Lustig, D. Donoho, and J. M. Pauly, “Sparse mri: The application of compressed sensing for rapid mr imaging,” Magn. Reson. Medicine: An Off. J. Int. Soc. for Magn. Reson. Medicine 58, 1182–1195 (2007). [CrossRef]

21. J. B. Fishkin and E. Gratton, “Propagation of photon-density waves in strongly scattering media containing an absorbing semi-infinite plane bounded by a straight edge,” JOSA A 10, 127–140 (1993). [CrossRef]

22. L. V. Wang and H.-i. Wu, Biomedical Optics: Principles and Imaging(John Wiley & Sons, 2012).

23. R. Aronson, “Boundary conditions for diffusion of light,” JOSA A 12, 2532–2539 (1995). [CrossRef] [PubMed]

24. V. G. Peters, D. Wyman, M. Patterson, and G. Frank, “Optical properties of normal and diseased human breast tissues in the visible and near infrared,” Phys. Medicine & Biol. 35, 1317 (1990). [CrossRef]

25. A. Beck and M. Teboulle, “A fast iterative shrinkage-thresholding algorithm for linear inverse problems,” SIAM J. on Imaging Sci. 2, 183–202 (2009). [CrossRef]

26. Y. Nesterov, et al., “Gradient methods for minimizing composite objective function,” (2007).

27. B. W. Pogue, M. S. Patterson, and T. J. Farrell, “Forward and inverse calculations for 3d frequency-domain diffuse optical tomography,” in Optical Tomography, Photon Migration, and Spectroscopy of Tissue and Model Media: Theory, Human Studies, and Instrumentation, vol. 2389 (International Society for Optics and Photonics, 1995), pp. 328–340.

28. F. Zhou, A. Mostafa, and Q. Zhu, “Improving breast cancer diagnosis by reducing chest wall effect in diffuse optical tomography,” J. Biomed. Opt. 22, 036004 (2017). [CrossRef]

29. C. Xu, H. Vavadi, A. Merkulov, H. Li, M. Erfanzadeh, A. Mostafa, Y. Gong, H. Salehi, S. Tannenbaum, and Q. Zhu, “Ultrasound-guided diffuse optical tomography for predicting and monitoring neoadjuvant chemotherapy of breast cancers: recent progress,” Ultrason. Imaging 38, 5–18 (2016). [CrossRef]

30. H. Vavadi and Q. Zhu, “A calibration method for diffuse optical tomography based on extracted target depth and size from ultrasound images,” in Optical Tomography and Spectroscopy, (Optical Society of America, 2018), pp. OF1D–4.

31. K. S. Uddin, A. Mostafa, M. Anastasio, and Q. Zhu, “Two step imaging reconstruction using truncated pseudoinverse as a preliminary estimate in ultrasound guided diffuse optical tomography,” Biomed. Opt. Express 8, 5437–5449 (2017). [CrossRef]

32. H. Dehghani, B. W. Pogue, J. Shudong, B. Brooksby, and K. D. Paulsen, “Three-dimensional optical tomography: resolution in small-object imaging,” Appl. Opt. 42, 3117–3128 (2003). [CrossRef] [PubMed]

33. Y. Xu, C. Xu, and Q. Zhu, “Clustered targets imaged by optical tomography guided by ultrasound,” J. Biomed. Opt. 16, 076018 (2011). [CrossRef] [PubMed]

34. M. Cope, “The application of near infrared spectroscopy to non invasive monitoring of cerebral oxygenation in the newborn infant,” Dep. Med. Phys. Bioeng.342 (1991).

35. A. Siegel, J. Marota, and D. A. Boas, “Design and evaluation of a continuous-wave diffuse optical tomography system,” Opt. Express 4, 287–298 (1999). [CrossRef] [PubMed]

36. B. Tavakoli and Q. Zhu, “Two-step reconstruction method using global optimization and conjugate gradient for ultrasound-guided diffuse optical tomography,” J. Biomed. Opt. 18, 016006 (2013). [CrossRef]

1:	input: initial guess ${\hat{x}}^{0} = 0$ , weight matrix W, step size τ, regularization parameters $λ$
2:	set: $s_{0} \leftarrow {\hat{x}}^{0}$ , $q_{0} \leftarrow 1$
3:	for $t = 1, 2, \dots$ do
4:	$z^{t} \leftarrow s^{t - 1} - τ \nabla D (s^{t - 1})$	$⊳$ (9)
5:	${\hat{x}}^{t} \leftarrow$ prox $_{R} (z^{t}, τ λ)$	$⊳$ (11)
6:	$q_{t} \leftarrow \frac{1}{2} (1 + \sqrt{1 + 4 q_{t - 1}^{2}})$
7:	$s^{t} \leftarrow {\hat{x}}^{t} + ((q_{t - 1} - 1) / q_{t}) ({\hat{x}}^{t} - {\hat{x}}^{t - 1})$
8:	$t \leftarrow t + 1$
9:	end for
10:	return: ${\hat{x}}^{t}$

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

Abstract

1. Introduction

2. Theory and methods

2.1. Dual-grid Born approximation

2.2. Gradient-based reconstruction

2.3. Finite difference photon-density wave estimation

2.4. Data acquisition

2.5. Experimental procedures

2.6. Summary of the proposed algorithm

3. Experimental results

3.1. Accuracy experiments

3.2. Two target resolution test

3.3. Patient study

4. Discussion and summary

Funding

Acknowledgments

Disclosures

References

Cited By

Figures (8)

Tables (1)

Equations (14)

Biomedical Optics Express