2.1.

Reconstruction Algorithm with Depth Correction

In our ultrasound-guided DOT image reconstruction, a dual-zone mesh scheme is used to segment the imaging volume into a lesion region, L, and a background region, B, with fine and coarse voxel sizes, respectively.²⁶ The forward model based on the diffusion equation, linearized by modified born approximation, is written as follows:

In Eq. 1, U _sd is the scattered field measured at each source and detector pair, [M _L, M _B] denotes the total absorption distribution, and [W _L, W _B] denotes the weight matrix of lesion and background regions, respectively. Weight matrix describes the distribution of diffusive wave in the homogenous medium and characterizes the measurement sensitivity to the absorption changes in both spatial dimensions and in depth. For computing the weight matrix, we have used a semi-infinite absorbing boundary condition at the surface and the corresponding analytic solution of the diffusion equation. In addition, the average optical absorption and reduced scattering coefficient of the background medium needed for the weight matrix computation are obtained from the measurements made with the homogenous intralipid solution for phantom experiments and the normal contralateral breast for clinical study. Finally, after computing the weight matrix, inversion is performed with the conjugate gradient iteration searching method to obtain the total absorption map from Eq. 1 and the total distribution was divided by different voxel sizes at the lesion and background regions to obtain the absorption distribution.

A proper scaling factor is required to balance the exponential decay of the weight matrix elements with depth. In general, singular value decomposition of a two-dimensional matrix, is intuitively interpreted as a sequence of three geometrical transformations: a rotation, a scaling, and another rotation. The scaling transformation is a diagonal matrix where the diagonal entries are matrix singular values. In addition, the Eckart–Young theorem states that the least-squares approximation in α dimension of a matrix of rank β can be found by replacing its smallest (α–β) singular values to zeroes.²⁷ It means that the large singular values contain most of the information in a data matrix, and this is principally the basis of a method for compressing the information in a matrix. On the basis of these properties, the maximum singular value of a matrix has the maximum information related to the scaling performance of its elements. As a result, for correcting the depth dependency of the weight matrix, the MSV of the weight elements of each layer is used.

Our modified algorithm for balancing the weight matrix using the MSVs of different depth layers, given the coregistered ultrasound information, is as follows. In the first step, the elements of weight matrix are grouped layer by layer from top to bottom in the lesion and background regions as [TeX:] $[W_{\rm L},W_{\rm B}]\break = [(W_{\rm L}^1,\ldots,W_{\rm L}^n),(W_{\rm B}^1,\ldots,W_{\rm B}^m)]$ $[W_{\mathrm{L}},W_{\mathrm{B}}]=[(W_{\mathrm{L}}^1,...,W_{\mathrm{L}}^n),(W_{\mathrm{B}}^1,...,W_{\mathrm{B}}^m)]$ , where [TeX:] $W_{\rm L}^i,i = 1,\ldots,n$ $W_{\mathrm{L}}^i,i=1,...,n$ represent the weights of voxels of the i’th layer of lesion region with n layers, and [TeX:] $W_{\rm B}^j,j = 1,\ldots,m$ $W_{\mathrm{B}}^j,j=1,...,m$ represents weights of voxels of the j’th layer of background region with m layers. Then a scaling diagonal matrix including the MSVs of the layers is composed as [TeX:] $[S_{\rm L},S_{\rm B}] = [{\rm diag}(\lambda _{\rm L}^n,\ldots,\lambda _{\rm L}^1),{\rm diag}(\lambda _{\rm B}^m,\ldots,\lambda _{\rm B}^1)]$ $[S_{\mathrm{L}},S_{\mathrm{B}}]=[\mathrm{diag}({\lambda }_{\mathrm{L}}^n,...,{\lambda }_{\mathrm{L}}^1),\mathrm{diag}({\lambda }_{\mathrm{B}}^m,...,{\lambda }_{\mathrm{B}}^1)]$ , where [TeX:] $\lambda _{\rm L}^i$ ${\lambda }_{\mathrm{L}}^i$ is the MSV of the i’th layer of the lesion region, and [TeX:] $\lambda _{\rm B}^j$ ${\lambda }_{\mathrm{B}}^j$ is the MSV of the j’th layer of the background. Note that the MSVs are rearranged inversely in the order from bottom to top layer in both regions to compensate the decrease of the photon density with depth at the next step. Then, the forward model is rewritten as follows:

