2.1.

Mie Scattering Theory Interpretation

Mie theory²⁶ provides an exact solution for elastic scatter, in the case of a perfect dielectric sphere of arbitrary size in a uniform background medium. Using Mie theory, the reduced scattering coefficient spectra of bulk homogeneous samples can also be expressed quantitatively, with the homogeneous medium equation as

The average scatterer size could be expressed as

A more empirical approach to Mie theory was first proposed by van Staveren ²⁵ in fitting the scattering spectrum of Intralipid. A number of groups have adopted the approach to characterize the spectrum of the reduced scattering coefficient observed in tissues. ^{37, 40, 41, 42} In these studies, the scattering spectra are considered to satisfy a power law relationship, similar to Rayleigh scattering where the scatter coefficient decreases as ${\lambda }^{-4}$ , but using a power value smaller than 4. Empirically, when there is a broad range of scattering particle sizes, this spectrum is described by a power law curve of the type

2.2.

Scattering from Tissue Spectra

In Eq. 2, the number density only influences the amplitude of the scattering spectra in the NIR range, and there are only two parameters to determine the shape of the scattering spectra including (i) the distribution function and (ii) the relative refractive index. The distribution function could be separated into two parts: the type of the distribution function and the average size. So there are three parameters to determine the shape of the scattering spectra including (i) the type of the distribution function, (ii) the average scatterer size, $⟨a⟩$ , and (iii) the relative refractive index. If the type of the distribution function and relative refractive index can be set by some known information, then there are only two parameters left to estimate in Eq. 2: average scatterer size and number density,²¹ which are called the effective scatterer size and number density. Considering that the effective scatterer size influences the shape of the scattering spectra in these two parameters, the effective scatterer size can be estimated by comparing the shape of the scattering spectra obtained with Mie theory. Then, using this effective scatterer size estimate, the effective number density can be calculated by comparing the magnitude of the scattering spectra with Mie theory.

2.3.

Image Reconstruction of Absorption and Scattering Coefficients

Reduced scattering coefficients were obtained with the clinical imaging system in use at Dartmouth (see Ref. 44 for details). Briefly, the system measures diffusely transmitted amplitude and phase shift at $100\mathrm{MHz}$ through up to $12\mathrm{cm}$ of tissue and is designed for breast imaging for characterization of tumors noninvasively. The measurements are taken at six wavelengths (661, 761, 785, 808, 826, and $849\mathrm{nm}$ ), the data is calibrated by a homogeneous phantom, and then fit to a finite element diffusion theory prediction of the data. The finite element model is iteratively fit to the measurements through a Newton minimization method, and images of absorption and reduced scattering coefficient are recovered at all wavelengths. In recent years, a direct reconstruction approach has been implemented that uses the known spectrum of the absorbing species in the tissue (hemoglobin, oxyhemoglobin, water) and the assumed model of scattering shown in Eq. 4. This inversion algorithm is described in detail in previous papers and is believed to provide more accurate prediction of the scattering parameters $A$ and $b$ from Eq. 4.^{45, 46} However, related to the fitting for particle size and density, it must be clearly understood that a sparse number of wavelengths cannot be used to accurately fit a highly detailed scattering spectrum, and in fact, oscillations in the spectrum have not been observed to date in bulk tissues. Thus the limited number of wavelengths is not sufficient for identifying oscillations in the scattering spectra that might be observed from larger Mie scattering particles, and these have not been observed with any certainty in diffuse spectra. However, oscillation in the diffuse spectra can be omitted if the dominant scattering particle sizes are thought to be smaller than $1000\mathrm{nm}$ .²¹

2.4.

Effective Scatterer Size and Number Density Estimation Method

