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 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.5 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 9 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.10 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,11 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.12 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,13 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.9 This original work by Corlu 9 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 , 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,16 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,17 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 .
The methods proposed by Nielsen and Brendal18 and Corlu 9 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.
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.2 For continuous wave light, the diffusion equation can be written asis the diffusion coefficient, , is the fluence rate, is the absorption coefficient, is the reduced scattering coefficient, , where is the transport scatter and is the anisotropic factor (set as 0.9 for soft tissue), and is an isotropic source modeled as one scattering distance inside the external boundary.
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.1
To solve the forward problem, the distribution of the optical parameters and must be known. However, a multispectral model can relate the wavelength dependence of the optical parameters with concentrations of oxygenated hemoglobin , deoxygenated hemoglobin (Hb), water fraction, and scattering properties.3, 9 The absorption coefficient can be calculated from the chromophore concentrations and published extinction coefficients .19 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 is to be estimated from the absorption spectra,is the number of chromophores and is the wavelength number.20 The reduced scattering parameter is given by the empirical formula, given by is the scattering amplitude and is the scattering power.4
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 and the calculated data (from the forward model) . The image reconstruction algorithm minimizes this objective function,is the optical parameters to be reconstructed and is the number of measurements. is known as the projection error and gives a value with which the convergence of the inverse model can be determined. However, the inverse problem is nonlinear and illposed due to the large number of unknowns. Assuming that the initial estimate of the optical properties is close to the actual optical properties, this inverse problem can be linearized by minimizing Eq. 4 with respect to the optical parameters and only keeping the first-order terms. This gives the iterative equation is the update optical parameter vector, is the update data vector, and is the regularization parameter scaled by the maximum of the diagonal of the Hessian and reduced systematically at each iteration. is the Jacobian matrix that relates the change in data to a change in optical parameters. In this work, Eq. 5 is solved iteratively until the projection error dose is reduced up to 0.2%, compared to the previous iteration.
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.14 Therefore by dimensionality, the update data vector has a length of by number of wavelengths, leading to the update vector 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.13 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 system13 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.18 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,13 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:is the figure of merit or norm of the residual that determines the uniqueness of this matrix problem, and is the extinction coefficient matrix normalized by the wavelength given by is the number of chromophores and is the total number of wavelengths considered in the image reconstruction. By the definition of Eq. 6, if two or more sets of parameters lead to the same solution, the inverse problem has a nonunique solution and the residual will be zero. Using this criterion, the inverse problem therefore has a unique solution if the residual is large. To achieve the uniqueness criteria as given by Eq. 6, if these three chromophores and the scattering amplitude are allowed to vary, then it has been suggested for cw data that at least four wavelengths must be used within a spectral model.9 For this analysis and the work presented, the scattering power is kept as a constant throughout the problem. This is because the separability of scatter power from other parameters is poor, and the crosstalk between chromophores is high using continuous wave data alone.9
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,8] has similar singular values, the condition number of the matrix is also small and each chromophore would contribute equally to the absorption.
Assuming that only six discrete wavelengths can be measured from the entire spectral range of , with a separation of (from a complete set of 71 wavelengths), there exists combinations of six wavelengths sets. Each set has a condition number and associated residual calculated using published values of extinction and scattering coefficients.20 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 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.
Figure 2 shows a histogram of the distribution of optimal wavelength sets with high residual and low condition number 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 . 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 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 and water at . The second, third, and fourth common wavelengths are found in the region of , while the fifth most common is taken from the region . 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.9
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 with a separation of , 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 ), 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 .
Residual, average residual (for the number of wavelength used), and condition number for corresponding wavelength sets.
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 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 . 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 , water content of 40%, deoxygenated hemoglobin concentration of , and oxygenated hemoglobin concentration of . Four circular anomalies of radius 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 , water content of 80%, deoxygenated hemoglobin concentration of , or oxygenated hemoglobin concentration of . The scatter power was set as a constant throughout the entire phantom at . 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.21 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 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,transforms the update vector back into a nodal basis. The Hessian was then much smaller (number of regions by number of reconstructed chromophores) which in this case is 20 (five regions and four chromophores).
The reconstructions were carried out using three different wavelength sets: 1. the full spectrum consisting of the range with separation; 2. the windows method consisting of the wavelength ranges , , , and at separation; and 3. a set of six wavelengths determined by a current realistic system using laser diodes at 660, 734, 760, 808, 826, and , which is clearly not an optimized set as described before, but limited by the typical availability of specific laser diodes.14
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.
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.
|Noise level||Number ofwavelengths||SA (mm−1)||HbO2 (mM)||Hb (mM)||H2O (%)|
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.
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 and length . A cylindrical anomaly of radius is placed from the center. The background hemoglobin concentration is , while the anomaly has hemoglobin concentration of . Wavelengths between 690 and with separation, excluding , 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 wavelengths, described elsewhere.22
Using the same method as previously defined with the new wavelength range, windows of wavelengths of , , , , and 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,1 and Fig. 7 shows the reconstructed images with the full spectrum , denoted by 29 wavelengths, the windows method, and a selection of five wavelengths from across the spectrum (700, 750, 800, 830, and ).
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 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 byis the calculated parameter and is the target value for each region of interest, summarized by Table 3 . Using region-based reconstruction, for both the windows method and the full spectrum the target values were found with an error less than 1.4%, with 5% added noise, whereas using six wavelengths the largest error was found with a maximum of 4.8%. However, the crosstalk between chromophores varied depending on the wavelength selection. Using both the windows method and the full spectrum gave the lowest crosstalk between each of the regions with the largest percentage error in the reconstruction of the background of 0.7 and 1.5%, respectively. However, using six wavelengths gave the largest crosstalk with the reconstructed background of 7%. The windows method has a much lower crosstalk at high noise in comparison to using just six wavelengths, and remains similar to the results of taking the full spectrum. As the noise increases, the qualitative accuracy attained by the reconstruction decreases for each wavelength set. However, even at 10% noise, the error in the target values for each chromophore is small for each selection of wavelengths, with a maximum error found for the reconstruction of the anomaly with six wavelengths of 6.1% and the anomaly with the windows method and the full spectrum of 6.8 and 3.9%, respectively. Importantly, the reconstructed values using the windows method are very similar to the reconstructed values found using the entire spectrum, and the crosstalk is much smaller than using the selection of six wavelengths, in Table 3. Similar results were found when reconstructing for domains with different background values (not shown). The method of selecting the optimum wavelength windows of the spectrum is entirely dependent on the number of chromophores considered, and thus the extinction coefficients of the available wavelengths. The optimum wavelengths are selected independent of the imaging volume, but since the problems of uniqueness and condition of the inverse problem have been considered objectively, the method is flexible to cope with a wide range of target and background values.
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% and 33% . With the windows method the error ranges from 18% to 30% . With six wavelengths, the range in error increases to between 24% (scattering amplitude) and 35% . 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 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 , the number of calculations is a function of and the memory of .23 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 to store the Hessian matrix, whereas the windows method (21 wavelengths) requires only . 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 , , 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 per iteration, while the full spectrum took . Thus the window method is computationally more efficient.
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.
This work has been sponsored by the National Cancer Institute through grants RO1CA78734 and PO1CA80139, as well as the Engineering and Physical Sciences Research Council, United Kingdom.