2.1.

Instrumentation

We have constructed a wide-field imaging system with a digital micro-mirror device for illuminating the tissue at different spatial frequencies and a spectrally-separated dual camera detection system, which records simultaneously bright-field reflectance images at the excitation wavelength of 808 nm and the fluorescence images at the emission wavelength of 850 nm. The setup is illustrated schematically in Fig. 1(a). The combined illumination and imaging system was mounted on a translation stage (MN10-015-E01-13, Velmex Inc., Bloomfield, New York) and positioned above the sample with the object distance of 32 cm. The Illumination was provided by a digital light projector (FC E4500 MKII, EKB Technologies Ltd. Bat-Yam Israel). Its intensity on the sample surface was $\sim 0.3\mathrm{mW}/{\mathrm{cm}}^2$ for uniform (spatially unmodulated) illumination, well below the American National Standards Institute (ANSI) limit for tissue damage. A digital micromirror device was electronically controlled to produce spatially-modulated (at different frequencies) illumination of the sample surface with wide-field light at the wavelength of 808 nm. The reflected light at the excitation and fluorescent wavelengths was directed to both camera channels and was further spectrally filtered to record the fluorescence image. The images recorded in the fluorescence camera are used for reconstruction of the depth and the transverse margin of the fluorescent inclusions.

Fig. 1

(a) Schematics of the dual-camera imaging system. SM, Spatial modulator of the excitation light; CL, collection lens; RL, relay lens; BS, 10:90 beam-splitter; BF, bright-field camera; EM, emission filter; FL, fluorescence camera. (b) Tissue-simulating liquid phantom and a 3D printed hollow cube ($1{\mathrm{cm}}^3$) containing a mixture of ICG and a scattering solution. For the single-target experiments, the cube was suspended in a tissue-simulating liquid phantom using two posts (blue), and the target depth was adjusted by varying the surface level of the background phantom. (c) For one set of two-target experiments, two identical targets were placed at the same depth below the surface; the transverse distance between the targets was controlled, and the target depth was adjusted as in (b). (d) For the second set of two-target experiments, two targets were separated both horizontally and vertically. The right target was held at 2 mm below the surface, and the depth of the left target was adjusted.

2.2.

Experimental Procedure

The background tissue was simulated using a mixture of intralipid (Intralipid 20%, Fresenius Kabi, Pune, Maharashtra, India), nigrosin (Acid Black, MP Biomedicals, Santa Ana, California), and water. The background has wavelength-dependent absorption coefficient ${\mu }_a$ and reduced scattering coefficient ${\mu }_s^{\prime }$. In choosing optical properties, we adopted generic tissue (breast, muscle, and brain) scattering coefficients, and we chose an absorption coefficient similar to breast tissue. Respectively, these coefficients are 0.004 and $0.8{\mathrm{mm}}^{-1}$ at the excitation wavelength of 808 nm, and 0.006 and $0.76{\mathrm{mm}}^{-1}$ at the emission wavelength of 850 nm. The hollow cube was 3D printed with inner side length of 10 mm and wall thickness of 1 mm (VeroWhitePlus, Stratasys Direct Manufacturing, Valencia, California). It was employed to simulate the fluorescent target(s); the cube was filled with a mixture of intralipid containing dissolved ICG powder (IC-GREEN, Akorn Inc., Lake Forest, Illinois). $6.5\mu \mathrm{M}$ of ICG concentration was chosen to maximize fluorescence yield. The absorption property of background tissue and the ICG concentration were both selected to accommodate our very low-intensity excitation light. The fluorescent contrast cubes were suspended in the liquid phantom tank using fishing lines and posts, as shown in Fig. 1 [Panels (b), (c), and (d)]. The depths of the targets were controlled by lowering or raising the surface of the background liquid phantom in the single-target experiments [Fig. 1(b)]. The image was recorded with 12-bit dynamic range, and the exposure time was increased until the maximum count of a pixel reached approximately 4000 for unmodulated incident light (setting it above 4000 occasionally led to saturation of a few of the sensors.) The fluorescence camera exposure time was adjusted from 800 ms for the 2 mm depth of the target to 5 s for the 10 mm depth to account for the difference in detected intensity. The exposure time of the bright-field camera was about 15 ms. Accordingly, the full acquisition time for the deepest target was four times longer than that of the shallowest target. Fluorescence emission data for single-target phantoms were acquired for target depths of 2, 4, 6, 8, and 10 mm.