To estimate effective scatterer size and number density, with the limited number of wavelengths available, any oscillation in the spectra must be omitted because the dominant scatterer sizes expected are small. The six experimental reduced scattering coefficients are fit to predictions of the power law expression in Eq. 4, first to get the scattering amplitude and scattering power, and then to generate a function of reduced scattering coefficient as a function of wavelength from 650 to $1000\mathrm{nm}$ . To extract the average particle size, this reduced scattering coefficient function is normalized at one wavelength $\left(800\mathrm{nm}\right)$ . This normalization eliminates the influence of number density, so that only average particle size needs to be fit in the first stage. Using a least-squares minimization, the average particle size can be effectively estimated. Then the number density can then be estimated by comparing the original unnormalized data with Mie theory [Eq. 2] using a second least-squares minimization method. In this fitting, the following error function is used:

2.5.

Effective Scatterer Size and Number Density Image Reconstruction Method

From the estimation method briefly presented above, the effective scatterer size, $⟨a⟩$ , and number density, $N$ , images are estimated from experimental measurements of the scattering power $b$ and amplitude $A$ . Data in Ref. 21 showed that out of $⟨a⟩$ and $N$ , it is only $⟨a⟩$ that influences the shape of the scattering spectra in the NIR region. Effective scatterer size $⟨a⟩$ can therefore be directly estimated from the experimental value of $b$ . When the scattering amplitude is normalized to $A=1$ at $800\mathrm{nm}$ , the effective number density $N=N_0$ is also, such that $N_0$ is also determined only by the scattering power $b$ . Then for arbitrary $A$ , the effective number density can then be estimated by $N=A{N}_0$ . So given $A$ and $b$ , both $⟨a⟩$ and $N$ can be estimated with this sequential process.

By simulating a series of $b$ values to get the results for $⟨a⟩$ and $N_0$ , an empirical polynomial power function was found for $⟨a⟩$ as a function of $b$ , and for $N_0$ as a function of $b$ . The equations for $⟨a⟩$ and $N_0$ are represented as polynomial expressions, listed in Table 1 . The fit to the data for $⟨a⟩$ was found to be optimal with a seventh order polynomial, whereas the fit to the data for $N_0$ required only a third order polynomial.

Table 1

Power function estimates for effective scatterer size ⟨a⟩ and normalized number density N0 , given the index of refraction values specified for outside (n1) and inside (n2) the spheres. For ⟨a⟩ , the power series was taken out to seven orders, whereas for N0 only three orders were required.

2.6.

Phantom Study

Data from the tomography system was used to validate this estimation approach by initially using tissue-simulating phantoms with well-characterized particle size and density. Heterogeneous gelatin phantoms were used, containing Intralipid inclusions at five different concentrations, ranging from 0.75 to 1.75% in steps of 0.25%. This range of inclusion concentrations provided a reduced scattering coefficient spectrum with the same particle size but linearly increasing number density. The data at each concentration was analyzed with the method described above with image reconstruction for the effective scatterer size and number density. By varying the concentration, each spectrum is expected to result in the same $⟨a⟩$ , with a theoretical value near $97(\pm 3)\mathrm{nm}$ , as measured by van Staveren ²⁵ This experiment provided an effective way to examine imaging inclusions within an otherwise homogeneous domain.

2.7.

Image Reconstruction of $⟨a⟩$ and $N$ for Human Breast Tissue

In the second step, the image reconstruction method was applied to cancer breast tissue images. The study was approved by the institutional committee for the protection of human subjects. NIR imaging studies were performed on symptomatic women recruited into the study following an abnormal mammogram and at the time of referral to biopsy. Informed consent was provided prior to the NIR imaging exam. Data from 12 patients were used in the analysis presented here, including 3 cancer patients and 9 benign tumor patients. For the diseased breast, both the cancer tissue and the background normal tissue were studied; the region of the tumor was delineated manually based upon a radiologist’s prediction of the location of the lesion from the x-ray mammograms. The regional average value was calculated for both the tumor and the normal tissue values. The contralateral breast was also analyzed separately as a whole value. Reconstruction of images for all patients was done with a direct estimate of the scattering parameters based upon an algorithm that incorporates all six wavelength data sets simultaneously and constrains the fitted images to recover hemoglobin, oxygen saturation, water fraction, and scattering spectra, fitted to $A$ and $b$ parameters.

2.8.

Study of Subject Images from the MRI-NIR System

