1.

Introduction

Near-infrared (NIR) diffuse optical tomography is a noninvasive imaging technique used to obtain pathophysiological images of biological tissue. Using this technique, NIR light between wavelengths of 650 and $930\mathrm{nm}$ is injected into the external surface of tissue and the emergent light is measured along the tissue boundaries. These “boundary” data are then used with a light propagation model to predict the internal chromophore concentrations of soft tissue such as hemoglobin, oxygen saturation, and water content, as well as scattering properties.^{1, 2, 3, 4} Although NIR light suffers from high scattering within tissue at these wavelengths, diffuse NIR optical imaging has shown its potential as a tool for quantifying brain function activity.⁵ as well as breast cancer characterization and detection.^{1, 6, 7, 8} The issue of which spectral bands and measurement data types to use has been a point of continuous debate in the field, and our study applies the original and unique method presented by Corlu ⁹ to broadband data selection to examine which spectral windows present the unique data needed for spectral image reconstruction.

There are three common methods used for the collection of NIR transmission data: continuous wave, time domain, and frequency domain.¹⁰ Of these methods, continuous wave (cw) systems are the most cost effective and easiest to implement, where only the intensity of the emergent NIR light is measured at a given wavelength. Using cw systems, no information about the path length of NIR light within tissue is measured. Time-resolved systems are the most expensive, since they rely on photon counting systems and suffer from long data integration times,¹¹ but can potentially provide absolute information regarding the total signal intensity as well as photon propagation time by means of measuring the temporal point spread function (TPSF). Alternatively, frequency domain systems take advantage of amplitude modulated signals to measure both the intensity and phase of the boundary data, providing information regarding the average path length of the measured photons in tissue.¹² Although it has been shown that the absorption and scattering properties at a given wavelength cannot be separated using continuous wave data alone without a priori information,¹³ recent work in spectral imaging has shown that by using spectral constraints using cw data, it may be possible to resolve chromophore concentrations directly, using an optimum set containing a finite number of wavelengths.⁹ This original work by Corlu ⁹ demonstrated that by using spectral reconstruction and wavelength optimization, it is possible to overcome the nonuniqueness problem observed in cw measurements.

In this work, a method is presented to examine the viability of using large cw datasets containing many wavelengths to achieve a unique solution with minimal crosstalk between chromophores. Other tomographic studies to date have restricted spectral reconstruction to using datasets of five or six wavelengths for image reconstruction due to technical limitations.^{14, 15} Although a set of five or six wavelengths will be computationally faster, the use of more wavelengths increases the amount of information available and is considered important in constraining the inverse problem further, such that crosstalk between chromophores can be reduced, which can be especially important for systems or wavelengths (such as 650 or $950\mathrm{nm}$ , where absorption by deoxyhemoglobin and water dominate) with higher measurement noise. The measurement of cw data over the entire range of NIR wavelengths is currently limited using most experimental systems available. The use of two or more wavelengths has been used within a time domain system,¹⁶ as well as five or six discrete wavelengths using frequency domain systems.^{7, 11} More recently, a cw system allowing the complete spectral measurements of NIR data for small animal imaging system has also been developed,¹⁷ allowing dual modality MRI and NIR data collection for the study of small animal brain physiology. However, to date, there has been little work into the investigation of an optimum set of wavelengths to provide a unique and well-conditioned solution, under the assumption that boundary data are available at all wavelengths between 650 and $930\mathrm{nm}$ .

The methods proposed by Nielsen and Brendal¹⁸ and Corlu ⁹ are used as the basis of the presented study. Rather than restricting the image reconstruction to the best five or six wavelengths that provide a better conditioned problem, a method is proposed that finds ranges, or “windows,” of wavelengths over the entire spectrum, which contribute to producing highly conditioned and unique images.

2.

Theory