Consequently, by applying the new weight matrix multiplied by the scaling matrix, W ^″, in the iterative inversion, depth dependency of the reconstructed layers is relatively compensated. Finally, the result of the inversion is scaled with the same scaling matrix to obtain the original total absorption map of [M _L, M _B] = [S _L, S _B][M _L, M _B]′′.

2.2.

Phantom and Clinical Experiments

The experiments were performed with our frequency domain near infrared (NIR) system that consisted of source and detection subsystems as well as a handheld probe.²⁶ The source has four laser diodes of wavelengths 740, 780, 808, and 830 nm. In this study, the measurements obtained from 780-nm laser diode was used because an optical absorption map was reconstructed to demonstrate the depth correction approach. Each – laser diode was modulated at 140 MHz and sequentially coupled into nine optical source fibers through 4×1 and 1×9 optical switches. The reflected light from tissue associated with each of the illumination source was detected with 14 optical fiber bundles of 3 mm diam, simultaneously. The detection fibers coupled the light into 14 parallel photomultiplier tubes and electronics channels. The configuration of the source and detector fibers distributed on the handheld probe with a coregistered ultrasound transducer in the middle, is shown in Fig. 1.

Fig. 1

Configuration of source and detector optical fibers distributed on the handheld probe. The source and detector locations are shown with stars and circles, respectively. A commercial ultrasound transducer is located in the middle open slot.

In phantom experiments, 0.8% intralipid solution emulated the background breast tissue optical properties with the calibrated absorption coefficient μ_a in the range of 0.02–0.03 cm⁻¹ and reduced scattering coefficient [TeX:] $\mu ^{\prime} _{\rm s}$ ${\mu }_{\mathrm{s}}^{\prime }$ in the range of 6.9–7.8 cm⁻¹. The solid plastisol phantoms with the calibrated optical properties of μ_a = 0.158 cm⁻¹, [TeX:] $\mu ^{\prime} _{\rm s} = 6.53\ {\rm cm}^{ - {\rm 1}}$ ${\mu }_{\mathrm{s}}^{\prime }=6.53{\mathrm{cm}}^{-1}$ , and μ_a = 0.075 cm⁻¹, [TeX:] $\mu ^{\prime} _{\rm s} = 5.6\ {\rm cm}^{ - {\rm 1}}$ ${\mu }_{\mathrm{s}}^{\prime }=5.6{\mathrm{cm}}^{-1}$ were used to model the high- and low-contrast lesions, respectively. The optical reflectance data of the phantom submerged in the intralipid at different depths were acquired with the NIR system, while the coregistered ultrasound provided the location information. The surface of the handheld probe was considered as the reference plane for the depth measurement. Figure 2 displays different types of phantoms used in the reported experiments, including homogenous high-contrast (HH), homogenous low-contrast (HL), top high-contrast (TH), and top low-contrast (TL) targets. Each of the homogenous and non-homogenous phantoms shown in Fig. 2 had two sizes of diameter, 2 and 3 cm. The imaging volume considered in the reconstruction was 8 × 8 × 4 cm³. The voxel sizes of 0.25 × 0.25 × 0.25 cm³, and 1.5 × 1.5 × 1.5 cm³ were applied to the lesion and background regions in x-y-z plane, respectively.

Fig. 2

Experimental phantoms of (a) homogenous high contrast (HH), (b) homogenous low contrast (HL), (c) top high contrast (TH), and (d) top low contrast (TL).

For the clinical studies, the protocol is approved by the local Institution Review Board and all patients signed the informed consent. Patients were scanned in a supine position while multiple sets of optical reflectance measurements were made with coregistered ultrasound images at the lesion location and the normal contralateral location of the same quadrant as the lesion. As mentioned in Sec. 2, the optical properties of the background required for computing the weight matrix were estimated from measurements obtained at the normal contralateral location.