In a final stage in the testing of the parameter estimation approach, the data processing method was applied to images from a magnetic resonance imaging (MRI)-coupled NIR system for breast tissue imaging. A clinical prototype hybrid imaging system, which combines NIR spectral tomography and MRI, was used to determine tissue structure with improved spatial resolution. ^{47, 48, 49, 50} Combining these two methodologies into a platform for simultaneous data acquisition allowed for excellent coregistration and data synergy. NIR reconstructions that use appropriate spatial constraints from a priori structural information about the boundary between adipose and fibroglandular tissue, lead to images with millimeter resolution and better represent the physiological differences between these two types of tissue. Thus, these images were used with the algorithm to predict effective scatterer size and number density to examine the typical values for adipose and fibroglandular tissue.

The imaging system used, described in detail by Brooksby, ⁴⁸ records measurements of NIR light transmission through a pendant breast in a planar, tomographic geometry. A portable cart houses the light generation and detection hardware subsystems. Six laser diodes (660 to $850\mathrm{nm}$ ) are amplitude modulated at $100\mathrm{MHz}$ . The bank of laser tubes is mounted on a linear translation stage, which sequentially couples the activated source into 16 bifurcated optical fiber bundles. The central seven fibers deliver the source light while the remaining fibers collect transmitted light and are coupled to photomultiplier tube (PMT) detectors located in the base of the cart. For each source, measurements of the amplitude and phase shift of the $100\mathrm{MHz}$ signal are acquired from 15 locations around the breast. The optical fiber bundles extend $13\mathrm{m}$ into a $1.5\mathrm{T}$ whole body MRI (GE Medical Systems), and the two data types (i.e., NIR and MRI) are acquired simultaneously. The patient lies on an open architecture breast array coil (MRI Devices) that also houses the MRI-compatible fiber positioning system. The plane of fibers spanning the circumference of a pendant breast can be positioned manually from nipple to chest wall if multiple planes of NIR data are desired.

Experience has shown that significant improvement in the stability and accuracy of the reconstruction process can be obtained by incorporating prior anatomical information as an input in the NIR parameter estimation problem. ^{48, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59} MRI priors are implemented through the incorporation of a regularization matrix, which has a Laplacian filter shape within the same tissue type, and effectively relaxes these smoothness constraints at the interface between different tissues, in directions normal to their common boundary.⁵⁰ This approach provides an inherent smoothing of a given tissue type, such as adipose or glandular tissues, but allows the two tissue types to be significantly different from one another. Based upon a series of phantom studies, this approach, combined with direct spectral reconstruction, has been shown to be the best estimate of the chromophore concentrations and scattering parameters of these normal breast tissues. The values of $A$ and $b$ in these images were then used to generate images of effective scatterer size and number density.

3.1.

Empirical Equations for Estimation of $⟨a⟩$ and $N$

In the imaging problem, the use of the empirical power functions listed in Table 1 were ideal to map the reconstructed images of scattering parameters onto the effective scatterer size and number density images. In all these studies, the size distribution function was assumed to be an exponentially decreasing function, and their corresponding relative refractive indexes were used. For the phantom study, $n_1=1.33$ , $n_2=1.47$ were used,²⁵ and for the cancer patient study, $n_1=1.36$ , $n_2=1.4$ were used.²² For the MRI-NIR patient studies, where the tissues could be uniquely segmented between adipose and fibroglandular tissues, $n_1=1.36$ and $n_2=1.4$ were used for the fibroglandular tissue, with the latter value being taken for bilipid membranes,²² while for the adipose tissues, $n_1=1.36$ and $n_2=1.49$ were used with the latter value being the measured value for the lipid.⁵⁰ The corresponding empirical equations for these three sets of refractive indexes were estimated and the coefficients are shown in Table 1. The fitting figures for the $n_1=1.36$ , $n_2=1.4$ set are shown in Fig. 1 .

Fig. 1

The relationship between scatter power, $b$ , and effective scatterer size, $⟨a⟩$ , is shown in (a) for simulated data (points) and the best fit polynomial curve (coefficients in Table 1). The relationship between effective normalized number density (the effective number density when $A=1$ ) and scatter power is shown in (b) for the same parameters, with the fitted polynomial curve shown (coefficients in Table 1).