For model-based image reconstruction, an accurate forward model is required that predicts the path of photons within tissue and a reliable inverse model, which gives a unique set of optical parameters based on some boundary measurements. The transport of light through tissue can accurately be described by the diffusion approximation to the radiative transport equation when the detectors are placed more than three scattering lengths away from the sources, and where light propagation is dominated by scattering events.² For continuous wave light, the diffusion equation can be written as

Given a set of optical parameters, the diffusion equation can be used to calculate light distribution within a domain and hence the boundary data. For complex geometries, the diffusion equation is difficult to solve analytically, so numerical solutions such as the finite element method (FEM) can be used as a flexible approach to the forward problem.¹

To solve the forward problem, the distribution of the optical parameters ${\mu }_a$ and $\kappa$ must be known. However, a multispectral model can relate the wavelength dependence of the optical parameters with concentrations of oxygenated hemoglobin $\left(\mathrm{Hb}{\mathrm{O}}^2\right)$ , deoxygenated hemoglobin (Hb), water $\left({\mathrm{H}}_2\mathrm{O}\right)$ fraction, and scattering properties.^{3, 9} The absorption coefficient can be calculated from the chromophore concentrations $\left(c\right)$ and published extinction coefficients $\left(\epsilon \right)$ .¹⁹ From Beer’s law, the absorption at a given wavelength is given by a linear summation of the contributions from each absorber, which has a known extinction spectra, and the concentration $c$ is to be estimated from the absorption spectra,

The inverse problem has the aim of recovering the distribution of optical parameters in tissue given a set of boundary data at a given wavelength. This is achieved by reducing the difference between the measured data ${\Phi }^M$ and the calculated data (from the forward model) ${\Phi }^c$ . The image reconstruction algorithm minimizes this objective function,

For a multispectral model, the image reconstruction is also an iterative problem and a modified version of Eq. 5 is used. However for multispectral reconstruction, the Jacobian is a matrix consisting of a number of measurements by a number of wavelength rows and number of FEM nodes by number of chromophores columns.¹⁴ Therefore by dimensionality, the update data vector has a length of $NM$ by number of wavelengths, leading to the update vector $\delta \mu$ corresponding to the total number of unknowns.

It has been shown that continuous wave measurements can result in identical data for different sets of optical parameters due to crosstalk between absorption and scattering. But the correct wavelength selection in a multispectra model can reduce the effect of this uniqueness problem.¹³ by directly reconstructing for the chromophore and scattering properties. Algorithms have been suggested that extend the image reconstruction process by choosing sets of optimized finite wavelengths for a continuous wave system¹³ to reduce the nonuniqueness from image reconstruction and reduce the crosstalk between the chromophores. A further criteria has also been suggested that considers the error in the assumed extinction coefficients of the chromophores and examines its role in the optimized wavelength set.¹⁸ However, these methods have been limited to small and finite wavelength sets (typically five wavelengths), restricting its use for modified systems where larger wavelength sets may be available, and hence allowing a better conditioning of the inverse problem.

Since different sets of optical parameters have been shown to give identical cw data due to crosstalk between absorption and scattering, the nonuniqueness concept for cw measurements,¹³ proposes that by solving a matrix expression that considers the spectral response of the chromophores directly, the ranges of wavelengths that provide unique measurement data can be found using:

The inverse problem must also be able to distinguish between chromophores, such that each chromophore has an equal contribution to the total measured absorption that can be derived from expressing Eq. 2 in matrix form,

Assuming that only six discrete wavelengths can be measured from the entire spectral range of $650\text{to}930\mathrm{nm}$ , with a separation of $4\mathrm{nm}$ (from a complete set of 71 wavelengths), there exists $1.4\times {10}^8$ combinations $\left(C_6^{71}\right)$ of six wavelengths sets. Each set has a condition number and associated residual calculated using published values of extinction and scattering coefficients.²⁰ In this work, three chromophores have been assumed to vary (hemoglobin, oxyhemoglobin, and water) as well as the scatter amplitude, while scatter power is held constant and is equal to 1. A scatter plot of the calculated residual and condition number for each of these $1.4\times {10}^8$ combination six wavelength sets have been calculated and are shown in Fig. 1 . Due to the large number of points, the results shown are limited to those sets with a residual above 0.35 and a condition number below 400. The sets of wavelengths of interest that provide the most adequate information about the problems are those in the proximity of the top left-hand corner (as highlighted by the shaded region in Fig. 1), as these contain the sets that have the highest residual and lowest condition number.