2.3.

Evaluation of Reconstructed Images

The performance of the depth-correction method is evaluated using phantom targets and clinical cases. To quantitatively compare the results without and with the correction method, the reconstructed absorption maps for each type of lesion shown in Fig. 2 are given in Sec. 3. The spatial map in the x-y plane is shown with 0.5- and 0.25-cm increments in depth for the phantom lesions and clinical cases, respectively. To further quantify the improvement with the correction method, maximum absorptions of the layers in the lesion region are extracted from the corresponding reconstructed absorption maps. The standard box and whisker plots of this data are presented, in which the first and third quartiles of data are at the ends of the box, the median is indicated with a vertical line in the interior of the box, and the maximum and minimum values are at the ends of the whiskers.²⁸

The location of the box within the whisker provides insight on the data distribution that implies the depth sensitivity in our case. Normally, the maximum reconstructed absorption is comparable to the true value. Therefore, improving depth sensitivity mainly raises the recovered values of deeper layers of a homogenous lesion. This effect is indicated in the plots by the increase of the median or shift of the box toward the maximum level. Furthermore, the box size is related to the data variation such that the smaller box size is due to the smaller absorption variation inside the lesion region. The other factor calculated for homogenous cases is the position error (PE), defined as the absolute difference between the depth of the maximum reconstructed absorption and the true depth of the lesion center.

To quantify the accuracy of characterization of the nonhomogenous phantoms, the ratio of the absorption is calculated in lesion region as defined by

The expected values of this ratio are R = 1 for the homogenous phantoms, R = 2 for TH contrast, and R = 0.5 for TL contrast phantoms.

3.1.

Homogenous Phantom Experiments

The absorption map of the high- and low-contrast homogenous phantoms of 2 cm diam, reconstructed without and with using the depth correction is shown in Fig. 3. The slices in each panel showing the spatial map in x-y plane at depths of 0.5–3 cm with 0.5-cm increments. The phantom center is located at depth of 2 cm, measured from the probe surface. Without depth correction, the top layer at depth of 1.25 cm has the maximum absorption of 0.1345 for the HH phantom and 0.104 for the HL one, which are about 84 and 72% of the expected values. With depth correction, the center layer at depth of 2 cm has the maximum absorption of 0.1356 for HH and 0.102 for HL phantom, which are about 85 and 73% of the expected values. Although the reconstructed maximum absorption values are similar before and after correction, the target mass is reconstructed at the center depth rather than at the top portion of the target.

Fig. 3

Reconstructed absorption map of a homogenous phantom of 2 cm diam with the center located at depth of 2 cm. (a, b) High-contrast phantom without and with depth correction, respectively; (c, d) low-contrast phantom without and with depth correction, respectively. Each image consists of six subimages. Each subimage is the spatial x-y image of 8×8 cm reconstructed at the depth marked on its title. The increment in depth is 0.5 cm. For the absorption maps reported in this paper, the spatial dimensions of the subimage were kept the same and the reconstruction depth might be different based on the target location.

The same experiment was performed for the HH and HL phantoms of 2 cm diam located at depths of 1.5, 2.5, and 3 cm, in addition to the HH and HL phantoms of 3 cm diam located at depths of 2, 2.5, and 3 cm. For all cases, the maximum absorptions of target layers were extracted from the reconstructed maps. Then, the box and whisker plots for that data are presented in Figs. 4, 5. These box plots show that the median of the reconstructed absorption inside the phantoms increased toward the maximum absorption after applying depth correction. This improvement is calculated for all phantoms, and the result is listed in Table 1. It is seen that the median is closer to the maximum for the HL phantoms in comparison to HH phantoms of the same size. It is also demonstrated in the box plots that the box size reduced due to applying depth correction. This confirms the decrease of the absorption variation inside the homogenous phantoms. The box-size reduction on average is about two times for the 2-cm-diam targets and 1.8 times for the 3-cm-diam ones. In general, the box size increases with increasing depth, while it is more consistent by applying the depth correction, indicating more uniform depth sensitivity. Overall, the boxes of the bigger phantoms are larger that demonstrate the increase of absorption variation inside the larger phantoms. It is also represented in Figs. 5b and 5d that, unlike the HH cases shown in Fig. 4, the depth sensitivity is increased for the bottom layers of HL phantoms located at depth of 3 cm.