3.2.

Phantom Study Results

Figure 2a shows the diagram of the phantom, and 2b shows the reconstructed images for $b$ (Scatt. Power) and $A$ (Scatt. Ampl.), from the phantom with an inclusion of Intralipid at concentration of 1.0%. From these two images and using the empirical equations, the reconstructed images of effective scatterer size (Effec. Size) and number density (Num. Density) were obtained as shown in 2b.

Fig. 2

In (a), a diagram of the phantom geometry is shown, and in (b), a set of reconstructed images of the phantom with 1.0% Intralipid concentration is shown. The background medium was a solid resin phantom used in many previous studies (Refs. 46, 62).

The physical region of interest (ROI) was taken as the size of the hole, and the average value over the ROI was estimated. For example, for the phantom with an Intralipid concentration of 1.0% in the inclusion, the effective scatterer size and number density average values in the ROI were $111\mathrm{nm}$ and $2.64\times {10}^{14}{\mathrm{m}}^{-3}$ .

The summary values for each of the five phantom concentrations of Intralipid is shown in Fig. 3 , with 3a showing the effective scatterer size. The circles indicate the results for the phantom data. The estimation results for the phantoms with different intralipid concentrations were almost all the same, with an average value of $106\mathrm{nm}$ , which is shown by the solid line in Fig. 3a, and standard deviation was $4\mathrm{nm}$ . This result agrees to within 9% of the expected value of $97\pm 3\mathrm{nm}$ from van Staveren, as determined by electron microscopy studies in their paper.²⁵

Fig. 3

The results are shown for different intralipid concentrations of (a) average effective scatterer size, $⟨a⟩$ , and (b) average effective number density, $N$ , using bulk values from within the inclusion region, as shown in Fig. 2.

Figure 3b shows the result for the ratio of effective number density to concentration, for varying concentrations of the Intralipid inclusion. The circles indicate the results, with the ratio of $N$ over concentration being nearly constant for these phantoms, having average and standard deviation values of $2.4\times {10}^{16}{\mathrm{m}}^{-3}\pm 0.7\times {10}^{16}{\mathrm{m}}^{-3}$ , which are shown by the solid line in Fig. 3b. Because these phantom results were consistent with what was expected, the approach was utilized on human data, which is discussed in Sec. 3.3.

3.3.

Patient Imaging Study Results

Figure 4a is a representative image example from a patient with an invasive ductal carcinoma tumor in the breast, showing an en face or craniocaudal slice of the breast as viewed toward the patient. From our reconstruction, the hemoglobin $\left({\mathrm{Hb}}_{\mathrm{T}}\right)$ , oxygen saturation (Oxy.), water fraction (Wat.), scattering power (Scatt. Power), and scattering amplitude (Scatt. Ampl.) images were recovered, and these are all shown in Fig. 4a in the plane of the tumor. Using the two empirical equations for $n_1=1.36$ , $n_2=1.4$ , the effective scatterer size and number density values were recovered at each pixel in the image, using the images of the scattering power and scattering amplitude. The resulting images are also shown in Fig. 4a. These images indicate that the cancer tissue has smaller effective scatterer size and larger effective number density than the background in this case. Figure 4b is a representative image example from a patient with a benign tumor in the breast.

Fig. 4

In (a), reconstructed NIR images are shown for a cancer patient, showing the plane of the tumor in the breast, sliced in a craniocaudal view. The panel of images shows total hemoglobin concentration $\left({\mathrm{Hb}}_{\mathrm{T}}\right)$ , oxygen saturation (Oxy.), water fraction (Wat.), scatter power (Scatt. Power) and amplitude (Scatt. Ampl.), as well as effective scatterer size (Effec. Size), and number density (Num. Density) images. In (b), the reconstructed NIR images for a patient with a benign tumor in the breast are shown.