Fig. 1

Scatter plot of the residual and condition number calculated for three chromophores, hemoglobin, deoxyhemoglobin, and water. Each point represents a different combination set of six wavelengths from the spectra $650\text{to}930\mathrm{nm}$ with $4\text{-}\mathrm{nm}$ separation. The region highlighted satisfies the criteria of uniqueness (large residual) and chromophore separability (small condition number).

Figure 2 shows a histogram of the distribution of optimal wavelength sets with high residual $(>0.47)$ and low condition number $(<200)$ criteria, as taken from the region highlighted in Fig. 1 as being those that best fulfill the criteria. The histogram is formed by counting the frequency where a particular wavelength appears for each of the six wavelength choices (or bins), normalized to the most frequently appearing wavelength for that specific bin. From Fig. 2, it is evident that six individual wavelengths are seen as having the highest count; however, these wavelength peaks themselves do not form the best optimized set, since it lies in the bottom right-hand corner of the highlighted region in Fig. 1. There are four clear spectral regions (or windows) that contribute to the high residual and low condition criteria, and these are distributed around 650, 736, 874, and $930\mathrm{nm}$ . This shows that the wavelengths within the optimized sets that meet the residual and condition calculation are confined to narrow spectral windows. For each optimized set, the wavelengths 650 and $930\mathrm{nm}$ are common, as there is only a single peak with a normalized count of unity due to the peak in the absorption by Hb at $650\mathrm{nm}$ and water at $930\mathrm{nm}$ . The second, third, and fourth common wavelengths are found in the region of $706\text{to}738\mathrm{nm}$ , while the fifth most common is taken from the region $846\text{to}882\mathrm{nm}$ . However, it is important to note that increasing the number of wavelengths in a set (from six wavelengths to perhaps 10 or 12 wavelengths per set) does not increase the width of these bands significantly, and the trends are similar.⁹

Fig. 2

Histogram of optimized wavelength sets that satisfy the high residual and low condition number criteria for three varying chromophores. Each wavelength count is normalized to the most frequent wavelength.

For a spectral model, a larger number of wavelengths confines the inverse problem by reducing the ill-posed nature. Taking all the wavelengths within the spectra $650\text{to}930\mathrm{nm}$ with a separation of $4\mathrm{nm}$ , a residual of 2 and a condition number of 246 is calculated. The condition number is comparable to each optimized set of six wavelengths, but the residual is approximately four times higher, implying that using 71 wavelengths for the inverse problem could give a better solution. However, increasing the amount of data also increases the size of the Jacobian, and thus the memory requirements and the computation time.

A large wavelength set can also be derived from the histogram in Fig. 2. Each of the wavelengths within the histogram contributes to a small condition number or a large residual within the optimized sets. Due to computational limits, it is difficult to find optimized 10, 15, or 40 wavelength sets without increasing the step size dramatically (to 10 or $20\mathrm{nm}$ ), which will inevitably reduce the spectral resolution from the absorption spectrum. Instead, by taking each wavelength that contributes to the residual and condition criteria, as shown in the histogram of Fig. 2, a wavelength set using windows of the spectrum can be formed. This set, in this instance, consists of 21 wavelengths, but importantly is only 28% of the size of the original full spectrum considered. The condition number is comparable to each of the other wavelength sets, as shown by Table 1 .

Table 1

Residual, average residual (for the number of wavelength used), and condition number for corresponding wavelength sets.

3.

Methods and Results