Fig. 4

The box and whisker plot of the maximum absorption of the layers in the lesion region. Homogenous high contrast of size 2 cm diam reconstructed (a) without and (b) with depth correction. Homogenous high-contrast lesion of size 3 cm diam reconstructed (c) without and (d) with depth correction.

Fig. 5

The box and whisker plot of the maximum absorption of the layers in the lesion region. Homogenous low-contrast lesion of size 2 cm diam reconstructed (a) without and (b) with depth correction. Homogenous low-contrast lesion of size 3 cm diam reconstructed (c) without and (d) with depth correction.

Fig. 10

The box and whisker plot of the maximum absorption of the layers in the lesion region. Malignant cases reconstructed (a) without and (b) with depth correction. Benign cases reconstructed (c) without and (d) with depth correction.

Table 1

Compare the median of the data in the lesion layers to its maximum absorption reconstructed without and with depth correction.

Lastly, the PE is evaluated for all cases. The error is in the range of the phantom radius, i.e., 0.75 cm for 2-cm-diam phantoms and 1.25 cm for the 3-cm-diam ones that decreases to the range of the reconstruction mesh sizes, i.e., 0.25 cm after using depth correction.

3.2.

Nonhomogenous Phantom Experiments

The TH and TL phantoms of 2 cm diameter, located at depth of 2 cm are also imaged. The maximum reconstructed absorption of TH phantom is 0.152 cm⁻¹ at 1.25-cm depth while it has moved to the depth of 1.75 cm with the value of 0.13 cm⁻¹ by using the depth-correction method. For the TL phantom, there is no information from the bottom half of the target when the reconstruction is performed without depth correction. The maximum absorption of 0.059 cm⁻¹ is measured at the top layer at a depth of 1.25 cm, whereas with depth correction, the maximum absorption is 0.073 cm⁻¹ at a depth of 2.25 cm. Figures 6a and 6b show the depth profiles from 0.5 to 3.5 cm, of four types of phantoms. All of them are 2 cm diam in size and are located at a depth of 2 cm. After applying depth correction, the maximum reconstructed absorption is at the center mass for the homogenous cases, whereas it is located at top of the center line for the TH phantom and below the center line for the TL one. Hence, it is clearly indicated in Fig. 6b that the phantoms can be more accurately characterized. In order to quantify this effect, the ratio R is calculated for the phantoms of 2 cm in size located at different depths and the results are shown in Table 2. It is seen that R is always >1 for the reconstructions without depth correction because the maximum absorption is always at the top layer. Whereas, by applying depth correction, R gets closer to the expected values for the targets located at depths up to 2.5 cm. At a depth of 3 cm, the bottom half value is close to the background value and R is still >1 for the TL case. Therefore, with the threshold of R = 1, different types of the phantoms are almost precisely characterized. The same factor, R, is calculated for the phantoms of 3 cm in size, and the results are listed in Table 3. Comparing the result without and with depth correction confirms the improvement in the ratio by applying the correction method for the targets located at a depth of 2 cm. At a depth of 2.5 cm, the ratio calculated for HH target is >1 and is very close to 1 for the TL phantom. Therefore, the threshold of R = 1 was not successful for these cases.

Fig. 6

Depth profiles of reconstructed absorption map (a) without and (b) with depth correction for top low-, homogenous low-, homogenous high-, and top high-contrast, 2-cm-diam phantoms located at a depth of 2 cm from left to right.

Table 2

Ratio of maximum absorption of the top half to the maximum absorption of the bottom half, R, calculated for phantoms of 2 cm diam, reconstructed without (w/o) and with (w) using depth correction.

Table 3

The ratio of maximum absorption of top half to the maximum absorption of bottom half, R, calculated for phantoms of 3cm diam that were reconstructed without (w/o) and with (w) using depth correction.