The ROI value is the average over the tumor area. The average value over the remaining breast tissue excluded the tumor ROI. The average value for the contralateral breast ROI was estimated at the same mirror location as in the symptomatic breast. Studies have shown that breast tissue is highly symmetric between breasts, such that comparisons between ROI values taken from mirror locations in the symptomatic and normal breasts are a good way to see differences from the normal condition.

Figure 5 is the summary data for the patients, including the cancer patients and the benign tumor patients. Figures 5a and 5b show the average effective scatterer size and number density values for the ROI, the background and contralateral breast over all the cancer patients included in this study, as well as the benign tumor patients. For the average effective scatterer size in the ROI, the $t$ -test value between of the cancer patients and the benign patients is 0.025. The data in Figs. 5a and 5b indicate that for the cancer patients, the ROI has smaller effective scatterer size and larger effective number density compared with both the background and the contralateral breast; while for the benign tumors, effective scatterer size and number density are all relatively close to that of the background and the contralateral breast values. But it needs to be mentioned that there are too few cancer patients to obtain a generic conclusion about tumors in this study. More patients will be studied in the future to reach such a conclusion. Figures 6a and 6b show the average difference over each kind of patient for effective scatterer size $⟨a⟩$ and number density $N$ between the ROI and background and between the ROI and the contralateral breast.

Fig. 5

The values for groups of subjects in terms of cancer and benign lesions are shown, for (a) average effective scatterer sizes and (b) average effective number density.

Fig. 6

The relative differences in average region values are shown for cancer and benign lesions for (a) average effective scatterer sizes and (b) average effective number density.

3.4.

MRI-NIR Subject Study Results

In the final stage of examining the predictions for effective scatterer size and number density values, the images reconstructed with combined MRI-NIR were used.⁶² The results for two MRI-NIR patients are shown in Figs. 7a and 7b . Figure 7a (Region) is the MRI image for the first patient. The inner part of the breast (darker in the image) is the fibroglandular tissue, and the outer part (gray in color) is the adipose or fatty tissue. The hemoglobin $\left({\mathrm{Hb}}_{\mathrm{T}}\right)$ , oxygen saturation (Oxy.), water fraction (Wat.), and scattering power (Scatt. Power) and amplitude (Scatt. Ampl.) images are shown in Fig. 7a. In the estimation of $⟨a⟩$ and $N$ images, separate values for $n_2$ were used for each of the two tissue types, as discussed in Sec. 2, and images of the breast interior were recovered. The estimated images for $⟨a⟩$ (Effec. Size) and $N$ (Num. Density) are also shown in Fig. 7a. These images indicate that the fibroglandular tissue has smaller effective scatterer size and larger effective number density than the fatty tissue. For the second patient, the same conclusion is obtained from these results.

Fig. 7

Reconstructed images for two of the MRI-NIR subjects (a,b) are shown, with the MRI image in (Region), and hemoglobin $\left({\mathrm{Hb}}_{\mathrm{T}}\right)$ , oxygen saturation (Oxy.), water fraction (Wat.), scatter power (Scatt. Power), and scatter amplitude (Scatt. Ampl.). Estimates of $⟨a⟩$ and $N$ are shown in (Effec. Size) and (Num. Dens.), as calculated from (Scatt. Power) and (Scatt. Ampl.). Image (b) is predominantly fatty tissue (gray in MRI image), whereas (a) is predominantly fibroglandular tissue (darker in MRI image).

In total, 10 sets of images from the MRI-NIR subjects were used in this evaluation and the summary result for these is shown in Fig. 8 . The summary results show that the fibroglandular tissue has smaller effective scatterer size and larger effective number density than the adipose tissue. For the average effective scatterer size, the $t$ test between of the fibroglandular and adipose tissue showed a significant difference with a $p$ value of 0.048. For the number density, the $t$ test showed a significant difference between the fibroglandular and adipose tissue with a smaller $p$ value of 0.004. Thus NIR tomography appears to be able to delineate fibroglandular tissue from fatty tissue more by the number density than by the average size.

Fig. 8

The average tissue values are shown pooled for all MRI-NIR subject images, for effective scatterer size (a), and number density (b), for adipose and fibroglandular tissues.