To compare the effects of wavelength selection, simulations were performed using a simulated tissue region. Using a known and defined numerical model has the benefit that the target locations and values are known exactly, allowing accurate analysis of the reconstructed images, and also allowing the effects of the wavelength selection to be investigated exactly. The simulated 2-D region consisted of a uniform circular mesh of radius $43\mathrm{mm}$ with 1785 nodes corresponding to 3418 linear triangular elements. 16 optical fibers placed equidistant from the center of the mesh were used for the data collection of amplitude with a modulation frequency of $100\mathrm{MHz}$ . Random noise of either 1, 5, or 10% is then added systematically to the cw intensity-only data. The background had scattering amplitude equal to $1{\mathrm{mm}}^{-1}$ , water content of 40%, deoxygenated hemoglobin concentration of $0.01\mathrm{mM}$ , and oxygenated hemoglobin concentration of $0.01\mathrm{mM}$ . Four circular anomalies of radius $7.5\mathrm{mm}$ are placed within the model in separate individual locations within the mesh, as shown in Fig. 3 (left-hand column). Within each anomaly, only one parameter was changed, corresponding to either scattering amplitude of $0.5{\mathrm{mm}}^{-1}$ , water content of 80%, deoxygenated hemoglobin concentration of $0.02\mathrm{mM}$ , or oxygenated hemoglobin concentration of $0.02\mathrm{mM}$ . The scatter power was set as a constant throughout the entire phantom at $1.0$ . To reduce the number of unknown parameters, a region mapping algorithm was used to transform the Jacobian from a nodal basis to a region basis.²¹ In this case, by using. “hard a priori” and reconstructing the values only for the regions, the differences in the reconstructed images will be due to the different selection of wavelengths alone and not due to the inverse problem. A transformation matrix $k$ that has dimensions of a number of unknowns (number of FEM nodes times number of chromophores) by a number of regions by a number of chromophores is applied to the Jacobian such that,

Fig. 3

Region-based spectral reconstruction of scattering amplitude, water, deoxyhemoglobin, and oxyhemoglobin using three different wavelength selections; the full spectrum $\left(650\text{to}930\mathrm{nm}\right)$ , denoted by every $4\mathrm{nm}$ , windows of wavelengths sampled from the full spectrum, and a set of six wavelengths derived from a typical measurement system. The cw data has 5% added noise.

The reconstructions were carried out using three different wavelength sets: 1. the full spectrum consisting of the range $650\text{to}930\mathrm{nm}$ with $4\text{-}\mathrm{nm}$ separation; 2. the windows method consisting of the wavelength ranges $650\mathrm{nm}$ , $706\text{to}738\mathrm{nm}$ , $846\text{to}882\mathrm{nm}$ , and $930\mathrm{nm}$ at $4\mathrm{nm}$ separation; and 3. a set of six wavelengths determined by a current realistic system using laser diodes at 660, 734, 760, 808, 826, and $850\mathrm{nm}$ , which is clearly not an optimized set as described before, but limited by the typical availability of specific laser diodes.¹⁴

Figure 3 (columns 2, 3, and 4) shows the reconstructed images to demonstrate the crosstalk for each of the three wavelength methods for 5% added noise. Similar results are found for the crosstalk between chromophores for 1 and 10% noise (not shown). The corresponding reconstructed and target values for each image are shown in Table 2 . Reconstructions were also carried out using different background values such as high oxyhemoglobin with a low target value. The results in these cases were similar to those presented here, and are not shown.

Table 2

Reconstructed anomaly values for the four regions of interest from a simple phantom using three different wavelength selections with 1, 5, and 10% added noise. Reconstructions using the whole spectrum with a 4-nm separation (denoted by 71 wavelengths), the windows method (denoted by 21 wavelengths), and a selection of six wavelengths.