Before applying the correction to the clinical cases, it is necessary to study its effect on the known inhomogeneous phantom. This experiment was performed for the plastisol phantom shown in Fig. 7a. The middle core had a low contrast half-sphere of 2.5 cm diam embedded in the high-contrast shell phantom of 5 cm diam and height of 2 cm. The absorption coefficient was μ_a = 0.25 cm⁻¹ for the outer shell and μ_a = 0.06 cm⁻¹ for the core. The distance between the bottom of the phantom and the probe surface was 2.3 cm. The result of reconstruction, shown in Figs. 7b and 7c, verified that the depth-correction method kept the inhomogeneous target profile while enforcing better sensitivity for deeper layers. It also shifted the reconstructed maximum absorption of 0.139 cm⁻¹ at a depth of 0.75 cm to the location of central mass at 1.25-cm depth with a maximum absorption coefficient of 0.11 cm⁻¹.

Fig. 7

(a) Schematic of the inhomogeneous phantom, (b, c) Absorption maps reconstructed without and with depth correction, respectively. Each image consists of eight subimages with the depth marked on the subimage title. The increment in reconstruction depth is 0.25 cm in this case.

3.3.

Clinical Studies

The absorption maps of benign and malignant large lesions reconstructed without and with depth correction are evaluated in this section. Figure 8a shows a coregistered ultrasound B-scan of a larger suspicious mass located at 10 o'clock position of the right breast of a 43-year-old woman. The large lesion region is marked with the white dashed line. The borders of the lesion were not clear, and the size estimated by ultrasound was about 3×3 cm in spatial dimensions and 2.0 cm in depth. The fitted background tissue optical properties were μ_a = 0.23 cm⁻¹ and [TeX:] $\mu ^{\prime} _{\rm s} = 2.12\ {\rm cm}^{ - 1}$ ${\mu }_{\mathrm{s}}^{\prime }=2.12{\mathrm{cm}}^{-1}$ . The core needle biopsy revealed that the lesion was an invasive ductal carcinoma. Figures 8b and 8d illustrate the reconstructed absorption maps at depth of 0.5–3 cm with 0.25-cm increments, and Figures 8c and 8e show the reconstructed absorption profiles in depth. Similarly, Fig. 9 shows a coregistered ultrasound and the absorption maps of a large suspicious lesion located at 8 o'clock position of the right breast of a 39-year-old woman. The borders of the lesion were not clear and the estimated dimensions by ultrasound were 5 × 5 cm in spatial dimensions and 2.2 cm in depth. The calibrated optical background was μ_a = 0.03 cm⁻¹ and [TeX:] $\mu ^{\prime} _{\rm s} = 5.93\;{\rm cm}^{ - 1}$ ${\mu }_{\mathrm{s}}^{\prime }=5.93{\mathrm{cm}}^{-1}$ . The core needle biopsy revealed that the lesion was a benign fibrocystic change with a moderate degree of intraductal hyperplasia. For the malignant case, the maximum absorption was 0.44 cm⁻¹ at a depth of 1.25 cm that was shifted to the depth of 1.75 cm with the value of 0.37 cm⁻¹. For the benign case, the maximum was 0.11 cm⁻¹ at a depth 1.25 cm that moved to 2.25 cm with the value of 0.118 cm⁻¹ after using depth correction.

Fig. 9

A benign lesion obtained from a 39-year old woman. (a) A coregistered ultrasound B-scan. (b, c) Reconstructed absorption maps in (b) spatial and (c) depth without depth correction. (d, e) Absorption maps with depth correction.

The box and whisker plot of the maximum absorption of the lesion layers for 10 larger lesions, including five benign and five malignant cases, is shown in Fig. 10. Table 4 provides patient information, lesion size, and center depth as measured by coregistered ultrasound. The box size reduction is ∼2.7 times for malignant cases and two times for the benign ones. The median is 76 and 87% of the maximum value that improves to 96 and 98% after depth correction used for the malignant and benign lesions, respectively. Note that, the benign lesions in general have smaller box size than the malignant ones.

Table 4

Patient information and lesion location size and center depth measured by coregistered ultrasound.