For two-target phantoms, two different sets of experiments were performed. In one set, two identical targets were placed at the same depth. This depth of the two targets, and the lateral separation between the targets, were varied as shown in Fig. 1(c). In a second set of two-target phantom experiments, the targets were laterally separated by either 20 or 40 mm. The depth of one target was fixed at 2 mm, whereas the depth of the other target was varied from 4 to 8 mm [Fig. 1(d)].

In all experiments, the data acquisition time for the fluorescence associated with the shallowest target (2 mm depth) was $\sim 2\min$. The acquisition time for fluorescence associated with the deepest target (10 mm depth) was $\sim 8\min$.

2.3.

Data Processing and Analysis

Theory of the spatial frequency-domain imaging has been exposed in Refs. 19 and 32. For clarity, we provide here a concise description of the relevant ideas, largely following Ref. 19. The symbols and corresponding units for all parameters and variables are summarized in Table 1.

Table 1

Descriptions of parameters and variables.

First, we describe the illumination scheme. The sample surface is illuminated by a sinusoidally-modulated intensity pattern of the form

Fig. 2

Example of the measurement in case of single target located at the depth of 2 mm ($f=0.04{\mathrm{mm}}^{-1}$). Scale bar indicates 10 mm.

For each illumination pattern $S_i({\mathbf{\rho }}_d;\mathbf{k})$, a corresponding fluorescent emission image ${\mathrm{\Phi }}_i({\mathbf{\rho }}_d;\mathbf{k})$ was recorded (the top row of Fig. 2). The images ${\mathrm{\Phi }}_0$, ${\mathrm{\Phi }}_1$, and ${\mathrm{\Phi }}_2$ were then combined to obtain the complex fluorescent emission signal, $\mathrm{\Phi }({\mathbf{\rho }}_d;\mathbf{k})$ according to

Image reconstruction relies on an integral relation between $\widehat{\mathrm{\Phi }}({\mathbf{q}}_d;\mathbf{k})$ and the heterogeneous fluorophore concentration in the sample, assuming uniform background optical properties and first-order approximation for diffuse wave propagation, which linearizes the inverse problem. The relation is of the form

The first step in the reconstruction algorithm is to estimate the target depth. To this end, we utilize a simplification of Eq. (4). Specifically, we assume that the fluorophores are localized at the depth of the top surface of each target, which we refer to as $z_{\text{target}}$, so that $\widehat{C}({\mathbf{q}}_d-\mathbf{k},z)=\widehat{C}({\mathbf{q}}_d-\mathbf{k})\delta (z-z_{\text{target}})$. In this case, Eq. (4) simplifies to the following form (with $\mathbf{q}={\mathbf{q}}_d-\mathbf{k}$):

The dependent variable, $y(\mathbf{k})$, is the normalized emission response at the modulation wave vector $\mathbf{k}$. The independent variable, $x(\mathbf{k})$, can be measured or estimated; it is, essentially, twice the reciprocal of the penetration depth, which depends on the spatial modulation frequency. The parameter $A=\epsilon \eta (\frac{{\ell }_{\mathrm{ex}}{\ell }_{\mathrm{em}}}{D_{\mathrm{ex}}{D}_{\mathrm{em}}})\widehat{C}(0)$ depends on the parameters of the background medium and the fluorophores, some of which may not be directly known or measured, but it is independent of the modulation wave vector $\mathbf{k}$. In a data set obtained by varying $\mathbf{k}$, we can view $A$ as a constant and as an adjustable parameter. Further, $y_o$ is an extra adjustable parameter added to account for the noise floor in the data since the fluorescent signal does not go to zero when $|\mathbf{k}|\to \infty$ (as predicted theoretically), perhaps as a result of background fluorescence and other noise, as well as imprecision in the theoretical model. Note that the adjustable parameters $A$ and $y_o$ are expected to be different in different data sets, but in each deta set they are independent of $\mathbf{k}$. We can now find $z_{\text{target}}$ by nonlinear fitting of the theoretical Eq. (8a) to the data in which $A$, $z_{\mathrm{target}}$ and $y_o$ are viewed as $\mathbf{k}$-independent adjustable parameters.

Once the depth of the target is determined, we can compute the lateral margins of the fluorophore concentration in the plane. Rearranging Eq. (7) for $\widehat{C}(\mathbf{q})$ and taking the inverse 2D Fourier transform, we obtain a simple reconstruction formula for the transverse distribution of the fluorophore target in the plane $z=z_{\text{target}}$. The obtained relation is parameterized by $\mathbf{k}$ and, theoretically, any value of $\mathbf{k}$ can be used to obtain the transverse distribution of the fluorophore. In the reconstructions, we have used the $\mathbf{k}=0$ for this purpose because the signal-to-noise ratio in the fluorescence images is best for unmodulated incident light. Therefore, we use the reconstruction formula