To test the robustness of the windows method, standard image reconstructions using a spectrally modified version of Eq. 5 have been carried out on the same simulated phantom as defined before, since in most real situations the exact locations of regions are not known. Reconstructions for 1, 5, and 10% added noise are shown in Figs. 4, 5, 6 respectively.

Fig. 4

Spectral image reconstruction of scattering amplitude, water, deoxyhemoglobin and oxyhemoglobin using a simulated phantom with 1% added noise for the three wavelength selections; the full spectrum, the windows sampled from the full spectrum, and a set of six wavelengths.

Fig. 5

Same as Fig. 4, but with 5% added noise.

Fig. 6

Same as Fig. 4, but with 10% added noise.

To test the feasibility of the windows method for wavelength selection for an experimental system, data have been used in reconstruction of a simple gelatine-based phantom. The cylindrical tissue simulation was designed with radius $43\mathrm{mm}$ and length $88\mathrm{mm}$ . A cylindrical anomaly of radius $12.5\mathrm{mm}$ is placed $25\mathrm{mm}$ from the center. The background hemoglobin concentration is $0.011\mathrm{mM}$ , while the anomaly has hemoglobin concentration of $0.027\mathrm{mM}$ . Wavelengths between 690 and $840\mathrm{nm}$ with $5\text{-}\mathrm{nm}$ separation, excluding $760\mathrm{nm}$ , were collected, giving a total of 29 wavelengths using a frequency domain system incorporated with a mode-locked Ti:sapphire laser for the source between 690- and $850\text{-}\mathrm{nm}$ wavelengths, described elsewhere.²²

Using the same method as previously defined with the new wavelength range, windows of wavelengths of $695\mathrm{nm}$ , $715\text{to}735\mathrm{nm}$ , $755\mathrm{nm}$ , $765\text{to}770\mathrm{nm}$ , and $800\text{to}810\mathrm{nm}$ can be calculated as the wavelength set (12 wavelengths in total). The residual per wavelength increases from 0.0065 using the whole spectrum to 0.0122 for the windows method, while the condition number for the windows method is 654 compared to 630 using the entire available spectrum. The inverse problem is carried out on the same circular mesh as before, after appropriate data calibration,¹ and Fig. 7 shows the reconstructed images with the full spectrum $\left(695\text{to}840\mathrm{nm}\right)$ , denoted by 29 wavelengths, the windows method, and a selection of five wavelengths from across the spectrum (700, 750, 800, 830, and $840\mathrm{nm}$ ).

Fig. 7

Experimental phantom images reconstructed with three different wavelength selections using cw data and reconstructing for scatter amplitude, total hemoglobin, oxygen saturation, and water.

4.

Discussions

Using a numerical model and region-based image reconstruction, the use of windows-based wavelength optimization was investigated to allow a direct analysis of results without image inaccuracy due to the inverse problem. Region-based images were reconstructed using either the complete spectrum at $4\text{-}\mathrm{nm}$ separation, a set of windows-based optimized wavelengths, or six defined wavelengths, as commonly available in most experimental systems, Fig. 3, and Table 2. For each method the target values were observed, and additionally, to allow a more quantitative analysis, error was also calculated by

Table 3

Calculated error (%) for the anomalous region (A) and background (BG) using two different reconstruction methods and three different wavelength selections with 5% added noise. A similar trend is seen at 1 and 10% added noise. Reconstructions using the whole spectrum with a 4-nm separation (denoted by 71 wavelengths), the windows method (denoted by 21 wavelengths), and a selection of six wavelengths.

