Weighting function effects in a direct regularization method for image-guided near-infrared spectral tomography of breast cancer

Jinchao Feng; Shudong Jiang; Brian W. Pogue; Keith Paulsen

doi:10.1364/BOE.9.003266

1. Introduction

Near-infrared spectral tomography (NIRST) uses near-infrared light (600-1000 nm) to image physiologically relevant optical properties of tissue, for breast cancer diagnosis [1–3] and functional brain mapping [4, 5]. However, NIRST alone suffers from low spatial resolution due to the strongly scattering nature of NIR light and leading to diffuse propagation in tissue [6]. To achieve high spatial resolution, the prior structural information provided by anatomical images such as X-ray/CT [7, 8], ultrasound [9, 10] or MRI [11–14], have been incorporated into NIRST reconstruction algorithms. The most common methods to combine anatomical images into NIRST reconstruction, are hard [15] or soft [16] priors based algorithms.

In hard-priors approach, anatomical images are segmented into several different regions with the different structure features, where each region is assumed to be optically uniform during NIRST reconstruction. With this approach, the number of unknown parameters in NIRST inverse problem is significantly reduced by lumping all nodes within these regions together into just a few homogenous regions [15]. This process has the peripheral benefit of significantly enhancing NIRST accuracy within the localized regions by reducing the ill-posedness of NIRST reconstruction. However, its stability is critically dependent on the accuracy of the structural priors derived from the co-registered images, and performance is degraded when incomplete or distorted structural priors are employed.

An alternative inclusion of prior information is so called “soft-priors” [16], which uses Laplacian-type or Helmholtz-type regularization structure to encode structural images into the inversion matrix regularization and allow to change optical properties within each region. Compared with the hard-priors approach, the soft-priors approach introduces some flexibilities in dealing with the correlation between anatomical prior and optical properties. However, the main shortcoming of the above two techniques is that both methods require manual segmentation to identify regions. This segmentation may lead to the objectivity in the process of combining images. Additionally, the segmentation step can be time-consuming, and requires sufficient segmentation experience to avoid segmentation bias or error.

To overcome the shortcoming of segmenting the MRI images, we have developed a NIRST reconstruction algorithm based on a direct regularization imaging (DRI) method [17, 18]. The simulation results have demonstrated its feasibility and effectiveness [17, 18]. Furthermore, statistical analysis of 24 MRI-guided patient images reconstructed by DRI have already demonstrated DRI’s success in distinguishing malignant from benign lesions [3]. In this method, a uniform weighting function which is also called Heaviside step function was used and works well due to its simple formulation, and, in some cases, has a high performance as compared to other weighting functions (see below).

To further improve DRI, the approach to weight the regularization matrix from the greyscale intensities was examined for the effects of different weighting function upon the reconstructed images, in terms of absolute bias error (ABE), mean square error (MSE), full width of half maximum (FWHM), and the contrast of the tumor to normal surrounding tissue (contrast). In this study, simulations, phantom measurements and clinical patient images were used to assess the outcome.

2. Methods

2.1 Weight functions

Under the assumption that light scattering dominates absorption in breast, light transport in breast can be modeled by the diffusion approximation. Because we utilize continuous-wave (CW) light illumination or dc data, the physical process of NIR light illuminating through a highly scattering medium can be approximated by the steady state diffusion equation (DE), given by [19, 20]

- \nabla \cdot κ (r) \nabla (r) + μ_{a} (r) Φ (r) = q_{0} (r) (r \in Ω)

where

Ω

is the imaged breast tissue,

Φ (r)

is the photon fluence rate at position

r

,

