Rapid quantitative imaging of tissue optical properties, namely absorption (${\mu }_{\mathrm{a}}$) and reduced scattering (${\mu }_{\mathrm{s}}^{\prime }$), has long been a challenge in the field of biomedical optics. The recent introduction of spatial frequency domain imaging (SFDI) provides a transformative approach capable of measuring optical properties over a large field-of-view.¹ In principle, SFDI relies on the analysis of the tissue response to structured illumination (patterns of light) in the spatial frequency domain in a multipixel manner over an entire field-of-view at once.^2,3 More specifically, the spatial-frequency-dependent response, called the modulation transfer function, is calibrated using a phantom with known optical properties to determine the tissue diffuse reflectance ($R_d$), which is then used to extract the optical properties using a light propagation model.

While the method is rapid, SFDI has until recently involved the measurement of several images to form maps of optical properties, typically a total of six images (two spatial frequencies at three phases), although a method employing only a single spatial frequency and three phases has also been described. Recent developments in acquisition methods reduced the number of images necessary to extract optical properties from 6 to 2,⁴ and even from 6 to a single image with the method called single snapshot of optical properties (SSOP).⁵ In their most advanced implementations, such methods enable profile-corrected measurements of optical properties in a single snapshot, in turn facilitating true real-time acquisition of optical properties.⁶

Unfortunately, while acquisition methods can be performed in real time, the processing is still commonly achieved postacquisition. The most rapid processing method employs a lookup table (LUT) approach, by which solutions for diffuse reflectance are generated from a light propagation model (diffusion, Monte Carlo) or empirically with various spatial frequencies and optical properties.^3,7,8 Such a method allows one to directly link a calibrated diffuse reflectance measurement at known spatial frequencies to a unique solution of optical properties (${\mu }_{\mathrm{a}}$ and ${\mu }_{\mathrm{s}}^{\prime }$). However, because it necessitates interpolations within the LUT for each pixel in the image, this method is still time-consuming (seconds), preventing the use of SFDI for true real-time imaging of optical properties.

In this letter, we present forms of LUTs allowing rapid extraction of optical properties from the measurement of calibrated diffuse reflectance. We propose two solutions that do not necessitate interpolations and thus dramatically reduce computation time. The first method consists of a hyper-dense linearized LUT, and the second invokes an analytical representation of the LUT. These methods are described, implemented, and compared to the standard LUT method in terms of precision, accuracy, and computation time. Combined with real-time acquisition, this letter facilitates real-time quantitative optical imaging of tissue properties.

The two methods described in this letter are based on a previously developed LUT.³ This standard LUT method uses “white” Monte-Carlo simulation (WMC, with zero absorbance) to model the spatially resolved impulse response, i.e., the steady-state diffuse reflectance $R_d(\rho )$, of a collimated point-source illumination for a given set of ${\mu }_{\mathrm{a}}$, ${\mu }_{\mathrm{s}}^{\prime }$, $n$, and $g$, which is then Fourier transformed to determine the spatial frequency response of the diffuse reflectance.^9,10 To generate the standard LUT, WMC was used to simulate ${10}^7$ photons into a homogenous medium with index of refraction $n=1.4$ and anisotropy factor $g=0.9$. The model used a detector with a numerical aperture of 0.22 and the radial bins were sampled in increments of $\mathrm{\Delta }\rho =0.09\mathrm{mm}$, allowing a maximum spatial frequency of over $5{\mathrm{mm}}^{-1}$. This process is repeated over several optical properties, and stored in a table that associates the diffuse reflectance at two spatial frequencies with a unique pair of optical properties. Recovering optical properties, therefore, involves measuring the diffuse reflectance of the sample and searching the table to find the corresponding optical property values (hence the name “lookup table”). With this method, it is important to understand that the optical properties are linearly sampled, so we refer to this method as the linear OP LUT. Using this method, the measured $R_d$ values are used to search through the nonlinear mapping of $R_d$ against the tissue properties in the linear OP LUT table, and then ${\mu }_{\mathrm{a}}$ and ${\mu }_{\mathrm{s}}^{\prime }$ are acquired using cubic spline interpolation. For our linear OP LUT, the resolution of the sampled optical properties was 0.001 and $0.01{\mathrm{mm}}^{-1}$ with a range of $[0,0.3]{\mathrm{mm}}^{-1}$ and $[0.3,3]{\mathrm{mm}}^{-1}$ for ${\mu }_{\mathrm{a}}$ and ${\mu }_{\mathrm{s}}^{\prime }$, respectively, with a table size of 81,571 points (see linear OP LUT in Fig. 1).