To evaluate the use of optimized wavelengths in a true tomographic setting, the data as used in the previous case were used to reconstruct spatially varying chromophore concentrations throughout the model (Figs. 4, 5, 6 and Table 3). With 5% added noise (Fig. 5), using the full spectrum gives errors in the reconstructed anomalies of between 20% $\left({\mathrm{H}}_2\mathrm{O}\right)$ and 33% $\left(\mathrm{Hb}{\mathrm{O}}_2\right)$ . With the windows method the error ranges from 18% $\left({\mathrm{H}}_2\mathrm{O}\right)$ to 30% $\left(\mathrm{Hb}{\mathrm{O}}_2\right)$ . With six wavelengths, the range in error increases to between 24% (scattering amplitude) and 35% $\left(\mathrm{Hb}{\mathrm{O}}_2\right)$ . The crosstalk also varied between methods. In each case the largest error was found in the background of the scattering amplitude at 7.1% using the full spectrum, 6.8% using the windows method, and 9.8% using six wavelengths. Using just six wavelengths was shown to be both quantitatively and qualitatively less accurate for all noise levels modeled, although there could be two reasons for this. First, the six wavelengths are based on a frequency modulated system currently in use for the collection of clinical data due to availability of laser diodes, and the selection is therefore not optimized. Second, the number of wavelengths is far fewer and the amount of information available is therefore limited. However, there are also computational benefits for choosing a smaller wavelength set than the full spectral range. The size of the Hessian matrix $\left({\mathbf{J}}^T\mathbf{J}\right)$ depends on the number of wavelengths chosen and has a computational limit dependent on the system. The size of the Hessian also has a direct impact on the computation time of the inverse problem. If the Hessian is of the size $N\times N$ , the number of calculations is a function of $N^3$ and the memory of $N^2$ .²³ The simulated phantom uses 16 sources and 15 detectors to model the data, which gives 240 measurements at each wavelength. Using the full spectrum (71 wavelengths) requires $2.3\mathrm{Gb}$ to store the Hessian matrix, whereas the windows method (21 wavelengths) requires only $203\mathrm{Mb}$ . The windows method is computationally more efficient than that of the full spectrum and more desirable due to the reduction in computation time.

To validate the results, a set of experimentally measured data from a gelatine-based phantom was used and images were reconstructed (Fig. 7). The reconstructed images using five wavelengths show an overestimate in the total amount of Hb by 30%, while the number of artifacts is also high for $\mathrm{S}{\mathrm{O}}_2$ , ${\mathrm{H}}_2\mathrm{O}$ , and scattering amplitude. This shows that a lack of information or nonoptimized wavelength sets do not give accurate results. The reconstructed images using the full spectrum have reconstructed the total Hb qualitatively well, while artifacts can be seen in the scattering amplitude. This could be because using the full spectrum in this case overconstrains the problem. The windows method reconstructed the peak in the Hb as well as reduced artifacts in the other images. Importantly, the windows method required $52\sec$ per iteration, while the full spectrum took $114\sec$ . Thus the window method is computationally more efficient.

5.

Conclusions

The use of large wavelength sets in NIR imaging can help optimize the inverse problem and is especially useful for systems where the magnitude of measurement noise is large. However, a large number of wavelengths also increases the size of the computational problem and thus the memory requirement, which becomes more important for large complicated models. Using an optimized set of a small number of wavelengths is computationally more efficient, and theoretically can provide unique images using continuous wave measurements. An objective method to allow the determination of a set of optimized wavelengths, which is a function of the chromophores being imaged, was presented. Using the criteria that maximizes this residual, while maintaining a low condition for the inverse problem, the number of wavelengths needed to provide adequate information for full spectral imaging can be reduced. Reconstructed images, based on either regional image reconstruction or those with tomographic images, show that by adequate optimization, the total number of required wavelengths can be reduced by more than 50%, without loss to image quality or accuracy, but with the added advantage of reducing both data collection and image reconstruction time. More importantly, using this optimized set of wavelengths, the total error within the reconstructed images, including the crosstalk between parameters, is reduced, as compared to the case where the entire spectrum is used. This conclusion is very important, suggesting that not all wavelengths within the NIR range are beneficial for spectral tomography. Future work will include the use of path-length information to allow the inclusion of scatter power within the wavelength optimization method, which will allow reconstructed images not only of absorption-based chromophores, but also scatter-dependent changes.