κ = 1 / 3 (μ_{a} + μ_{s}^{'})

is the diffusion coefficient (mm⁻¹),

μ_{a}

is the absorption coefficient (mm⁻¹),

μ_{s}^{'}

is the reduced scattering coefficients (mm⁻¹), and

q_{0} (r)

is the source term.

Here, the boundary condition used for the Eq. (1) is Robin-type condition, which can be expressed as [19]:

Φ (r) + \frac{D}{α} \hat{n} \cdot \nabla Φ (r) = 0 (r \in \partial Ω)

where

\partial Ω

is the boundary of

Ω

,

α

is a term which incorporates reflection as a result of refractive index mismatch at the boundary, and

\hat{n}

is the outer normal on

\partial Ω

.

The goal of the NIRST algorithm is to recover optical properties in biological tissue using measurements of light fluence from the tissue surface. This is a typically inverse problem. And this inversion can be achieved using a Tikhonov minimization. If the measured fluence at the tissue surface is represented by $Φ_{m}$ and the calculated data $f (x)$ for a given µ_a and $κ$ distributions for all source/detector combinations (NM), then the objection function is given by:

χ^{2} (x) = {‖ Φ_{m} - f (x) ‖}_{2}^{2} + λ {‖ L x ‖}_{2}^{2}

where

x

are the optical properties

(μ_{a}, μ_{_{s}}^{'})

,

λ

is the regularization parameter, and

L

is the regularization weight matrix generated from MRI images, acting on the solution

x

.

f (x)

can be solved based on the finite element method with our open source software Nirfast [21].

Minimizing $χ^{2} (x)$ with respect to $x$ , which is achieved by setting $\frac{\partial χ^{2} (x)}{\partial x} = 0$

\frac{\partial χ^{2} (x)}{\partial x} = J^{T} δ - λ L^{T} L x = 0

where

δ = Φ_{m} - f (x)

, termed as the data-model misfit, and

J

the Jacobian matrix. Rewriting Eq. (4) for the kth iteration leads to

J_{k}^{T} δ_{k} - λ_{k} L^{T} L x_{k} = 0

Taking account into $δ_{k} = δ_{k - 1} - J_{k} Δ x_{k - 1}$ [22], and substituting $δ_{k}$ into Eq. (5) results in

[J_{k}^{T} J_{k} + λ_{k} L^{T} L] Δ x_{k} = J^{T} δ_{k}

where

Δ x_{k}

is the update for the optical properties (µ_a, µ_s^’) in the kth iteration. Based on Eq. (6), the update equation for

Δ x_{k}

can be expressed as

Δ x_{k} = {[J_{k}^{T} J_{k} + λ_{k} L^{T} L]}^{- 1} J_{k}^{T} δ_{k}

Since CW measurements are used to reconstruct optical properties in our experiments, we assumed that µ_s^’ is known and uniform, and only the absorption coefficient µ_a is recovered. For the regularization parameter $λ$ at the kth iteration, it was setting as $λ_{k} = 10 * \max (d i a g (J_{k}^{T} J_{k}))$ .

The weight matrix, L, has the form of:

L_{i j} = {\begin{cases} 1 i = j \\ - \frac{1}{M_{i}} \exp (- \frac{{| γ_{i} - γ_{j} |}^{2}}{2 σ_{g}^{2}}) g (d_{i j}) o t h e r w i s e \end{cases}

where

γ_{i}

is the grayscale value in the MRI images mapped to the node i in the finite element mesh,

σ_{g}

is the characteristic grayscale difference over which to apply regularization, and

M_{i}

is a factor chosen for the ith row, and satisfies

\sum L_{i j} = 0

. The function,

g (d_{i j})

, is a weighting function about the distance

d_{i j}

between the nodes i and j, which determines the local weight applied to the ith node of the finite element mesh. As for the criterion of choosing the weighting functions, the value of the weighting function should be non-negative and inversely weighted by the distance

d_{i j}

, i.e., the maximum value of the weighting function should be at or approach the zero distance, and the function should decay as the distance increases for eliminating discontinuities in reconstructed NIRST images. Table 1 shows the nine different functions of

g (d_{i j})

used in Eq. (4) for comparison their effects on reconstructed NIRST images.

Table 1. The nine different functions used in the study.

View Table | View all tables in this article

A simple weighting function is the uniform function (Function 1), which has been adopted into the previous study [18]. In this case, the function value of $g (d_{i j})$ is constant when the distance $d_{i j}$ is smaller than a threshold. The weights were ignored in this function.

Function 2 is the exponential function, and it has infinite extent. A related Function 3 is a Gaussian function. This function also has infinite extent. Functions 4 and 5 are the simple weighting function that just raise the distance to a negative power. The magnitude of the power determines the rate of drop off of the weight with distance. These two weighting functions go to infinity as the node i approaches the node j. Here, we just consider p = 1 (Function 4) and p = 2 (Function 5).

Function 6 is an inverse distance function and its alternative is the linear weight function (Function 7). The quadratic weight function (Function 8) is also a commonly used function, and Function 9 is a tricube weight function.

For the purposes of objectively evaluating the performances of different weighting functions, the distance $d_{i j}$ was normalized in this study, as shown in Fig. 1.

Fig. 1 Plots of the 9 weighting functions used in Eq. (4) with normalized distance (x-axis) investigated in this study. a. u., arbitrary unit.

Download Full Size | PDF

2.2 Quantitative NIRST image comparison

We conducted a series of numerical experiments to access the performances of these functions. The absolute bias error (ABE), mean square error (MSE) [23, 24], peak-signal-to-noise ratio (PSNR) [25], contrast-to-noise ratio (CNR) [26], FWHM and contrast were used to quantitatively compare the NIRST image reconstructed with 9 weighting functions. These parameters were defined as

A B E = \frac{\sum_{i = 1}^{N} | μ_{i} - μ_{r e c o n} |}{N}

V a r = \frac{{\sum_{i = 1}^{N} | μ_{i} - {\bar{μ}}_{r e c o n} |}^{2}}{N}

M S E = A B E^{2} + V a r

P N S R = 10 \log_{10} (\frac{{(M a x (μ))}^{2}}{M S E})

C N R = \frac{μ_{r o i} - μ_{b a c k}}{{(ω_{r o i} σ_{{_{μ}}_{_{r o i}}}^{2} + ω_{b a c k} σ_{{_{μ}}_{_{b a c k}}}^{2})}^{\frac{1}{2}}}

where

μ_{i}

and

μ_{r e c o n}

are the true and reconstructed absorption coefficients at the finite node i respectively,

{\bar{μ}}_{r e c o n}

is the average value of reconstructed absorption coefficients,

N

is the total number of finite element nodes,

μ_{r o i}

and

μ_{b a c k}

represent the mean of the reconstructed absorption coefficient, in the region-of-interest (ROI) and the background, respectively. The

σ_{μ}_{_{r o i}}

and

σ_{μ}_{_{b a c k}}

represent the standard deviations of the reconstructed absorption coefficient in the ROI and background, respectively. Here the weights

ω_{r o i} = \frac{A_{r o i}}{A_{t o t a l}}

and

ω_{b a c k} = \frac{A_{b a c k}}{A_{t o t a l}}

represent the ratio of areas between the background and total area as well as ROI and total area, respectively. The PSNR is used to compare the restoration of the images, without depending strongly on the image intensity scaling. The CNR indicated whether the inclusion could be clearly seen in the reconstructed image. We expect lower ABE, and MSE, and higher PSNR, CNR and contrast, when there is better image quality.

Note that the expected/target absorption coefficient values were not feasible to obtain in some phantoms and in all patient cases, and therefore the CNR and contrast were used as metrics of success in these cases.

2.3 Numerical simulation

The 2D circular simulation phantoms that have a diameter of 80 mm are used for ease of comparison. The absorption coefficient (µ_a) and the reduced scattering coefficient (µ_s’) of the phantom were 0.01mm⁻¹ and 1.0 mm⁻¹, respectively. A total of 16 sources and 16 detectors were evenly arranged along the circumference of the phantom. For each source illumination, data was collected at 16 detector locations which lead to a total of 256 measurements.

Figure 2 shows the geometries and MRI images of all simulation phantoms. In the first simulation (study 1), a single inclusion with a diameter of 15 mm was added into the phantom to show the effect of the weighting function on the accuracy of the reconstructed absorption images, as shown in Fig. 2(a). In this and following numerical simulations, the optical properties of the inclusion were set to be µ_a = 0.02 mm⁻¹ and µ_s^’ = 1.0 mm⁻¹. As shown in the bottom row of Fig. 2(a), the intensities of the inclusion and the background were set to be 80 and 50, respectively to simulate the type of DCE-MRI contrast commonly observed [17].

Fig. 2 The geometries of the phantom (top row) and the corresponding MRI images (bottom row) used for the simulation studies 1 and 2. The absorption and scattering coefficients of a homogeneous phantom in study 1(a), and simulation study 2(b) were 0.01 mm⁻¹ and 1.0 mm⁻¹ in background and inclusions were 0.02 mm⁻¹ and 1.0 mm⁻¹, respectively.

Download Full Size | PDF

In the second simulation (study 2), two inclusions with the same diameter of 15 mm were embedded at (50, 50) and (70, 50) with an edge-to-edge distance of 10 mm as shown in Fig. 2(b). The corresponding MRI image was also shown in the bottom row of Fig. 2(b). The intensities of the MRI image in the inclusion region and the background were the same as that in study 1 and the contrast was still 1.6.

To generate simulated data, the simulation phantoms were discretized each with a finer mesh with 7909 nodes and 15503 triangular elements. 5% random noise was added to the measurement data. To avoid the so called ‘inverse crime’ of using simulation data with no noise or no change in meshes, NIRST reconstruction was performed on a coarser mesh with 2001 nodes and 3867 triangular elements. The pixel basis used in the reconstruction was 30 × 30. The NIRST reconstruction stops if the difference between the forward data and the reconstructed data do not decrease by more than 1% for successive iterations or the maximum number of iteration (50) is reached. The initial regularization parameter is set to be 10.

2.4 Phantom study

To further evaluate the effect of weighting function on reconstructed images, a three-layer gelatin phantom study has been performed. The outer layer and middle layer were gelatin with different concentration of blood and the inner layer was a 25 mm-diameter cylindrical hole with a blood solution. The contrast in total hemoglobin (HbT) concentration from inner layer to outer layer was about 2:1:0.5 and the water concentration in the solution contrast was around 95%. Gadolinium (Omniscan gadodiamide) was used as a MR contrast agent, which was added to the middle layer of the phantom to create contrast in the MR images. The optical data acquisition plane was marked with a MR sensitive fiducial marker, as show in Fig. 3(a). CW measurements were taken from 700 nm to 900 nm at 13 wavelengths (700, 720, 740, 750, 765, 780, 790, 800, 820, 840, 860, 880 and 900 nm) using a spectrometer tomography system with 16 sources and 16 detectors, indicated with ’o’ and ’x’ in Fig. 3(b) [27]. An FEM mesh with 1780 nodes and 3404 triangles was generated from the T1-weighted spin echo MRI image shown in Fig. 3(a) for image reconstruction. And the structural information provided by the MRI image was encoded into NIRST reconstruction.

Fig. 3 (a) The MRI T1 image of the phantom and (b) a schematic of the sources and detectors setup. The red arrows represent fiducial markers, and ‘o’ and ‘x’ denote sources and detectors, respectively.

Download Full Size | PDF

2.5 Patient imaging study

Patient data was collected through an imaging protocol for human subject studies approved by the Institutional Review Board (IRB) for the Committee for the Protection of Human Subjects at Dartmouth College and at Xijing Hospital. Written consent was obtained after the protocol was explained to a subject. NIRST data and MRI images were concurrently acquired by a NIRST/MRI imaging system developed at Dartmouth College. The details of our imaging system can be found in our previous publications [13, 28]. After data acquisition, MRI images were processed by the open source software NIRFASTSlicer [29] to generate a uniform mesh, which is discretized from the T1-weighted MRI volume. Then MRI DCE images was incorporated into NIRST reconstruction to guide optical reconstruction. Since it is impossible to obtain true optical properties, so we use HbT contrast as a metric for comparison. To calculate contrast, we defined the region-of-interest (ROI) manually based on MRI DCE images.

3. Results

3.1 Simulation

Figure 4(a) shows the reconstructed contrasts between the inclusion and the background using different weighting functions with different σ_g in the case of a single inclusion study. The average reconstructed contrast of µ_a (the ratio of average value of recovered µ_a of inclusion to background) for all weighting functions, decreased from 1.93 to 1.23 when σ_g increased from 0.001 to 10. Considering when σ_g = 0.001_, the reconstructed µ_a were overestimated significantly, while µ_a were underestimated at larger σ_g (0.1, 1 and 10), the optimal σ_g = 0.01 was used in the following experiments.

Fig. 4 The reconstructed contrast with varied σ_g for 9 functions (a) and the profiles from the reconstructed images through the center of the inclusion and along the x-axis when σ_g = 0.01, (b) in simulation study 1.

Download Full Size | PDF

Figure 4(b) shows the corresponding cross-section profiles through the center of the inclusion and along the x-axis. The results demonstrate that the weighting function has an effect on the reconstructed absorption coefficients, and the reconstructed µ_a in the inclusion region are overestimated compared to their true values, with any of the 9 functions. The reconstructed average value of µ_a in the inclusion with each of the 9 weighting functions varied from 0.019 to 0.021mm⁻¹, which is within ± 5% of the true value. The maximum average value of reconstructed µ_a is achieved using the Function 5 with 0.021 mm⁻¹ and the minimum average value is using the Function 1 with 0.019 mm⁻¹. The corresponding quantitative comparison results are given in Table 2. As it shown in Table 2, it can be seen that Functions 4 and 5 are better than other functions by about higher PSNR and CNR. For FWHM, and contrast, the differences between any two functions are within ± 3%. However, the values of ABE and MSE of the Function 5 are about 11% & 24% higher than those of the other functions, as it is evident from the profile plot that the variation from the true distribution is far more over estimated when compared with that of the other functions.

Table 2. Quantitative comparisons in the case of single inclusion.

View Table | View all tables in this article

Figure 5 shows the reconstructed images of absorption coefficient with σ_g = 0.01 in the case of two inclusions. And Fig. 6 shows the corresponding cross-section profiles through the center of the inclusions and along the x-axis. The two inclusions are observable and reconstructed with their centers at the correct positions. Despite the comparatively small differences produced by all 9 functions, the accuracy of the reconstructed absorption coefficient can indeed be influenced by using different functions. The peak values of reconstructed µ_a using the 8 weighting functions are all on 0.021 mm⁻¹ of both central and near the boundary inclusions (except for the Function 5 (0.022 & 0.023 mm⁻¹)), which is about 5% higher than the true target values. As it shown in Table 3, the average values of reconstructed µ_a using all 9 weighting functions are all 0.018 mm⁻¹ for both the central and near the boundary inclusions, which is 10% lower than that of the true targeted values. Similar to the previous one inclusion simulation result, the Functions 4 and 5 can have a PSNR gain of about 0.2dB over other seven functions while the higher values of CNR are obtained. However, the Function 5 produced the largest bias error and MSE because the absorption coefficient has again been overestimated with a peak value of 0.023 mm⁻¹, especially for the inclusion near the boundary.

Fig. 5 Reconstructed images for different weighting functions. (a) - (i) are the reconstructed absorption images using different weighting functions (Function 1 to Function 9) in the case of two inclusions.

Download Full Size | PDF

Fig. 6 The profiles of reconstructed absorption coefficient from the reconstructed images in simulation study 2, which are through the centers of the inclusions and along the X-axis.

Download Full Size | PDF

Table 3. Quantitative results in the case of two inclusions.

View Table | View all tables in this article

3.2 Phantom study

Figure 7 shows the reconstructed results of the gelatin phantom using different weighting functions. The reconstructed images shown that the location of inclusion of HbT can be recovered very accurately, and the trend of HbT concentration changes was well recovered for all three layers. The difference between inner layer and middle layer could be differentiated and the HbT contrast between the inner layer and middle layer are summarized in Table 4. In this case, although the Function 5 produces highest HbT contrast, which is very close to true value, artifacts are observed just below the inclusion. Function 4 reduces the artifacts and has offered more than 13% improvement in CNR, compared with those of other forms. Considering the contrast, with the imaging artifact and FWHM, the weighting functions 1 or 4 are the better solution for the reconstruction.

Fig. 7 Resulting HbT images from a gelatin phantom with one inclusion was used for evaluation. The images were reconstructed with different weighting functions, as shown where (a) - (i) are the results using each of the Function 1 to Function 9, respectively.

Download Full Size | PDF

Table 4. The quantitative results of reconstructed HbT using different weighting functions for a gelatin phantom study.

View Table | View all tables in this article

3.3 patient study

Case 1: a 58-year-old woman with an undiagnosed 15x25x42 mm³ lesion in her left breast, which had a pathologically confirmed diagnosis of malignancy. Her pre-biopsy BIRADS score (from combined all MRI sequences) was 5. The breast was imaged with MR-guided NIRST, and Fig. 8(a) displays the 3D rendering result using NIRFASTSlicer from T1 volume image where the fiber locations are evident from the tissue depressions of the breast surface and the fiducial markers. Figure 8(b) shows a representative MR image slice from the standard T1 sequence, and Fig. 8(c) is a DCE-MR image, clearly showing the lesion location is about at the center of the breast. Figure 8(d)-(l) show reconstructed NIRST HbT images overlaid on the corresponding T1 scans guided by MRI DCE images with 9 different functions. The results reveal that the DRI method can localize the tumor accurately in spite of the functional form. And the reconstructed contrasts between the lesion and background using different weighting functions are shown in Table 5. The estimated HbT contrasts in the suspicious region obtained for all functions were higher than the surrounding normal tissue, and suggested that the tumor was malignant according to prior reports [28]. In this case, the weighting Function 5 had the highest lesion contrast of 4.4. The smallest HbT contrast was 1.6 when using either of Functions 1, or 3 or 9, which is consistent with the results obtained from simulation and phantom studies. The improvement in the reconstruction with the Function 4 or 5 is obvious in terms of the CNR.

Fig. 8 The first patient example shown by: (a) 3D volume rendering; (b) T1 MRI; (c) DCE MRI; (d) - (l) are reconstructed HbT images overlaid on the T1 MRI cross-section using the Functions 1 to 9, respectively. The red arrow in (c) indicates the tumor.

Download Full Size | PDF

Table 5. HbT contrasts and CNR values are listed for patient case 1.

View Table | View all tables in this article

Case 2: a 24-year-old woman with an undiagnosed 23x40x70 mm³ lesion in her left breast with later pathologically confirmed malignant, was imaged with MR-guided NIRST. Her BIRADS score (from combined all MRI sequences) was 5. Figure 9(a)-(c) display the 3D rendering, MRI T1 image and MRI DCE images, respectively. Figure 9(d)-(l) show NIRST HbT reconstructed images overlaid on the corresponding T1 scans guided by MRI DCE images based on different weighting functions. The HbT images show that the lesion was well localized for all functions. The HbT contrasts in the suspicious region obtained for all functions were also higher than 1.0, which indicates that the abnormality was malignant according to prior studies [28]. In this case, the Function 4 had the highest lesion contrast and the contrast was 2.4. The smallest HbT contrast was obtained for the Function 5 (1.5) because the maximum value of HbT was obtained near the position of boundary and the HbT for the lesion which was far away from the boundary was suppressed. Similar with previous studies, the Function 4 yielded more than 5% improvement in CNR as compared with those of other weighting functions.

Fig. 9 The second patient example shown by: (a) 3D volume rendering; (b) T1 MRI; (c) DCE MRI; (d) - (l) are reconstructed HbT images overlaid on the T1 MRI cross-section using the Functions 1 to 9, respectively. The red arrows in (c) indicate the tumor location.

Download Full Size | PDF

4. Discussion and conclusions

The motivation for a direct regularization imaging methodology based within NIRST arises from the need to reduce computational complexity by eliminating the need to segment tissue areas. However, the direct regularization imaging method depends on a pre-specified weighting function. To the best of our knowledge, the only existing functional form in published DRI guided NIRST reconstruction algorithm is the uniform function [3, 18]. However, this form may not be the best choice since it does not optimize the weights. To our best knowledge, the effect on quantitative accuracy and recovered contrast of NIRST images based on different weighting functions has not been systematically studied. In order to better understand whether there are other weighting functions that can gain better performance in this NIRST reconstruction, in the manuscript, the criteria of choosing weighting function was presented for the first time. And the 9 weighting functions in this study are created based on our experience. Although we believe that there are numerous other weighting functions which can be used for DRI based NIRST reconstruction, the 9 weighting functions that we studied are more common. Both phantom and clinical results have demonstrated that a significant improvement on reconstructed images will be gained when an appropriate weighting function is adopted. This is shown both quantitatively and qualitatively in the results. Based on the simulations, phantoms and patient results shown in this manuscript, a guideline for choosing a weighting function is provided for DRI based NIRST breast imaging.

For example, if only a single inclusion needs to be recovered in reconstruction, and one wants to obtain the best quantitative accuracy, a uniform function should be encouraged. The Function 4 or 5 are better in terms of PSNR, CNR and contrast, as compared to the other functions. But if the imaging problem contains two or multiple inclusions or when the tumor heterogeneity can be seen from high resolution MRI, the use of the Function 5 may not be optimal (Fig. 9), as it produces higher value close to the boundary (compared with other weighting functions), and suppresses the peak absorption coefficient of other lesions. In particular, with respect to the phantom and patient image recovery results, the Function 4 achieves the best overall performance in terms of CNR and contrast (Figs. 7–9, Tables 4–6).

Table 6. HbT contrasts and CNR values are listed for patient case 2.

View Table | View all tables in this article

The tumor size effect was also studied. The tumor size varied from 5 mm to 25 mm with a step size of 5 mm. When the tumor size was 5 mm, the Function 4 improves the values of ABE, MSE, and CNR more than 5.3%, 8.1%, and 4.1%, respectively, compared with other weighting functions. In this case, the Function 5 had the highest PSNR of 27.5, and the Function 4 had a higher PSNR of 26.9. When the tumor size is between 10 and 20 mm. These results are agreed with the results as that when tumor size of 15 mm (Fig. 4(b)). However, when tumor size is larger than 20 mm, the average values of PSNR decreased 0.7%, and 2.8% for Functions 4 and 5, respectively while the other performances of Functions 4 and 5 are consistent with that when tumor size of <20 mm. The representative results are shown in Fig. 10.

Fig. 10 The plots of (a) ABE, (b) MSE, (c) PSNR and (d) CNR with increased target size in the single inclusion simulation experiment are shown.

Download Full Size | PDF

Note that the distance $d_{i j}$ between the nodes i and j were normalized by the maximum distance, and the distance $d_{i j}$ varied in the range of [0, 1]. However, the distance $d_{i j}$ was be truncated when it became smaller than a threshold value. We tested the effect of truncated distance $d_{i j}$ on the reconstructed images with the phantom experiment, and the corresponding results were shown in Fig. 11. When $d_{i j}$ was truncated with a threshold value of 0.3, there was a 6% improvement in the HbT contrast and 11% improvement in the CNR for the Function 4. And the FWHM was reduced from 25.3 mm to 24.8 mm, for a true value of 25 mm. Similar results were also observed from patient cases. For example, the HbT contrast was improved from 2.7 to 3.0 and CNR was improved from 7.6 to 7.8 with the truncated threshold of 0.3 for the first patient study.

Fig. 11 The plots of CNR and HbT contrast with the truncated threshold value in the phantom experiment are shown, indicating an optimal threshold value at 0.3.

Download Full Size | PDF

Even though the cases considered here were limited, the aim in this work was to show that the function used for this image-guidance weighting can have an effect upon the NIRST image recovery. Note that as the quantitative comparison of phantom and clinical results requires true information about the tissue optical properties, and the HbT contrast was used as the metric of success here. Future work could include extending this study to more patient data cases to show their effect on differentiating the malignant versus benign tumors.

Funding

National Natural Science Foundation of China (81370038); National Cancer Institute (R01 CA176086 and R01 CA069544).

Disclosures

The authors declare that there are no conflicts of interest related to this article.

References and links

1. S. D. Konecky, G. Y. Panasyuk, K. Lee, V. Markel, A. G. Yodh, and J. C. Schotland, “Imaging complex structures with diffuse light,” Opt. Express 16(7), 5048–5060 (2008). [CrossRef] [PubMed]

2. E. A. Lim, J. E. Gunther, H. K. Kim, M. Flexman, H. Hibshoosh, K. Crew, B. Taback, J. Campbell, K. Kalinsky, A. Hielscher, and D. L. Hershman, “Diffuse optical tomography changes correlate with residual cancer burden after neoadjuvant chemotherapy in breast cancer patients,” Breast Cancer Res. Treat. 162(3), 533–540 (2017). [CrossRef] [PubMed]

3. J. Feng, J. Xu, S. Jiang, H. Yin, Y. Zhao, J. Gui, K. Wang, X. Lv, F. Ren, B. W. Pogue, and K. D. Paulsen, “Addition of T2-guided optical tomography improves noncontrast breast magnetic resonance imaging diagnosis,” Breast Cancer Res. 19(1), 117 (2017). [CrossRef] [PubMed]

4. M. A. Franceschini and D. A. Boas, “Noninvasive measurement of neuronal activity with near-infrared optical imaging,” Neuroimage 21(1), 372–386 (2004). [CrossRef] [PubMed]

5. C. W. Lee, R. J. Cooper, and T. Austin, “Diffuse optical tomography to investigate the newborn brain,” Pediatr. Res. 82(3), 376–386 (2017). [CrossRef] [PubMed]

6. S. R. Arridge, “Optical tomography in medical imaging,” Inverse Probl. 15(2), R41–R93 (1999). [CrossRef]

7. 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(1), 89–97 (2011). [CrossRef] [PubMed]

8. R. Baikejiang, W. Zhang, and C. Li, “Diffuse optical tomography for breast cancer imaging guided by computed tomography: A feasibility study,” J. X Ray Sci. Technol. 25(3), 341–355 (2017). [CrossRef] [PubMed]

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

10. 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(1), 5–18 (2016). [CrossRef] [PubMed]

11. V. Ntziachristos, A. G. Yodh, M. Schnall, and B. Chance, “Concurrent MRI and diffuse optical tomography of breast after indocyanine green enhancement,” Proc. Natl. Acad. Sci. U.S.A. 97(6), 2767–2772 (2000). [CrossRef] [PubMed]

12. 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. U.S.A. 103(23), 8828–8833 (2006). [CrossRef] [PubMed]

13. M. A. Mastanduno, J. Xu, F. El-Ghussein, S. Jiang, H. Yin, Y. Zhao, K. E. Michaelsen, K. Wang, F. Ren, B. W. Pogue, and K. D. Paulsen, “Sensitivity of MRI-guided near-infrared spectroscopy clinical breast exam data and its impact on diagnostic performance,” Biomed. Opt. Express 5(9), 3103–3115 (2014). [CrossRef] [PubMed]

14. 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(23), 6006–6015 (2014). [CrossRef] [PubMed]

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

16. P. K. Yalavarthy, B. W. Pogue, H. Dehghani, C. M. Carpenter, S. Jiang, and K. D. Paulsen, “Structural information within regularization matrices improves near infrared diffuse optical tomography,” Opt. Express 15(13), 8043–8058 (2007). [CrossRef] [PubMed]

17. L. Zhang, Y. Zhao, S. Jiang, B. W. Pogue, and K. D. Paulsen, “Direct regularization from co-registered anatomical images for MRI-guided near-infrared spectral tomographic image reconstruction,” Biomed. Opt. Express 6(9), 3618–3630 (2015). [CrossRef] [PubMed]

18. J. Feng, S. Jiang, J. Xu, Y. Zhao, B. W. Pogue, and K. D. Paulsen, “Multiobjective guided priors improve the accuracy of near-infrared spectral tomography for breast imaging,” J. Biomed. Opt. 21(9), 090506 (2016). [CrossRef] [PubMed]

19. S. R. Arridge and J. C. Hebden, “Optical imaging in medicine: II. modelling and reconstruction,” Phys. Med. Biol. 42(5), 841–853 (1997). [CrossRef] [PubMed]

20. A. P. Gibson, J. C. Hebden, and S. R. Arridge, “Recent advances in diffuse optical imaging,” Phys. Med. Biol. 50(4), R1–R43 (2005). [CrossRef] [PubMed]

21. H. Dehghani, M. E. Eames, P. K. Yalavarthy, S. C. Davis, S. Srinivasan, C. M. Carpenter, B. W. Pogue, and K. D. Paulsen, “Near infrared optical tomography using NIRFAST: algorithm for numerical model and image reconstruction,” Commun. Numer. Methods Eng. 25(6), 711–732 (2009). [CrossRef] [PubMed]

22. P. K. Yalavarthy, B. W. Pogue, H. Dehghani, and K. D. Paulsen, “Weight-matrix structured regularization provides optimal generalized least-squares estimate in diffuse optical tomography,” Med. Phys. 34(6), 2085–2098 (2007). [CrossRef] [PubMed]

23. B. W. Pogue, X. Song, T. D. Tosteson, T. O. McBride, S. Jiang, and K. D. Paulsen, “Statistical analysis of nonlinearly reconstructed near-infrared tomographic images: Part I-Theory and simulations,” IEEE Trans. Med. Imaging 21(7), 755–763 (2002). [CrossRef] [PubMed]

24. X. Song, B. W. Pogue, T. D. Tosteson, T. O. McBride, S. Jiang, and K. D. Paulsen, “Statistical analysis of nonlinearly reconstructed near-infrared tomographic images: part II-experimental interpretation,” IEEE Trans. Med. Imaging 21(7), 764–772 (2002). [CrossRef] [PubMed]

25. A. P. Cuadros and G. R. Arce, “Coded aperture optimization in compressive X-ray tomography: a gradient descent approach,” Opt. Express 25(20), 23833–23849 (2017). [CrossRef] [PubMed]

26. S. C. Davis, B. W. Pogue, H. Dehghani, and K. D. Paulsen, “Contrast-detail analysis characterizing diffuse optical fluorescence tomography image reconstruction,” J. Biomed. Opt. 10(5), 050501 (2005). [CrossRef] [PubMed]

27. J. Wang, S. Jiang, K. D. Paulsen, and B. W. Pogue, “Broadband frequency-domain near-infrared spectral tomography using a mode-locked Ti:sapphire laser,” Appl. Opt. 48(10), D198–D207 (2009). [CrossRef] [PubMed]

28. M. A. Mastanduno, J. Xu, F. El-Ghussein, S. Jiang, H. Yin, Y. Zhao, K. Wang, F. Ren, J. Gui, B. W. Pogue, and K. D. Paulsen, “MR-guided near-infrared spectral tomography increases diagnostic performance of breast MRI,” Clin. Cancer Res. 21(17), 3906–3912 (2015). [CrossRef] [PubMed]

29. http://www.dartmouth.edu/~nir/nirfast/.

Function	Formulation
1	$g (d_{i j}) = 1$
2	$g (d_{i j}) = \exp (- d_{i j})$
3	$g (d_{i j}) = \exp (- d_{i j}^{2})$
4	$g (d_{i j}) = \frac{1}{d_{i j}}$
5	$g (d_{i j}) = \frac{1}{d_{i j}^{2}}$
6	$g (d_{i j}) = \frac{1}{1 + d_{i j}}$
7	$g (d_{i j}) = 1 - d_{i j}$
8	$g (d_{i j}) = {(1 - d_{i j}^{2})}^{2}$
9	$g (d_{i j}) = {(1 - d_{i j}^{3})}^{3}$

	Function 1	Function 2	Function 3	Function 4	Function 5	Function 6	Function 7	Function8	Function 9
ABE	3.8e-4	3.8e-4	3.8e-4	3.8e-4	4.2e-4	3.8e-4	3.8e-4	3.8e-4	3.8e-4
MSE	9.6e-7	9.7e-7	9.7e-7	1.0e-6	1.2e-6	9.7e-7	9.9e-7	9.9e-7	9.9e-7
PSNR (dB)	28.4	28.5	28.5	28.6	28.6	28.5	28.5	28.5	28.5
CNR	10.2	10.2	10.2	10.4	10.8	10.2	10.2	10.2	10.2
FWHM	13.8	13.8	13.8	13.6	13.4	13.8	13.8	13.8	13.8
Contrast	1.9	1.9	1.9	1.9	2.0	1.9	1.9	1.9	1.9

	Function 1	Function 2	Function 3	Function 4	Function 5	Function 6	Function 7	Function 8	Function 9
ABE	1.7e-3	1.7e-3	1.7e-3	1.7e-3	1.8e-3	1.7e-3	1.7e-3	1.7e-3	1.7e-3
MSE	6.9e-6	6.9e-6	6.9e-6	7.0e-6	7.3e-6	6.9e-6	7.0e-6	7.0e-6	7.0e-6
PSNR (dB)	22.6	22.6	22.6	22.8	23.6	22.6	22.7	22.7	22.7
CNR	5.6	5.6	5.6	5.7	6.1	5.6	5.6	5.6	5.6
FWHM 1	15.8	15.8	15.8	15.8	15.4	15.8	15.8	15.8	15.8
FWHM 2	16.2	16.2	16.2	16.0	15.6	16.2	16.2	16.2	16.2
Inclusion1 (mm⁻¹)	0.018	0.018	0.018	0.018	0.018	0.018	0.018	0.018	0.018
Inclusion2 (mm⁻¹)	0.018	0.018	0.018	0.018	0.019	0.018	0.018	0.018	0.018

Function	Formulation
1	$g (d_{i j}) = 1$
2	$g (d_{i j}) = \exp (- d_{i j})$
3	$g (d_{i j}) = \exp (- d_{i j}^{2})$
4	$g (d_{i j}) = \frac{1}{d_{i j}}$
5	$g (d_{i j}) = \frac{1}{d_{i j}^{2}}$
6	$g (d_{i j}) = \frac{1}{1 + d_{i j}}$
7	$g (d_{i j}) = 1 - d_{i j}$
8	$g (d_{i j}) = {(1 - d_{i j}^{2})}^{2}$
9	$g (d_{i j}) = {(1 - d_{i j}^{3})}^{3}$

	Function 1	Function 2	Function 3	Function 4	Function 5	Function 6	Function 7	Function8	Function 9
ABE	3.8e-4	3.8e-4	3.8e-4	3.8e-4	4.2e-4	3.8e-4	3.8e-4	3.8e-4	3.8e-4
MSE	9.6e-7	9.7e-7	9.7e-7	1.0e-6	1.2e-6	9.7e-7	9.9e-7	9.9e-7	9.9e-7
PSNR (dB)	28.4	28.5	28.5	28.6	28.6	28.5	28.5	28.5	28.5
CNR	10.2	10.2	10.2	10.4	10.8	10.2	10.2	10.2	10.2
FWHM	13.8	13.8	13.8	13.6	13.4	13.8	13.8	13.8	13.8
Contrast	1.9	1.9	1.9	1.9	2.0	1.9	1.9	1.9	1.9

Weighting function effects in a direct regularization method for image-guided near-infrared spectral tomography of breast cancer

Abstract

1. Introduction

2. Methods

2.1 Weight functions

2.2 Quantitative NIRST image comparison

2.3 Numerical simulation

2.4 Phantom study

2.5 Patient imaging study

3. Results

3.1 Simulation

3.2 Phantom study

3.3 patient study

4. Discussion and conclusions

Funding

Disclosures

References and links

Cited By

Figures (11)

Tables (6)

Equations (13)

Biomedical Optics Express

	Function 1	Function 2	Function 3	Function 4	Function 5	Function 6	Function 7	Function 8	Function 9
HbT Contrast	1.7	1.6	1.6	1.8	1.9	1.7	1.5	1.5	1.4
FWHM	32.6	33.8	34.3	25.3	35.1	32.8	36.6	37.1	38.6
CNR	3.2	3.0	3.0	3.6	3.0	3.1	2.7	2.6	2.5

	Function1	Function 2	Function 3	Function 4	Function 5	Function 6	Function 7	Function 8	Function 9
HbT contrast	1.6	1.7	1.6	2.7	4.4	1.7	1.7	1.7	1.6
CNR	6.0	6.4	6.2	7.6	8.1	6.3	6.5	6.4	6.2