2.4.

Nonlinear Fitting

Prior to fitting, the data were normalized to the maximum value. That is, both sides of Eq. (8a) were divided by $y(\mathbf{k}=0)$. This allowed us to compare errors of the fit for different data sets quantitatively. We then obtained an estimate of the fluorophore target depth $z_{\text{target}}$ from the best fit of the theoretical Eq. (8a) to the data points [$x(\mathbf{k}),y(\mathbf{k})/y(0)$]. We removed the first data point (with the smallest $|\mathbf{k}|$) from the fitting procedure for the reasons explained in Sec. 3.1. For the fitting, we used a nonlinear regression model algorithm, fitnlm (MATLAB 2022a, The Mathworks Inc., Natick, Massachusetts), and chose the initial values for $A$, $z_{\text{target}}$, and $y_o$ to be 5.0, 1.0, and 0.1, respectively. For verification, we also used the fitting function of Gnuplot. Very similar results were obtained (these data are not shown below).

3.1.

Single-Target Experiments

The depth of the fluorescence contrast target was varied from 2 to 10 mm, and the corresponding exponential fitting results are displayed in the left column of Fig. 3. Note that, as mentioned above, the first data point corresponding to $k=0$ was not used for fittings. If used, this data point can skew the decay constant to larger values and increase the fitting errors significantly. Removal of this point can be justified on grounds that the fluorescent targets employed in our experiments are not infinitely thin layers but cubes and therefore have finite depth of 10 mm. Owing to the high absorption coefficient of ICG and the small Stokes shift, we perform the reconstructions under the assumption that the fluorescent signal comes predominantly from the top faces of the cube. This assumption, however, disregards penetration of the incident photons in the interior of the fluorescent targets. As a result of this effect, the inverse problem of reconstructing the ICG concentration is nonlinear, and replacing it with a linearized problem, as is done in our theoretical analysis, is an approximation. The quality of this approximation is expected to be especially poor at $k=0$, which explains the anomalous position of the respective data point in Fig. 3. At higher values of $k$, the effects of nonlinearity appear to be weaker and can be neglected. However, a more precise corroboration of this picture will require additional (future) experiments with smaller fluorescent targets or smaller ICG concentration; if successful, inclusion of the $k=0$ (and other data points with relatively small values of $k$) will be possible and is expected to increase the dynamic range of measurements, which is useful for more accurate depth estimation.

Fig. 3

Single target experiment results. The target depth is (a) 2 mm, (b) 4 mm, (c) 6 mm, (d) 8 mm, and (e) 10 mm. The left column shows the exponential fitting of the normalized data used to estimate the target depth. The estimated depth is marked on top of the plot along with its standard deviation. The middle column shows the reconstructed transverse images. The right column shows the transverse margin from the image and the true margin. The yellow scale bar corresponds to 10 mm.

After the target depth and its uncertainty are estimated, we determine the transverse margins of the fluorophore concentration by using Eq. (9). The resulting reconstructions depend on the regularization parameter $\sigma$. If $\sigma$ is too large, the resulting images typically contain significant random noise, which renders the target unrecognizable [see Fig. 4(a)]. On the other hand, if $\sigma$ is too small, the reconstructed image shows only a blob [see Fig. 4(c)]. To determine the optimum $\sigma$, we apply the weak positivity constraint on $C(\mathbf{\rho },z=z_{\text{target}})$. Since $C$ is intrinsically positive, we require that the mean value of the normalized concentration over all pixels is at least one standard deviation larger than zero; this roughly corresponds to the condition that the fraction of pixels holding negative values is less than 16%. The optimal $\sigma$ is the maximum $\sigma$ satisfying the constraint. Note, the regularization parameter, $\sigma$, does not impose any assumptions about the spatial variation of fluorophore concentration nor symmetry of the inclusions. The optimized image is shown in Fig. 4(b), and the transverse images of the middle column of Fig. 3 are produced from the same positivity constraint.

Fig. 4