To avoid the time-consuming searching and interpolation steps required for the standard linear OP LUT, an LUT linear with respect to $R_d$ (linear $R_d$ LUT) was made by interpolating linearly sampled $R_d$ values from 0 to 1 by 500 points within a high-resolution linear OP LUT (0.0001 and $0.001{\mathrm{mm}}^{-1}$ resolution with a range of $[0,0.3]{\mathrm{mm}}^{-1}$ and $[0.3,3]{\mathrm{mm}}^{-1}$ for ${\mu }_{\mathrm{a}}$ and ${\mu }_{\mathrm{s}}^{\prime }$, respectively, with a table size of 8,105,701 points). The resolution of the $R_d$ linear sampling is chosen so that measured $R_d$ values can be found directly in the table by rounding their value to the nearest increment of 0.02, without interpolation (see linear $R_d$ LUT in Fig. 1). Given a sample’s $R_d$ values, the LUT’sorresponding indices can be calculated by parameterizing the linear $R_d$ grid and rounding to the nearest vertex to acquire ${\mu }_{\mathrm{a}}$ and ${\mu }_{\mathrm{s}}^{\prime }$.

Using the linear $R_d$ LUT, each optical property table was fit to a two-dimensional (2-D) function using a nonlinear least squares solver (“fit” using MATLAB)

where $a_1$ through $d_4$ are optimized parameters. Using these functions, optical properties can be directly deduced from the diffuse reflectance at low spatial frequency (noted DC here) and the diffuse reflectance at high spatial frequency (noted AC here). This is referred to as the 2-D fit LUT (see 2-D fit LUT in Fig. 1).

Figure 1 summarizes the LUT formations and workflows. The white arrows describe the formation process from Monte-Carlo simulation to the linear OP LUT (arrow 1), to the linear $R_d$ LUT (arrow 2), and finally to the 2-D fit LUT (arrow 3). The black arrows indicate the three options for workflow when using the LUTs to extract optical properties. Diffuse reflectance measurements at two spatial frequencies (low: DC, and high: AC) are used as inputs to the LUTs to extract the optical properties.

The custom imaging system and associated processing have been described extensively and results have been published, including a clinical trial.^11,12 Briefly, the system utilizes a digital micromirror device (DMD) to project patterns of 670-nm laser illumination which is cross polarized with the collection optics to minimize specular reflections in images collected on a CCD camera.

An array of tissue-like silicone phantoms was made using India ink (Blick Art Materials, Boston, Massachusetts) and titanium dioxide (Atlantic Equipment Engineers, Bergenfield, New Jersey) for absorption and scattering, respectively. The array was made to span the optical properties ${\mu }_{\mathrm{a}}=[0.01,0.1]{\mathrm{mm}}^{-1}$ and ${\mu }_{\mathrm{s}}^{\prime }=[0.5,2]{\mathrm{mm}}^{-1}$. The array was imaged using SFDI and the inversion from $R_d$ to optical properties was then evaluated on all three LUT methods and compared to a Monte Carlo (MC) least-squares solver over square ROIs of $15\times 15\text{pixels}$.⁹

