Spatial frequency domain imaging (SFDI) is a wide-field diffuse optical imaging modality that has attracted considerable interest in recent years. Typically, diffuse reflectance measurements of spatially modulated light are used to quantify the optical absorption and reduced scattering coefficients of tissue, and with these, chromophore concentrations are extracted. However, uncertainties in estimated absorption and reduced scattering coefficients are rarely reported, and we know of no method capable of providing these when look-up table (LUT) algorithms are used to recover the optical properties. We present a method to generate optical property uncertainty estimates from knowledge of diffuse reflectance measurement errors. By employing the Cramér-Rao bound, we can quickly and efficiently explore theoretical SFDI performance as a function of spatial frequencies and sample optical properties, allowing us to optimize spatial frequency selection for a given application. In practice, we can also obtain useful uncertainty estimates for optical properties recovered with a two-frequency LUT algorithm, as we demonstrate with tissue-simulating phantom and in vivo experiments. Finally, we illustrate how absorption coefficient uncertainties can be propagated forward to yield uncertainties for chromophore concentrations, which could significantly impact the interpretation of experimental results.
© 2018 Optical Society of America under the terms of the OSA Open Access Publishing Agreement
1 February 2018: A typographical correction was made to Ref. 19.
Spatial frequency domain imaging (SFDI) is a wide-field, noncontact diffuse optical imaging technique that has attracted considerable interest for a variety of applications. In clinical settings, it has been used to evaluate skin and breast lesions, while in mouse models, it has been employed to monitor the response to chemotherapy regimens, drug delivery to the brain, and the progression of Alzheimer’s disease [1–8]. Typically, diffuse reflectance measurements of spatially modulated light are used to quantify the optical absorption and reduced scattering coefficients (µa and ) of tissue, and with these, chromophore concentrations of interest are extracted (e.g., hemoglobin). However, uncertainties in estimated absorption and reduced scattering coefficients ( and ) are rarely reported, and we know of no method capable of providing such uncertainties when look-up table (LUT) inversion algorithms are used to recover the optical properties. Quantifying these uncertainties would have several important benefits. For example, they could be propagated forward to yield uncertainties in estimated chromophore concentrations, which could significantly affect the interpretation of experimental results. They could also be employed to help guide the selection of spatial frequencies used for SFDI measurements, given the requirements of the specific application.
In the diffuse optics community, theoretical performance prediction for system design or optimization is usually accomplished with Monte Carlo simulation [6, 9, 10] or with singular value decomposition (SVD) [11–14]. The former can be very time- and resource-intensive, and the results are valid only for the particular inversion algorithm simulated. Methods that rely on SVD avoid these difficulties, as they require only the sensitivity matrix of the forward problem. While evaluating the singular value spectrum or condition number of the sensitivity matrix is conceptually straightforward, this has a number of limitations. It ignores the information contained in the singular vectors of the sensitivity matrix, which can lead to errors in performance predictions [15, 16]. It also cannot incorporate the noise of the input data into the analysis. Finally, it can only provide qualitative estimates of performance (i.e., this system configuration is better than the others), not quantitative ones. In recent years, other alternatives to SVD have been proposed [17, 18], but to our knowledge, none of these are capable of providing quantitative estimates of performance.
In this work, we make three novel contributions. First, we show how knowledge of the precision of diffuse reflectance measurements from an SFDI instrument (i.e., diffuse reflectance uncertainty estimates, ) can be quickly and efficiently transformed to yield quantitative estimates of optical property uncertainties using the Cramér-Rao bound (CRB), a concept from the field of statistical signal processing. The CRB shares many of the advantages of SVD-based methods, but it can sometimes be difficult to compute. We leverage prior work in statistical signal processing for Gaussian data models and provide an expression for the CRB that applies to SFDI with an arbitrary number of spatial frequencies. We use this to evaluate trade-offs in theoretical system performance as a function of spatial frequency combinations and sample optical properties. Second, we demonstrate that a CRB-dervied metric can provide optical property uncertainty estimates when a two-frequency LUT inversion algorithm is used to recover µa and . We validate this with both simulated and phantom (experimental) data. We make use of this metric and LUT algorithm for our third contribution: an end-to-end example using in vivo data that culminates in uncertainties for chromophore concentration estimates.
The remainder of this paper is organized as follows. In Sec. 2, we describe our procedure for generating a diffuse reflectance error model by characterizing the precision of our SFDI instrument and present those results therein. Sections 3–5 contain the methods and results for the three contributions described above. We discuss and summarize our results in Sec. 6, and we present conclusions and a description of future work in Sec. 7. Custom Matlab (The MathWorks, Inc., Natick, MA) code developed for Sec. 3 can be downloaded from OSA’s figshare repository .
2. Diffuse reflectance error model
In order to estimate uncertainties in recovered optical property estimates, we need to know the uncertainties associated with the diffuse reflectance measurements. We assumed that our diffuse reflectance measurements were unbiased (i.e., systematic errors had been eliminated) and could therefore be modeled as
In Table 1, we summarize the relevant properties of our SFDI instrument and data acquisition and processing parameters. A detailed description of spatial frequency domain instrumentation and image formation can be found in [20, 21]. Briefly, the instrument projects spatially modulated 1-D sinusoidal patterns of light, of varying spectral wavelengths and spatial frequencies, onto the sample, and the resulting diffuse reflectance at each location is measured by a CCD camera. The data is demodulated with a conventional communications algorithm and converted into calibrated reflectance using a calibration phantom with known optical properties. In our case, the measured reflectance of the calibration phantom was compared to that predicted by a Monte Carlo simulation . This yields calibrated diffuse reflectance values over the surface of the sample (i.e., per pixel).
We characterized the precision of the instrument under these operating conditions by measuring 16 tissue-simulating phantoms with optical properties shown in Fig. 1(a). The optical properties are also listed in tabular form (Table 3) in Appendix A. The phantoms were homogeneous silicone blocks with Nigrosin and TiO2 added to control µa and , respectively, and they ranged in size from 13 × 20 × 3 cm3 to 7 × 11 × 3 cm3. Their optical properties spanned a range of values representative of typical applications in our lab, e.g., human breast and small animal (mouse) imaging. We made 10 measurements of each phantom at each wavelength and spatial frequency. The data was downsampled by a factor of two (2 × 2 binning) to increase the signal-to-noise ratio (SNR). After converting to calibrated diffuse reflectance, a circular region of interest (ROI) with a ~40 mm diameter was selected, and the mean and standard deviation of pixel values in this region for all 10 repeat measurements was computed (∼2 × 105 pixels per phantom). The phantoms were divided into two groups: 8 in the training set, and 8 in the test set, as shown in Fig. 1(a), such that the training set spanned the range of optical properties of interest. The diffuse reflectance error model was constructed by averaging the ratio of the standard deviation to the mean value (at each wavelength and spatial frequency) of the phantoms in the training set. The data in the test set was used to validate the model (see Sec. 4 for results).
The diffuse reflectance error model is shown in Fig. 1(b). We note that the noise as a function of spatial frequency was rather flat, while the mean diffuse reflectance decreased with increasing spatial frequency, as expected. This results in the increasing trends of for all four spectral wavelengths, i.e., the SNR worsens with increasing spatial frequency. The differences due to spectral wavelengths are likely caused by the particular properties of the LEDs in our instrument.
We conducted additional analyses to determine if the Gaussian noise assumption specified in Eq. (1) was appropriate for our diffuse reflectance measurements. In Fig. 2, we display histograms and quantile-quantile plots of calibrated diffuse reflectance from the ROIs of three different phantoms. The data in Figs. 2(a) and 2(d) shows excellent agreement with the Gaussian noise assumption. However, some data does deviate from the Gaussian model, exhibiting some kurtosis (Figs. 2(b) and 2(e)) and skewness (Figs. 2(c) and 2(f)). Nevertheless, for the majority of the data, excess kurtosis and skewness are within ±0.5. (Both values are 0 for the normal distribution.) So overall it appears that our diffuse reflectance data is well described by the Gaussian noise model.
3. Optical property uncertainty estimates using the Cramér-Rao bound
To transform diffuse reflectance precision estimates into precision estimates for the absorption and reduced scattering coefficients ( and ), we employ the CRB . The CRB considered here is a lower bound that defines the best theoretical precision (i.e., lowest variance) of any unbiased estimator for a given data model. We distinguish here between the concept of a “data model” (e.g., how optical properties determine measured Rd) and an “estimator” or “algorithm” (e.g., how measured Rd can be used to estimate optical properties). Loosely speaking, the data model defines the forward problem, while the estimator or algorithm solves the inverse problem. The CRB is often used in the statistical signal processing community, especially in the sonar and radar signal processing communities, to perform feasibility studies and system design. The basic idea is that the curvature of the log-likelihood function of the data model determines the precision with which an unknown (e.g., µa, ) may be estimated. The more sharply peaked the log-likelihood function is in the neighborhood of the unknown, the more precisely the value of the unknown may be determined. The CRB provides a measure of this curvature. For a vector of unknowns, it is computed by first calculating the Fisher information matrix, which requires second derivatives and first derivative cross products of the log-likelihood function with respect to each unknown, and then taking its inverse.
Because it uses information only about the data model, the CRB is independent of the particular algorithm employed to estimate the unknown(s). It is therefore ideally suited to explore trade-offs in performance as a function of the data, which, for our purposes, includes data with different spatial frequency combinations and sample optical properties. We note that the CRB has an additional advantage over SVD-based methods: it includes the impact of the input data noise.
3.1.1. Cramér-Rao bound for spatial frequency domain imaging
Let us organize our data, which consists of calibrated diffuse reflectance measurements for a given pixel (described by Eq. (1)) at K spatial frequencies, where K ≥ 2, into a column vector:23]. We have suppressed the dependency of on λ in Eq. (5) for clarity. As we seek to estimate both µa and , i and j will take on values of 1 and 2, so F(θ) will be a 2 × 2 matrix.
We computed the vectors and matrices of derivatives,Eqs. (6a)–(6d) can be computed by numerically differentiating each slice of the cube in the x and y directions, respectively, and evaluating the derivative at θ:
The CRB for θ is obtained by taking the inverse of F(θ):
We can define a second quantity that employs only the first term of the Fisher information matrix:Eq. (5) as Eqs. (4), Eq. (6), and (7) into Eq. (5). It reveals that when , the difference between the rCRB (computed with the first term) and CRB (computed with both terms) will be very small, and hence, for most cases of practical interest to us, the CRB is well approximated by the simpler expression of the rCRB. We make use of the rCRB in Secs. 4 and 5. Custom Matlab code that computes the CRB and rCRB is available in Code 1,  >.
3.1.2. Diffuse reflectance cube
We computed the diffuse reflectance cube in two different ways. Initially, we used a Monte Carlo simulation of diffuse reflectance for a semi-infinite homogeneous medium that was subsequently scaled for different values of µa and , as described in . We then took the Hankel transform to obtain results in the spatial frequency domain. However, we noticed some artifacts when computing the derivatives (Eqs. (7a)–(7d)) at spatial frequencies higher than 0.4 mm−1. These are likely caused by the lower SNR at higher spatial frequencies coupled with the fact that differentiation amplifies noise. So we created a second cube employing the solution for the diffusion equation , and we used it for the results presented in Figs. 3, 4, and 6, which involved frequencies ≥ 0.4mm−1.
For SFDI, the diffusion equation approximation is considered valid when fx ≪ µtr, where . For two of the four pairs of optical properties considered in these figures, this approximation should start to break down at fx ~ 0.2 mm−1 . However, the differences between CRB results generated with the diffusion-based cube and the Monte Carlo-based one were less than 3% for all spatial frequencies below 0.4 mm−1. This led us to believe that the diffusion-based cube was capable of adequately capturing the trends as a function of spatial frequency that we were looking to explore.
In addition to providing us with a convenient way of evaluating Eqs. (7a)–(7d), diffuse reflectance cubes can serve as the core of an LUT algorithm that relates values of Rd at two or more spatial frequencies to unique pairs of µa and . We implemented our LUT algorithm with Matlab’s “griddata” function, which accepts as input a cube and measured Rd values and returns the best fit optical properties, where “best fit” is determined by linear interpolation of the values in the cube. Although it may seem that we are linking the CRB computation—said to be algorithm-independent earlier—with our optical property recovery algorithm by using the same cube in both, notice that the CRB is agnostic as to how the information in the cube is used by the algorithm. For example, an LUT algorithm could employ linear, nearest neighbor, or cubic interpolation, or perform polynomial fitting . The cube could also be used by a maximum likelihood estimator to create a chi-squared surface that could be minimized in a number of different ways. All of these particulars can significantly affect algorithm performance, but none of them appear in the CRB calculation. Indeed, the cube is just a numerical model that captures how µa, , and fx determine Rd, and it is not surprising that optical property recovery algorithms would make use of that information somehow.
We used the CRB (Eq. (8)) to explore performance as a function of spatial frequency combinations and sample optical properties. We show results using the diffuse reflectance error model for 659 nm, but the other three wavelengths yielded very similar patterns. In Fig. 3, we address the question of which spatial frequency to pair with DC (0 mm−1) for the determination of four sets of optical properties relevant to our mouse studies . We plot the uncertainties in recovered µa and as a function of the spatial frequency combination used. The sweet spot near fx = [0, 0.1] mm−1 is likely caused by the confluence of two opposing trends: (1) performance improves as the gap between the two frequencies increases (as discussed in ); and (2) performance worsens as the SNR decreases with increasing spatial frequency (see Fig. 1(b)). An interesting pattern that emerged is that and are consistently lower for the lower value of , irrespective of absorption.
In Fig. 4, we explore the potential advantage of using more than two frequencies to recover the sample optical properties, as well as the impact of not including DC. With respect to absorption uncertainty (Figs. 4(a) and 4(b)), we see that a well chosen spatial frequency pair (e.g., [0, 0.1] mm−1) performs almost as well as 3 or even 13 spatial frequencies. The diminishing returns of adding more spatial frequencies was noted in , albeit under rather different circumstances. We also see that not including DC significantly degrades the precision of the µa estimate, an effect that is especially pronounced for the lower value of µa. Since a low spatial frequency like DC is known to propagate more deeply into tissue , it is perhaps not surprising that it is needed to adequately sample the rarer absorption events that correspond to very low absorption coefficients. The impact on scattering uncertainty, however, is much less severe, as none of the errors exceed 9% (Fig. 4(c)). We see the same pattern here as in Fig. 3(b): is consistently lower for the lower value of , irrespective of absorption.
Finally, in Fig. 5 we consider the impact of sample optical properties on the precision of µa and estimates obtained using the spatial frequency pair [0, 0.1] mm−1. Here we expand the range of optical properties considered well beyond that of Figs. 3 and 4. We see that while performance is reasonable for optical properties in the neighborhood of our previously considered set of four, higher variances are obtained for increasing values of µa and (Figs. 5(a) and 5(b)). The variances are not only larger in magnitude but also highly correlated (Fig. 5(c)). Thus the expected range of optical properties for a given application, along with desired precision, should inform the choice of spatial frequency combinations. These CRB “maps” are a quick and efficient way to explore this trade space. The maps in Fig. 5, which contain approximately 6 million combinations of µa and , were generated in ~5 min on a desktop (Intel Core i7-6700 Processor, 16 GB RAM) using non-optimized Matlab code.
4. Optical property uncertainty estimates for a two-frequency look-up table inversion algorithm
While the CRB has a number of advantages for feasibility studies and system design, there is no guarantee that the bound can be attained by a realizable estimator . Under some circumstances, it is possible to prove that the CRB will be attained, at least asymptotically, as with maximum likelihood estimators . But in general, there is no way to know if a particular algorithm’s performance can be described by this bound. In this section, we establish that the rCRB can be used to predict the performance of a two-frequency LUT inversion algorithm. This algorithm, which uses the Matlab function “griddata,” was described in Sec. 3.1.2. We characterized its performance with simulated diffuse reflectance data first, and then evaluated it with phantom data.
We created simulated diffuse reflectance data according to Eq. (1). The diffusion-based cube (Sec. 3.1.2) was used to define Rd, while the diffuse reflectance error model at 659 nm (Sec. 2) provided the components of . We generated 104 samples for each of the optical property combinations shown in Fig. 3 at the 13 spatial frequencies listed in Table 1. These samples were fed to the LUT algorithm, using DC paired with one other spatial frequency at a time, and values for the mean and standard deviation of the recovered optical properties were computed. The diffusion-based cube was used to compute the CRB and rCRB (as described in Sec. 3.1).
We also tested the algorithm with calibrated diffuse reflectance data from the tissue-simulating phantoms described in Sec. 2. Since these are homogeneous phantoms, the optical properties should be the same over the surface of the sample, so estimates of optical property precision can be obtained by computing the standard deviation of recovered µa and values. We randomly selected 104 samples from the ROI of each phantom and recovered the optical properties using two spatial frequency combinations: [0, 0.1] and [0.05, 0.1] mm−1. We chose these two pairs based on earlier results that revealed them to be good ([0, 0.1]) and less ideal ([0.05, 0.1]) combinations for optical property recovery in the range spanned by our phantoms (e.g., see Fig. 4(a)). We used the Monte Carlo-based cube in the LUT algorithm and in the CRB and rCRB computations, as all frequencies considered are well below 4 mm−1.
While the CRB and rCRB computations yielded very similar values, we found that the rCRB was, on average, a better match for LUT algorithm performance than the CRB (results not shown). Hence, this became our method of choice when the task was to estimate LUT algorithm uncertainties, and we used this quantity in the performance comparisons described in this section. As the rCRB was derived from the CRB, which assumed an unbiased estimator, we first checked the LUT algorithm results for bias by comparing the means of the recovered optical properties to the true values used in the simulations. For the lowest spatial frequency pair ([0, 0.025] mm−1), we found evidence of moderate bias, but for all other pairs, the average absolute deviation from the true values was quite small: 0.3% for µa and 0.2% for . This pattern recurs in Figs. 6(a) and 6(b), where we see excellent agreement between the rCRBs (lines) and the standard deviations of the recovered optical properties (symbols) for all but the lowest spatial frequency pair. Excluding this pair, the difference is ~0.1%, on average, for both and .
We note that for spatial frequency combinations that are close together and have large CRB values (e.g., [0, 0.025] mm−1), the LUT algorithm sometimes fails to produce a solution. This is because the inherent noise results in a value of diffuse reflectance that lies outside of the LUT. These cases bias the optical property estimates and produce values of and that do not match the rCRBs or CRBs.
In Table 2, we summarize the results obtained with phantom data. We have excluded from our statistics 6 cases (out of 96 total) for which the LUT algorithm failed to work reliably. Recall that the phantom data was divided into two groups: the training set, which was used to create the diffuse reflectance error model, and the test set. We see no significant difference in performance between these two groups, which suggests that an error model can be successfully created by measuring a limited number of phantoms that span the optical properties of interest. In comparing the results at 659 nm for the spatial frequency pairs [0, 0.1] and [0.05, 0.1] mm−1, we see that the mean absolute difference between LUT algorithm performance and the rCRB is slightly larger for the second pair, especially for . However, at roughly 2%, this value is still quite low and indicates good performance.
For a “good” spatial frequency pair like [0, 0.1] mm−1, agreement between the rCRB and LUT algorithm performance is generally excellent. Averaged over all four spectral wavelengths, the mean absolute difference is ~ 1% for both and . This is in spite of the fact that some phantom data does not satisfy the Gaussian noise assumption very well, as discussed in Sec. 2. In Fig. 7, we display the individual data points summarized in the second line of the table. It appears that the rCRB is able to provide useful uncertainty estimates for our LUT algorithm.
5. Uncertainties for chromophore concentration estimates
We have shown how optical property uncertainties can be obtained from diffuse reflectance uncertainties. Here we take it one step further and demonstrate, with an arterial cuff occlusion experiment, how uncertainties for chromophore concentration estimates can be obtained from absorption coefficient uncertainties.
5.1.1. Expression for chromophore concentration uncertainties
In near-infrared spectroscopy, absorption coefficients are usually related to chromophore concentrations via a linear model:23]: Eq. (14)) is optimal . This means that is equal to the CRB for this problem and could be used to optimize the choice of spectral wavelengths, as was done for spatial frequencies in Sec. 3.
5.1.2. Arterial cuff occlusion experiment
We performed an arterial cuff occlusion test on a 30-year-old healthy male. This experiment was designed to generate a well-known “signal” that could be used to help us evaluate the impact of chromophore concentration uncertainties along with some of the other claims made in this paper. The cuff was applied to the upper arm, and we measured the inner forearm with our SFDI instrument. We took measurements at three spatial frequencies (0, 0.05, 0.1 mm−1) and four spectral wavelengths (659, 691, 731, 851 nm). Measurements were made every ~8.6 s for a total of 7 min. During the first 2 min, no pressure was applied to the cuff. At t = 2 min, the pressure was rapidly increased to ~ 200 mmHg. At t = 4 min, we released the pressure and continued measuring for another 3 min.
The data was converted to calibrated reflectance as described in Sec. 2. As with the phantom data, we performed 2 × 2 binning to increase the SNR. The optical properties were recovered with the two-frequency LUT algorithm (Sec. 4) using [0, 0.1] mm−1 and [0.05, 0.1] mm−1, and the corresponding uncertainties were estimated with the rCRB (Sec. 3). The Monte Carlo-based cube was used in the LUT algorithm and rCRB calculations. To further reduce the noise in our measurements, we averaged the recovered optical properties over a small ROI (diameter ~4 mm). We chose the ROI to be large enough to give us significant variance reduction but small enough that the optical properties within it could be considered uniform. This gave us one pair of µa and for each time point (50 total during the 7 min experiment).
The absorption coefficients at the four spectral wavelengths were converted into chromophore concentrations with Eq. (14), and the chromophore uncertainties were computed with Eq. (15). We focused on measuring changes (relative to the first measurement) in oxy- (O2Hb) and deoxyhemoglobin (HHb) concentrations, so we solved only for these two chromophores. The chromophore extinction coefficients came from . Because we took the average of µa over an ROI, (assumed constant over the ROI) should be divided by the number of degrees of freedom to properly reflect the uncertainty in the averaged value. If the pixels were independent, the number of degrees of freedom would be equal to the number of pixels in the ROI. However, pixels in an SFDI image are typically not statistically independent, as the diffuse nature of light propagation results in spatial correlations. We can get an idea of the extent of the spatial correlations by taking the Fourier transform of the modulation transfer function in the spatial frequency domain (i.e., Rd as a function of spatial frequency), which gives us the point spread function (PSF) of the medium in space. We did this, assuming average values for µa and of 0.01 and 1.86 mm−1, respectively, for all wavelengths over the course of the experiment, and arrived at a full-width half-maximum value of ~1 mm for the PSF. Given the long tails of the PSF, we decided to use 2 mm as our estimate of the PSF diameter. The ratio of the ROI to PSF area is therefore , so we used 4 as the effective degrees of freedom in our computation of for the chromophore concentration uncertainties. We note that this is only a rough estimate, and we are currently working to develop methods that can better quantify this value.
In Fig. 8, we display the results of our cuff occlusion experiment with data at two pairs of spatial frequencies: [0, 0.1] mm−1 and [0.05, 0.1] mm−1. The data in Fig. 8(a) unambiguously shows the expected increase in HHb during the occlusion, coupled with the smaller rate of decrease in O2Hb. The latter may be caused by an incomplete occlusion of the arterial blood supply. After releasing the cuff pressure, HHb decreases rapidly and O2Hb increases. Both chromophore concentrations look like they are returning to baseline when the experiment ends 3 min later. On average, the uncertainty in O2Hb concentration was ±1.8 µM, while the uncertainty in HHb concentration was ±0.73 µM. The behavior of O2Hb and HHb concentrations during the experiment is consistent with previously reported results [14, 27]. Interestingly, the average uncertainty of the less noisy signal, HHb, is roughly a factor of two smaller than that of O2Hb. Both uncertainties make it clear that we are seeing significant changes in chromophore concentrations as a result of the occlusion.
This case is much harder to make with the data in Fig. 8(b), which is obviously a lot noisier. It is difficult to confirm the expected trend in HHb or to see any pattern at all in O2Hb. And the much larger uncertainty of the latter suggests that one cannot make any meaningful claims about O2Hb concentration with this data. This stark difference in performance with the two spatial frequency pairs is not very surprising, given the CRB results of Fig. 4. Here we saw that values of for [0, 0.1] mm−1 ranged from 5 to 11%, while for [0.05, 0.1] mm−1, they ranged from 19 to 89%. This points to the importance of choosing the right spatial frequency combination for a given application and confirms the ability of the CRB to provide useful guidance for in vivo measurements.
6. Summary and discussion
In this work, we employed the CRB to estimate uncertainties in recovered optical properties using knowledge of diffuse reflectance measurement errors. The CRB is ideally suited to explore how theoretical SFDI performance varies as a function of sample optical properties and the spatial frequencies used to recover them, since it is algorithm-independent and can take into account data noise. We made two choices that allowed us to compute the CRB for our application: (1) diffuse reflectance data was modeled with a multivariate normal distribution; and (2) derivatives of Rd were computed by numerical differentiation. We validated the first choice, to the extent possible, with data from tissue-simulating phantoms. With respect to the second choice, evaluating the derivatives numerically proved not only convenient, but it gives one the flexibility to use Monte Carlo simulations of Rd, which may lead to more accurate results in some cases.
Although popular in other communities, the CRB has not received much attention in the diffuse optics community, in part likely due to the difficulty of applying it to ill-posed inverse problems (as discussed in ). But for many SFDI applications, where the goal is not tomography, the recovery of optical properties with data at two or more spatial frequencies is not an ill-posed problem. We can apply the same approach to other diffuse optics problems that are not ill-posed, e.g., the determination of optical properties from frequency domain measurements of amplitude and phase. If we can make the same two choices as above for the relevant data (e.g., amplitude and phase), computing the CRB for those applications could be accomplished with rather straightforward modifications of the code that we have made available.
We used the CRB to explore theoretical SFDI performance. One rather surprising result was that for optical properties typical of human breast or mouse imaging, there is not much benefit in using multiple spatial frequencies versus a well chosen pair. Selecting DC (0 mm−1) as one member of that pair led to better µa estimation performance. We also saw that for a given set of spatial frequencies, and can vary significantly depending on the sample’s optical properties. We note that these results rely on the diffuse reflectance error model of Sec. 2. With a different error model (e.g., different instrument or operating conditions), these conclusions may not hold, so characterizing is an indispensable prerequisite for this analysis.
From the CRB derivation, we arrived at a simpler quantity (the rCRB) capable of providing uncertainties when optical properties are recovered with a two-frequency LUT inversion algorithm. We characterized the performance of this algorithm with simulated diffuse reflectance data first, and then evaluated it with (experimental) phantom data, which at times deviated from the assumed data model. In both cases, the rCRB was able to successfully predict the uncertainties in estimated optical properties. With experimental data, the difference between predicted and measured and was ~1%, on average. This gives us a quick and efficient way to generate uncertainties for an LUT algorithm. It opens the door to the possibility of pre-computing an LUT of optical property uncertainties—in essence, the maps of Fig. 5—which could be a big asset for applications requiring near real-time performance. It also reveals that in the case of two spatial frequencies, a simple LUT algorithm is capable of achieving near-optimal performance.
Finally, we provided an end-to-end example, using in vivo data, that illustrated how absorption coefficient uncertainties could be propagated forward to yield uncertainties for chromophore concentration estimates. Since a linear model relates chromophore concentrations to absorption coefficients, obtaining uncertainties for the former is straightforward once the covariance matrix of the absorption coefficients is known. Our example reveals that chromophore concentration uncertainties can have an important impact on the interpretation of experimental results. While our method provides uncertainties at each pixel, it is not uncommon to average over an ROI when analyzing in vivo data, as we have done. To compute uncertainties when spatial averaging over an ROI is performed, we must develop a robust way to estimate the number of degrees of freedom, as pixels in an SFDI image are usually not statistically independent. This is an area of ongoing research in our group.
7. Conclusions and future work
We have developed a method to generate optical property uncertainty estimates from knowledge of diffuse reflectance measurement errors. By relying on the CRB, we are able to do this quickly and efficiently with non-optimized Matlab code that takes minutes to run on a desktop computer. Our method is ideally suited to characterize theoretical SFDI performance as a function of variables like spatial frequency combinations and sample optical properties. In practice, it can also provide useful uncertainty estimates for optical properties recovered with a two-frequency LUT inversion algorithm, as was demonstrated with phantom and in vivo data. We also showed how absorption coefficient uncertainty estimates could be propagated forward to yield uncertainties for chromophore concentration estimates. To fully realize this potential, however, a robust way of estimating the effective degrees of freedom when spatial averaging is performed over an ROI is needed. This work is currently in progress.
One important aspect of SFDI performance that we did not explore here is the impact of penetration depth. As mentioned in Sec. 3.2, it is known that penetration depth decreases with increasing spatial frequency. However, SFDI data processing assumes that the frequencies employed interrogate volumes with the same optical properties. This is true when measuring a homogeneous phantom, but it is likely not the case for heterogeneous biological tissue. To minimize the potential sampling mismatch, one might be tempted to choose spatial frequencies that are close together. But our CRB analysis shows that it is beneficial to have some separation between spatial frequencies, e.g., [0, 0.1] mm−1 performs better than [0, 0.025] mm−1. Developing a framework that can reconcile both of these concerns would be a valuable step forward.
This work was funded by a grant from the Dept. of Defense (Award No. W81XWH-15-1-0070).
We thank Yanyu Zhao, Rachita Chaudhury, Alyssa Torjesen, and Claire Fu for making the phantoms measured in this study.
The authors declare that there are no conflicts of interest related to this article.
References and links
1. A. Yafi, T. S. Vetter, T. Scholz, S. Patel, R. B. Saager, D. J. Cuccia, G. R. Evans, and A. J. Durkin, “Postoperative quantitative assessment of reconstructive tissue status in a cutaneous flap model using spatial frequency domain imaging,” Plast. Reconstr. Surg. 127(1), 117–130 (2011). [CrossRef] [PubMed]
2. D. J. Rohrbach, N. C. Zeitouni, D. Muffoletto, R. Saager, B. J. Tromberg, and U. Sunar, “Characterization of nonmelanoma skin cancer for light therapy using spatial frequency domain imaging,” Biomed. Opt. Express 6(5), 1761–1766 (2015). [CrossRef] [PubMed]
3. R. B. Saager, A. N. Dang, S. S. Huang, K. M. Kelly, and A. J. Durkin, “Portable (handheld) clinical device for quantitative spectroscopy of skin, utilizing spatial frequency domain reflectance techniques,” Rev. Sci. Instrum. 88(9), 094302 (2017). [CrossRef] [PubMed]
4. S. Gioux, A. Mazhar, B. T. Lee, S. J. Lin, A. M. Tobias, D. J. Cuccia, A. Stockdale, R. Oketokoun, Y. Ashitate, E. Kelly, M. Weinmann, N. J. Durr, L. A. Moffitt, A. J. Durkin, B. J. Tromberg, and J. V. Frangioni, “First-inhuman pilot study of a spatial frequency domain oxygenation imaging system,” J. Biomed. Opt. 16(8), 086015 (2011). [CrossRef] [PubMed]
5. C. M. Robbins, G. Raghavan, J. F. Antaki, and J. M. Kainerstorfer, “Feasibility of spatial frequency-domain imaging for monitoring palpable breast lesions,” J. Biomed. Opt. 22(12), 121605 (2017). [CrossRef]
6. S. Tabassum, Y. Zhao, R. Istfan, J. Wu, D. J. Waxman, and D. Roblyer, “Feasibility of spatial frequency domain imaging (SFDI) for optically characterizing a preclinical oncology model,” Biomed. Opt. Express 7(10), 4154–4170 (2016). [CrossRef] [PubMed]
7. R. P. Singh-Moon, D. M. Roblyer, I. J. Bigio, and S. Joshi, “Spatial mapping of drug delivery to brain tissue using hyperspectral spatial frequency-domain imaging,” J. Biomed. Opt. 19(9), 096003 (2014). [CrossRef]
8. A. J. Lin, M. A. Koike, K. N. Green, J. G. Kim, A. Mazhar, T. B. Rice, F. M. LaFerla, and B. J. Tromberg, “Spatial frequency domain imaging of intrinsic optical property contrast in a mouse model of Alzheimer’s disease,” Ann. Biomed. Eng. 39(4), 1349–1357 (2011). [CrossRef] [PubMed]
9. 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. Imag. 21(7), 755–763 (2002). [CrossRef]
10. A. B. Milstein, J. J. Stott, S. Oh, D. A. Boas, R. P. Millane, C. A. Bouman, and K. J. Webb, “Fluorescence optical diffusion tomography using multiple-frequency data,” J. Opt. Soc. Am. A 21(6), 1035–1049 (2004). [CrossRef]
11. J. P. Culver, V. Ntziachristos, M. J. Holboke, and A. G. Yodh, “Optimization of optode arrangements for diffuse optical tomography: A singular value analysis,” Opt. Lett. 26(10), 701–703 (2001). [CrossRef]
12. F. Leblond, H. Dehghani, D. Kepshire, and B.W. Pogue, “Early-photon fluorescence tomography: spatial resolution improvements and noise stability considerations,” J. Opt. Soc. Am. A 26(6), 1444–1457 (2009). [CrossRef]
13. X. Zhang and C. Badea, “Effects of sampling strategy on image quality in noncontact panoramic fluorescence diffuse optical tomography for small animal imaging,” Opt. Express 17(7), 5125–5138 (2009). [CrossRef] [PubMed]
14. A. Mazhar, S. Dell, D. J. Cuccia, S. Gioux, A. J. Durkin, J. V. Frangioni, and B. J. Tromberg, “Wavelength optimization for rapid chromophore mapping using spatial frequency domain imaging,” J. Biomed Opt. 15(6), 061716 (2010). [CrossRef]
15. A. J. Chaudhari, F. Darvas, J. R. Bading, R. A. Moats, P. S. Conti, D. J. Smith, S. R. Cherry, and R. M. Leahy, “Hyperspectral and multispectral bioluminescence optical tomography for small animal imaging,” Phys. Med. Biol. 50(23), 5421–5441 (2005). [CrossRef] [PubMed]
16. F. Leblond, K. M. Tichauer, and B. W. Pogue, “Singular value decomposition metrics show limitations of detector design in diffuse fluorescence tomography,” Biomed. Opt. Express 1(5), 1514–1531 (2010). [CrossRef]
17. D. Karkala and P. K. Yalavarthy, “Data-resolution based optimization of the data-collection strategy for near infrared diffuse optical tomography,” Med. Phys. 39(8), 4715–4725 (2012). [CrossRef] [PubMed]
18. R. W. Holt, F. L. Leblond, and B. W. Pogue, “Methodology to optimize detector geometry in fluorescence tomography of tissue using the minimized curvature of the summed diffuse sensitivity projections,” J. Opt. Soc. Am. A 30(8), 1613–1619 (2013). [CrossRef]
19. V. Pera, “Cramér-Rao bounds for spatial frequency domain imaging,” figshare (2017) [uploaded 20 Dec 2017], https://doi.org/10.6084/m9.figshare.5715100.
20. D. J. Cuccia, F. Bevilacqua, A. J. Durkin, and B. J. Tromberg, “Modulated imaging: quantitative analysis and tomography of turbid media in the spatial-frequency domain,” Opt. Lett. 30(11), 1354–1356 (2005). [CrossRef] [PubMed]
21. D. J. Cuccia, F. Bevilacqua, A. J. Durkin, F. R. Ayers, and B. J. Tromberg, “Quantitation and mapping of tissue optical properties using modulated imaging,” J. Biomed. Opt. 14(2), 024012 (2009). [CrossRef] [PubMed]
22. M. Martinelli, A. Gardner, D. Cuccia, C. Hayakawa, J. Spanier, and V. Venugopalan, “Analysis of single Monte Carlo methods for prediction of reflectance from turbid media,” Opt. Express 19(20), 19627–19642 (2011). [CrossRef]
23. S. M. Kay, Fundamentals of Statistical Signal Processing: Estimation Theory (Prentice-Hall, 1993).
25. T. A. Erickson, A. Mazhar, D. Cuccia, A. J. Durkin, and J. W. Tunnell, “Lookup-table method for imaging optical properties with structured illumination beyond the diffusion theory regime,” J. Biomed. Opt. 15(3), 036013 (2010). [CrossRef] [PubMed]
26. S. Prahl, “Tabulated molar extinction coefficient for hemoglobin in water,” (Oregon Medical Laser Center, 1998), http://omlc.ogi.edu/spectra/hemoglobin/summary.html.
27. F. Teng, T. Cormier, A. Sauer-Budge, R. Chaudhury, V. Pera, D. Chargin, S. Brookfield, N. Y. Ko, and D. Roblyer, “Wearable near-infrared optical probe for continuous monitoring during breast cancer neoadjuvant chemotherapy infusions,” J. Biomed. Opt. 22(1), 014001 (2017). [CrossRef]