Transverse image of the fluorescent inclusion at 4 mm depth with different values of $\sigma$ using the same data as in Fig. 3(b). (a) The top and bottom rows show the cases of weak regularization corresponding to $\sigma =3.77{\mathrm{mm}}^{-1}$ and (b) strong regularization corresponding to $\sigma =0.11{\mathrm{mm}}^{-1}$, The middle row was obtained using the optimal regularization (as defined in the text) with $\sigma =0.46{\mathrm{mm}}^{-1}$. In the optimized images, the number of pixels with negative values is relatively small ($<16\%$ of all pixels). The left column shows the reconstructed images obtained from Eq. (9) and the right column provides normalized concentration distribution histogram. The scale bar corresponds to 10 mm.

The lateral margin of the target was set to be the full-width-at-half-maximum (FWHM) contour line of this distribution (the last column of Fig. 3). To measure the discrepancy between the reconstructed and true margins, we computed the relative width—defined as the square root of the ratio of transverse target area (area confined by the blue line in the last column of Fig. 3 to the true target area ($100{\mathrm{mm}}^2$). For the single-target experiments, the reconstructed depth and the relative width are fairly accurate and are summarized in Fig. 5.

Fig. 5

Reconstruction results of single-target experiments. (a) Depth-estimation results. The black dotted line has a unity slope. (b) The ratio between reconstructed image width and the true target width as a function of depth; margin estimate error varies from 20% to 33%.

3.2.

Two-Target Experiments

We first placed two identical targets in the same plane but with varying lateral separation. The estimated depths and transverse margins were then computed. A few examples are shown in Fig. 6(a) The two targets are clearly distinguished in most cases. As expected, the ability to distinguish between the targets is reduced as the depth is increased. However, for the range of separations used in our experiments, we clearly observe two targets at every depth. A summary of depth sensitivity and relative width is given in Figs. 6(b) and 6(c).

Fig. 6

Reconstruction results of two targets in the same vertical plane. (a) Reconstructed transverse margin. When the separation increases, the cross-talk between the targets is decreased. (b) Depth-estimation results. (c) Relative width estimation results. For clarity, we spread the data horizontally so that a “cluster” of points arises at each target depth.

The second, more challenging problem is to reconstruct two targets that are separated both vertically and laterally. We fixed the lateral separations to be either 20 or 40 mm. The vertical separations were set to be 2, 4, and 6 mm; the depth of Target 2 was fixed at 2 mm below the surface and the depth of Target 1 was varied taking the values 4, 6, and 8 mm. Since the targets were well separated transversely, even when they are located at different depths, we can safely segment the data to process each target separately (Fig. 7). We use the minimum in the peak intensity profile along the horizontal (dotted) line to define the segmentation regions for the image as is shown by the vertical (dashed) lines in figure. After the depth of each target was independently estimated, we used these estimates to process the data further for each target. We obtained the transverse fluorophore distributions by utilizing the same algorithm as in the single target reconstructions.

Fig. 7

Raw data showing substantially different emission signals from targets 1 and 2. The more deeply located target is target 1, and the more superficial target is target 2. The intensity profile is plotted along the dotted line A. The vertical dashed line indicated by B defines the two segmented regions. The scale bar corresponds to 10 mm.

Reconstructions with 20 and 40 mm separations between the targets are presented in Fig. 8(a). The depth-sensitivity and depth estimates for these targets [Fig. 8(b)] are not as good as our results for single targets and for two-targets in the same horizontal plane [Figs. 5(b) and 6(b)]. This can be understood by noting that both the experimental acquisition time and the resulting computationally processed data are different for different target depths. As described in Sec. 2, the camera exposure times needed to obtain reliable fluorescence images for deeper targets are significantly larger than those for shallow targets. In our instrument and data collection scheme, it is not possible to use large exposure times when a shallow target is present due to the possibility of overexposure. We were, therefore, forced to use the exposure times appropriate for the more shallow target. Thus, when target 1 is located at the depth of 10 mm, its fluorescent signal was not detectable due to the insufficient exposure time. Cross-talk of the targets located at different depths can also cause artifacts. As expected, the errors of the transverse margins of the target are larger for smaller horizontal separations [Fig. 8(c)]. When the transverse separation is 20 mm, the cross-talk from target 2 in the segmented image of target 1 becomes significant and produces an artifact near the edge. When this effect is severe, the transverse margins of the deeper target are corrupted [see last column of the first row of Fig. 8(a)]. In the case of the middle column of the first row of Fig. 8(a), the overestimation of the depth is significant. As a result, the two-point Green’s function kernel [Eq. (5)] decays faster as a function of $q$, which forces strong regularization for transverse margin reconstruction. The strong regularization suppresses the bleeding artifact and produces an image of a nearly featureless large blob.

Fig. 8

Results of the two-target experiments with vertical separations of the targets. (a) 3D rendering of the results. Columns and rows represent the depth of target 1 and the transverse separation of the targets, respectively (Target 2 is fixed at the 2 mm depth). Estimated depths are marked below each margin slice. The yellow scale bar represents a 10 mm length in the transverse plane. The vertical scale is exaggerated for better visibility. (b) Results for the reconstructed depth of targets 1 and 2. The dashed line has the unity slope. (c) Results of the relative width of targets 1 and 2.

4.1.

Depth Sensitivity

The simple algorithm proposed in this paper successfully estimated the depth of fluorescent inclusions to within $\sim 1\mathrm{mm}$ of their true depth in both single target and the two target experiments in which the targets were at the same depth but laterally separated. Prior published work by Konecky et al.¹⁹ and by Li et al.²⁴ also employed the spatial frequency domain method. However, both works used an image reconstruction algorithm based on the pseudoinverse technique. In this algorithm, noise is suppressed by regularization, which often results in a low-contrast image in which the reconstructed fluorophore target stretches further than its true size; for this reason, the accuracy of reconstructed target depth and transverse margins is limited. By contrast, our technique uses a simple exponential relationship between the fluorescence emission signal $k$-dependence (i.e., dependence on modulation frequency) and the target depth; in this approach, the noise floor [$y_o$ in Eq. (8a)] is decoupled from depth information in the exponential decay. Thus, our results are superior to the previous SFD-FDOT simulation results by Li et al.²⁴ wherein the accuracy of reconstructed simulated data degraded with increased target depth. (Note, our experimental depth sensitivity results are comparable to simulated data results reconstructed by neural-network-aided FDOT.^27,29) Interestingly, our estimation error (both absolute and especially fractional) for single targets was minimum at the largest depth (Fig. 9). The finite target thickness introduces nonlinearities into the reconstruction problem that we do not account for; these effects, briefly discussed in Sec. 3.1, were ameliorated by removing the data point $k=0$ in the fitting procedure, but it does not fully resolve the effect for all other values of $k$. However, as the target moves away from the surface the nonlinear effects of target thickness diminish, leading to better depth estimation. In the case of two targets at different depths, the modulation response from more deeply located targets severely suffered cross-talk from the superficially located target. This, and the experimental limitations related to the available exposure times, produced a larger error in terms of depth estimation in this set of experiments, typically, of the order of 20% to 40% of the true depth.

Fig. 9

Mean absolute errors for all depth estimation results.

Clearly, some technical improvements are possible in the future work. Theoretically, better accounting for the effect of finite target thickness will potentially allow more precise fitting of depth. This will increase the dynamic range of useful data and reduce the systematic errors of the theoretical model. In addition to reducing depth estimation errors of shallow targets, such improvements can enable imaging of targets located deeper than 10 mm below the surface.

4.2.

Transverse Margins

Another important advance is the improved ability to constrain the transverse margins of the fluorescent inclusions (tumors). From the single target experiments, we compared the transverse margin estimates of our reconstructions of fluorophore concentration versus those of traditional 2D fluorescence projection images. (The transverse margins were set to be the FWHM of the image contours.) The relative transverse margin width obtained by each technique is plotted in Fig. 10 as function of target depth. Notice, the traditional 2D fluorescence projection images increasingly overestimate the transverse margins as the target depth becomes larger. For example, at $\sim 1\mathrm{cm}$ target depth, the error of the 2D projection image margins was roughly 2.5 times larger than of the margins obtained by our method. The transverse margins of the single targets, or two targets in the same plane, mostly overestimated the true margins by 30% or less. For the two-target experiments with targets at different depth, however, these margin errors were larger, in part because the larger error in depth estimation propagates to the error in transverse margin (Fig. 8(b,c)).

Fig. 10

Transverse margin comparison between our reconstruction approach and estimation based on 2D-projection images of fluorescence based on the single-target experiments.

4.3.

Spatial Smoothing (Regularization) Parameter

Figure 11 shows the average value of the optimized regularization parameter $\sigma$ as a function of the target depth for single-depth reconstructions. It can be seen that the regularization parameter tends to decrease as the target depth increases. This is because the signal-to-noise ratio for the detected fluorescent intensity decreases as the pathlength of light increases and the high spatial frequencies in the detected intensity distribution become dominated by noise. Empirically, we note that the value of $k$ for which $y(\mathbf{k})$ settles into the noise plateau is roughly proportional to the optimal $\sigma$ (roughly, $k\sim \sigma$). Therefore, the value of $k$ settling into the noise plateau can be considered as a guideline for choosing the optimal $\sigma$.

Fig. 11

Optimized regularization parameter $\sigma$ for single-depth reconstruction as a function of the target depth.

4.4.

Stability of the Estimates

The proposed two-step reconstruction method utilizes Green’s function formalism, which is affected by the optical properties of the background. Therefore, errors in the background optical properties are expected to propagate to errors in estimation of the depth and transverse margin of the fluorescent target. We tested the stability of the reconstructions by introducing a $\pm 10\%$ error into ${\mu }_a$ and ${\mu }_s^{\text{'}}$. The resulting absolute error of depth estimation was $<0.8\mathrm{mm}$, which does not significantly impact the overall depth-sensitivity. Similarly, the change in transverse margins for the targets at the depths of 2, 4, and 6 mm was $<2\%$. Thus, while it is desirable to use the most accurate background optical properties, our method appears to provide fairly stable results for optical property variations of order 10% or less.

4.5.

Data Acquisition and Reconstruction Times

Data acquisition time introduces the longest temporal delays in our experiments, typically, a few minutes. This time can be reduced. For example, the 31 distinct spatial frequencies currently used are probably more than necessary for good quality imaging. If we remove every second data point from the left panel of Fig. 3, the quality of nonlinear fit is not significantly affected. In this scenario, the image quality is the same, and the data acquisition time is reduced by a factor of two. On the other hand, in some situations additional data points can help improve signal-to-noise. Future work is needed to optimize the number of spatial frequencies for specific instruments and instrumentation parameters.

Since our technique does not attempt to reconstruct full 3D-shapes, the time needed to estimate the depth and transverse margins ranged from 1 to 2 s (Intel(R) Core(TM) i5-8600K CPU 3.60 GHz, 6 Cores). Clearly, the approximate approach is different from the traditional analytic inversion and nonlinear image reconstruction techniques, which require computation times ranging from minutes to hours depending on complexity (i.e., for traditional DOT techniques). Recently, advances based on more complex techniques, such as confocal time-of-flight FDOT and neural network aided FDOT, have demonstrated improvement in reconstruction runtime (from a few milliseconds to seconds) and spatial resolution (millimeter scale). While these advances are excellent and may surpass our reconstruction speed, their clinical suitability have to be evaluated.

4.6.

Limitations

The proposed SFD-FDOT estimation methodology is promising albeit with the aforementioned limitations. Notably, reconstruction of absolute fluorophore concentration was not pursued, but absolute concentration is not a feature currently employed for image guidance or diagnosis. Per image guidance, the analysis assumed homogeneous background optical properties, semi-infinite geometry, and a thin slab geometry of the fluorescence inclusion. In practice, heterogeneous tissue optical properties in brain or lung tissue could generate errors in estimates of depth and transverse margin. Even though this assumption is commonly used in diffuse optics analysis, more in vivo work needs to be done to fully characterize these limitations. Per the semi-infinite geometry, we note that this approximation has proven adequate for many human and animal studies in the diffuse optics field. Therefore, we expect that the proposed technique will be applicable in many open surgery cases for neuro- or thoracic surgery. In principle, this instrumentation can be modified to concurrently measure surface profiles and thereby deduce the magnitude of the deviation from the assumption; this will also provide concrete data that could be employed to modify the current approach perturbatively to include corrections due to surface curvature. Per the thin slab approximation for the fluorescence inclusion, in practice, the high absorption coefficient of ICG will cause the top surface of the target to absorb most of the excitation light propagating downward, and consequently the fluorescence emission will also be dominated by the signal from the top surface. The thin slab approximation is useful, as long as the user realizes that the target depth corresponds to the top surface of the fluorescent region. Note also, even if tumor tissue is not flat on the top surface, the proposed technique will provide a useful estimate of depth that is skewed toward the shallowest part of the tumor.

Finally, it is desirable to make measurements quickly in the operating room. In our current setup, the longest part of the procedure is data acquisition (up to 8 min for the deepest occlusions). In the future, this can be improved by not using the unnecessary spatial modulation frequencies (at high-$q$) and increasing the excitation light intensity. Notably, the current light intensity on the sample surface was smaller than the ANSI limit by the factor of 3 (see Appendix A of Ref. 34).