Finally, a Yorkshire pig was used as an in vivo model to evaluate each LUT method during realistic surgical conditions (Beth Israel Deaconess Medical Center Institutional Animal Care and Use Committee approved institutional protocol #034-2013). Abdominal skin flaps were elevated on a pair of perforator vessels, the venous pedicle was occluded, and released to restore tissue perfusion.

All computational processing and analyses were done on a 64-bit Dell Optiplex 9020 (Dell, Round Rock Texas) with an Intel Core i7-4770 CPU (Intel, Santa Clara, California). $R_d$ inversions were repeated 100 times for each LUT type and averaged for comparison.

Figure 2 shows the results for the phantom sample array. Agreement is seen among all three LUTs for both absorption and reduced scattering maps [Fig. 2(a)]. ROIs were chosen to quantify the results for each phantom well [Fig. 2(a)] and to compare the LUT methods [Fig. 2(b)]. All methods show low relative error compared to the MC results $({\mu }_{\mathrm{a}}\%,{\mu }_{\mathrm{s}}^{\prime }\%)$: linear OP LUT (1.3, 0.5), linear $R_d$ LUT (0.2, 0.09), 2-D fit LUT (1.9, 2.8). Table 1 summarizes the processing times for each method on this sample, from a single pixel to the entire image. The linear $R_d$ and 2-D fit LUTs exhibit better performances over the entire range. Importantly, processing shows a 100-fold improvement for an entire image ($572\times 612\text{pixels}$), for both absorption and reduced scattering.

Table 1

Processing speeds for LUT inversion of Rd to both μa and μs′ property maps.

The in vivo results are shown in Fig. 3. All LUT methods show strong agreement for both absorption and reduced scattering maps [Fig. 3(a)]. A line profile was taken across the sample [dotted lines in Fig. 3(a)] to quantitatively compare the LUT results [Fig. 3(b)]. All LUTs show agreement with the MC results $({\mu }_{\mathrm{a}}\%,{\mu }_{\mathrm{s}}^{\prime }\%)$: linear OP LUT (4.1, 2.1), linear $R_d$ LUT (0.9, 0.6), 2-D fit LUT (1.5, 3.2).

Both newly developed LUT methods, the linear $R_d$ LUT and 2-D fit LUT, decrease the processing time by a factor of $\sim 100\times$ by avoiding searching and interpolating with two different approaches. The linear $R_d$ LUT parameterizes a linear grid using the table’s sampling resolution and range and rounds to the nearest vertex that matches the sample’s $R_d$ input. The 2-D fit method avoids this discrete sampling by fitting a continuous function to each surface of ${\mu }_{\mathrm{a}}$ and ${\mu }_{\mathrm{s}}^{\prime }$ as a function of diffuse reflectance at two spatial frequencies.

Each LUT method has its own unique relationship among inversion speed, memory space, and accuracy. The standard linear OP LUT must use interpolation for practically sized tables, which increases the required processing time. Because this table samples $R_d$ in a nonlinear manner, increasing the sampling grid greatly increases the searching time for each input. In comparison, the linear $R_d$ LUT is sampled with a parameterized grid, thus its speed is minimally affected by increased grid size because no searching occurs. However, because the table is linearly sampled, the density required to maximize accuracy at high gradients in the table is used for the entire table, resulting in a large grid that requires considerable computer memory. This could potentially be mitigated by choosing a sampling function whose density increases for higher gradients of the table. The 2-D fit LUT requires less computer memory space but has variable accuracy across the table due to under-fitting, though this accuracy could potentially be weighted toward optical property values that are more likely to occur, e.g., separate LUTs for brain, liver, or skin tissues. While the presented 2-D fit LUT was fit to the entire linear $R_d$ LUT, accuracy was weighted to prioritize optical properties reasonably expected for skin tissues (${\mu }_{\mathrm{a}}=[0.00350.192]{\mathrm{mm}}^{-1}$, ${\mu }_{\mathrm{s}}^{\prime }=[0.32.281]{\mathrm{mm}}^{-1}$). Within this range of optical properties, errors in comparing the 2-D fit LUT with ground truth does not exceed $0.004{\mathrm{mm}}^{-1}$ in absorption and $0.1072{\mathrm{mm}}^{-1}$ in scattering.

The new frameworks introduced can be applied to other LUT processes with possibly more dimensions. For profile corrected SFDI¹¹ or three-dimensional-SSOP,⁶ there is the potential to include height/phase dependence among the table’s input parameters for optical property mapping. In addition, it is important to note that to make multispectral measurements, each wavelength is processed separately with the same LUT. Because wavelengths can be independently processed in parallel, processing time is minimally affected. Given a known set of wavelengths, one could also go directly from $R_d$ input values to hemodynamic values such as oxygen saturation. Generating these calculations after calibrating the imaging system but before the sample acquisition makes real-time feedback for surgical guidance possible.

In this letter, we introduced a new framework for LUT formation and evaluated its performance on tissue-mimicking phantoms and in vivo in comparison with a standard OP LUT.³ Overall, the new techniques were accurate compared to MC simulations, i.e., within 0.9% and 0.6% for linear $R_d$ LUT ${\mu }_{\mathrm{a}}$ and ${\mu }_{\mathrm{s}}^{\prime }$, respectively, and within 1.9% and 3.2% for 2-D fit LUT ${\mu }_{\mathrm{a}}$ and ${\mu }_{\mathrm{s}}^{\prime }$, respectively, and are $100\times$ faster than the standard linear OP LUT. These techniques help enable real-time image-guidance feedback for spatial frequency domain techniques.

The authors would like to thank Steve McClure for his helpful discussion and support. This work was supported by R01-DE-022820 (BL), K01-DK-093603, France Life Imaging, University of Strasbourg and ICube Laboratory (S.G.), and NIH/NIDDK Award Number F31-DK-105839 (J.A.).