3.1.

Phase Space Representation

A discretized phase space is required to evaluate the radiances using the discrete segments of the photon trajectories generated by the MC simulation. However, increasing the granularity of the phase space may not only increase the variance of the radiance estimates, due to a reduction in the sample size for a given bin, but also requires prohibitive amounts of memory. To alleviate the latter concern, we restrict our determination of the forward-adjoint radiances to a ROI that occupies a subvolume of the tissue domain in which the photons propagate. The implementation of the ROI allows for the computation and coupling of the forward-adjoint light fields with high granularity, while allowing photons to propagate throughout a much larger domain that need not be discretized.

Figure 1 depicts a ROI in the form of a rectangular solid, located within a volume bordering the surface of an optically thick tissue. The optical probe employed to interrogate the tissue provides light from a normally incident source and uses two detectors located at the same location but with different orientations. The ROI is composed of a uniform 5-D spatial-angular mesh with volume elements (voxels) and solid angle elements of dimension $\mathrm{\Delta }x\mathrm{\Delta }y\mathrm{\Delta }z=\mathrm{\Delta }{l}^3$ and $\mathrm{\Delta }\mu \mathrm{\Delta }\varphi$, respectively. The total number of discrete bins contained in the ROI is $N_{\mathrm{tot}}=N_x{N}_y{N}_z{N}_{\mu }{N}_{\varphi }$, where $N_{(·)}$ is the number of elements for each phase space variable. The volume of the ROI is $l_x{l}_y{l}_z$, where the lengths are given by $(l_x,l_y,l_z)=(N_x\mathrm{\Delta }l,N_y\mathrm{\Delta }l,N_z\mathrm{\Delta }l)$, and the subvolume of tissue that the ROI encloses is further defined by the minimum vertex ${\mathit{r}}_{\min }$. We denote each voxel in the ROI as $\mathrm{\Delta }\mathit{r}=(i,j,k)\in \{[1,N_x],[1,N_y],[1,N_z]\}$, which encloses the spatial locations $\mathit{r}=(x,y,z)\in \{[x_{\min }+(i-1)\mathrm{\Delta }l,x_{\min }+i\mathrm{\Delta }l],[y_{\min }+(j-1)\mathrm{\Delta }l,y_{\min }+j\mathrm{\Delta }l],[z_{\min }+(k-1)\mathrm{\Delta }l,z_{\min }+k\mathrm{\Delta }l]\}$. The photon direction is conveniently mapped to the two angular variables specified on the unit sphere as described in Sec. 2, where $\mathrm{\Delta }\mu =2/N_{\mu }$ and $\mathrm{\Delta }\varphi =2\pi /N_{\varphi }$. Each solid angle element is then specified as $\mathrm{\Delta }\widehat{\mathit{s}}=(u,v)\in \{[1,N_{\mu }],[1,N_{\varphi }]\}$, which spans the photon propagation directions with angular variables $\mu \in [(u-1)\mathrm{\Delta }\mu ,u\mathrm{\Delta }\mu ]$ and $\varphi \in [(v-1.5)\mathrm{\Delta }\varphi ,(v-0.5)\mathrm{\Delta }\varphi ]$. The boundaries of the $\varphi$ bins are shifted by $-0.5\mathrm{\Delta }\varphi$, so that $\varphi$ values that are co-planar with the $x\text{-}z$ and $y\text{-}z$ planes occupy the center of the angular bin. To preserve this feature, we restrict the value of $n_{\varphi }$ to multiples of 4.

Fig. 1

Illustration of the space-angle mesh contained in the region of interest (ROI) that is located within the tissue volume (a). Each volume element (voxel) if further discretized in angle, where $\mathrm{\Delta }\widehat{s}=\mathrm{\Delta }\mu \mathrm{\Delta }\varphi$ is a solid angle element, which is depicted in the closeup view of a single voxel (a). Illustration of the trajectory elements that photon packet $m$ might experience during transport through the tissue (b).

Using the spatial-angular mesh, both the discrete forward-adjoint radiances, $L(\mathrm{\Delta }\mathit{r},\mathrm{\Delta }\widehat{\mathit{s}})$ and $L^{\star }(\mathrm{\Delta }\mathit{r},\mathrm{\Delta }\widehat{\mathit{s}})$, respectively, can be estimated within the ROI. The absorption sensitivity at any given discrete location $\mathrm{\Delta }{\mathit{r}}_{i,j,k}$ is then obtained by coupling these radiances and integrating them over the solid angle elements $\mathrm{\Delta }\widehat{\mathit{s}}$ via a discrete analog of Eq. (4)

In the following section, we provide the algorithmic details of our approach to evaluate the radiances.

3.2.

Radiance Evaluation

To understand how the radiance is evaluated from the detailed information of a CAW-MC simulation, we must consider the specific characteristics of ${\mathbb{T}}_{m,n}$, which records information pertaining to the $n$’th segment of the $m$’th photon packet and specifies its location ${\mathit{r}}_{m,n}$, weight $W_{m,n}$, and direction ${\widehat{\mathit{s}}}_{m,n}$ from which the photon moves a length ${\ell }_{m,n}$ to location ${\mathit{r}}_{m,n+1}$ with weight $W_{m,n+1}=W_{m,n}\exp (-{\mu }_a{\ell }_{m,n})$. It is from the location specified by ${\mathbb{T}}_{m,n+1}$ that the photon begins the next segment. Using the trajectories that result from $N_{h\nu }$ photon packets, the radiance is tallied for each segment and the photon traces within the ROI. For photon segments that take place within a single voxel, the radiance tally can be computed using

According to Eq. (7), the $n$’th segment of the $m$’th photon packet deposits a weight fragment $\mathrm{\Delta }{W}_{m,n}=W_{m,n}-W_{m,n+1}=W_{m,n}[1-\exp (-{\mu }_a{\ell }_{m,n})]$ in the space-angle element $[\mathrm{\Delta }{\mathit{r}}_{i,j,k},\mathrm{\Delta }{\widehat{\mathit{s}}}_{u,v}]$, and the ratio of the deposited weight to the absorption coefficient is tallied to the radiance within the corresponding spatial.

However, Eq. (7) does not consider photon segments ${\mathbb{T}}_{m,n}$ that span multiple voxels. This case is important because accurate MC estimation of the radiance distribution requires a spatial discretization with voxel dimension $\mathrm{\Delta }l$ that is comparable to or less than the scattering length ${\ell }_s=1/{\mu }_s$.³¹ In CAW, the segment length that the photon packet traverses between its $n$’th and $(n+1)$’th collisions is sampled from ${\ell }_{m,n}=-\ln (\xi )/{\mu }_s$, where $\xi$ is a random value between 0 and 1. Thus, it is likely that many photon segments ${\ell }_{m,n}$ will span multiple volume elements. To obtain radiance estimates in this case, we must deposit the appropriate weight fragments $\mathrm{\Delta }{W}_{m,n,p}$ in each of the $p$ voxels traversed by photon segment ${\mathbb{T}}_{m,n}$. The consideration of this detail requires great care when using CAW-MC to obtain accurate internal RTE quantities. To account for these events, Eq. (7) is reformulated as follows:

Equation (8) provides the radiance by properly considering each voxel that is traversed by the photon segment ${\mathbb{T}}_{m,n}$. However, this requires a method to compute each of the individual weight fragments $\mathrm{\Delta }{W}_{m,n,p}$ and the corresponding voxel indices $(i,j,k)$ for each subsegment. Moreover, we must map the photon direction ${\widehat{\mathit{s}}}_{m,n}$ for segment ${\mathbb{T}}_{m,n}$ to the proper solid angle bin $(u,v)$. To efficiently evaluate $\mathrm{\Delta }{W}_{m,n,p}$ required for Eq. (8), we utilize the intersection(s) of the photon segment ${\mathbb{T}}_{m,n}$ with the spatial mesh of the ROI. Using these segment-mesh intersection locations ${\mathit{r}}_{m,n,p}$, we evaluate a series of cumulative segment lengths ${\ell }_{m,n,p}=\Vert {\mathit{r}}_{m,n,p}-{\mathit{r}}_{m,n}\Vert$, relative to ${\mathit{r}}_{m,n}$. This allows the photon weight at the segment-mesh intersections to be determined as $W_{m,n,p}=W_{m,n}\exp (-{\mu }_a{\ell }_{m,n,p})$. Finally, the evaluation of the deposited weight fragments $\mathrm{\Delta }{W}_{m,n,p}$ requires the successive photon weights, that is, $W_{m,n,p+1}$ and $W_{m,n,p}$. Therefore, if the lengths ${\ell }_{m,n,p}$ are determined sequentially, the calculation of the deposited weight fragments follows easily.

To efficiently compute both ${\ell }_{m,n,p}$ and the unique $\mathrm{\Delta }{\mathit{r}}_{i,j,k}$ voxel that correspond to weight fragment $\mathrm{\Delta }{W}_{m,n,p}$, we developed an algorithm that directly computes a complete, but unordered, list of cumulative segment lengths ${\ell }_{m,n,q}$. The ordering of this list then provides the proper sequence of cumulative segment lengths ${\ell }_{m,n,p}$. Additionally, details are extracted from each segment-mesh intersection to provide a map of the voxel indices traversed by segment ${\mathbb{T}}_{m,n}$. This approach provides for a more efficient procedure as compared to determining the segments of ${\mathbb{T}}_{m,n}$ by tracing the trajectory of the photon segment through the spatial mesh of the ROI.

Our procedure begins by leveraging the voxel indices that correspond to collision locations ${\mathit{r}}_{m,n}$ and ${\mathit{r}}_{m,n+1}$. These indices are evaluated as

This procedure is then repeated for both the $y$ and $z$ spatial variables to complete the list of ${\ell }_{m,n,q}$. However, as previously described, the evaluation of the cumulative lengths alone does not provide sufficient information to tally the radiance. Therefore, the length ${\ell }_{m,n,q}$ is stored in a structured array along with the unit vector collinear with the photon propagation direction and the mesh intersection spatial dimension. These values provide both the spatial dimension and direction of the segment-voxel intersection. For example, $x$-voxel intersection has a collinear unit vector $[\mathrm{sgn}({\widehat{s}}_{x,m,n}),0,0]$. Next, the qsort algorithm³³ is employed to produce the ordered list ${\ell }_{m,n,p}$ by sorting the structured data array according to ascending cumulative segment length ${\ell }_{m,n,q}$.

Figure 2 provides a 2-D representation of this process in the $x\text{-}z$ plane. Beginning with the photon segment ${\mathbb{T}}_{m,n}$, the unordered cumulative segment lengths ${\ell }_{m,n,q}$ are evaluated along with their collinear unitvectors. These data are first sorted to provide the contiguous cumulative segments ${\ell }_{m,n,p}$. The corresponding weight fragments are then determined and tallied to the appropriate spatial-angular bin. The first weight fragment is denoted as $\mathrm{\Delta }{W}_{m,n,p=0}=W_{m,n}-W_{m,n,p=1}$ and contributes to $\mathrm{\Delta }\mathit{r}({\mathit{r}}_{m,n})$. The voxel indices for the $p=1$ subsegment is then determined by addition of the collinear unitvector corresponding to ${\ell }_{m,n,p=1}$ to $\mathrm{\Delta }\mathit{r}({\mathit{r}}_{m,n})$. The weight fragment $\mathrm{\Delta }{W}_{m,n,p=1}=W_{m,n,p=1}-W_{m,n,p=2}$ contributes to this updated voxel index, and the process of updating the voxel index and tallying the weight fragments $\mathrm{\Delta }{W}_{m,n,p}$ is repeated until $p=N_n-1$. The last subsegment where $p=N_n$ assigns the weight fragment $\mathrm{\Delta }{W}_{m,n,p=N_n}=W_{m,n,p=N_n}-W_{m,n+1}$ to $\mathrm{\Delta }\mathit{r}({\mathit{r}}_{m,n+1})$.

Fig. 2

Illustration of the algorithm used to account for multiple voxel intersections, shown in 2-D for convenience. Using the segment length ${\ell }_n$, the unsorted cumulative segment lengths ${\ell }_{n,q}$ are determined for each spatial dimension. This list of lengths is sorted to provide the ordered cumulative segment lengths ${\ell }_{n,p}$. Moreover, the details regarding each mesh intersection direction are used to provide a map that traverses the voxel indices for efficient tallying, illustrated by the color-coded intersection points.

This algorithm efficiently evaluates the radiance using CAW, through the proper determination of the subsegment weight fragments $\mathrm{\Delta }{W}_{m,n,p}$ and the assignment to the appropriate spatial-angular element $(\mathrm{\Delta }{\mathit{r}}_{i,j,k},\mathrm{\Delta }{\widehat{\mathit{s}}}_{u,v})$. In cases where the source-detector pair lies in a single plane that is orthogonal to the tissue surface, the forward-adjoint radiances exhibit even symmetry about this plane, and we can exploit this symmetry to reduce the memory requirements of the cFAMC simulation. In this study, we consider cases where the source-detector pair lies in the $x\text{-}z$ plane. In all these cases, the radiances exhibit even symmetry with respect to $y$, that is,

This relationship allows us to replace the physical ROI with a memory-reduced computational ROI, while simultaneously increasing the quantity of photon tallies that contribute to the radiance in any given bin.

The proper utilization of this symmetry relationship reduces the memory requirements by nearly a factor of 2 and improves overall computational efficiency.

While this approach to evaluate the radiance applies to both forward-adjoint photon transfers, it is necessary to distinguish between forward-adjoint simulations. While a separate MC code could be developed for adjoint photon propagation from adjoint sources, this is not necessary. Rather, we leverage the reciprocity relationship of the forward-adjoint radiances²⁴

Following the determination of both forward-adjoint radiances for a given source-detector pair, the spatial distribution of the sensitivity is determined using Eq. (6). Moreover, the radiance resulting from a single forward simulation can be reused to couple with multiple adjoint radiances provided by different detector configurations, that is, changing the detector orientation and/or position.

4.1.

Simulation Details

We compute the forward-adjoint radiances and the resulting absorption sensitivity distributions using four levels of angular discretization. These levels and descriptive names are provided in Table 1. We consider the highest level of angular discretization (fine) to provide the most accurate radiances and associated sensitivity distributions. Thus, we will compare the three other levels of angular discretization to the fine case. The case of no angular discretization is an important one as this corresponds to a calculation where the forward-adjoint light fields are modeled as a scalar density function without any angular variation. To obtain the results shown below, we launched ${10}^8$ photons or ${10}^9$ photons for each forward-adjoint simulation in the cases of highly (${\mu }_s^{\prime }/{\mu }_a=100$) or moderately scattering media (${\mu }_s^{\prime }/{\mu }_a=1$), respectively. More photons were required in the moderately scattering case since the increased absorption and decreased scattering led to a reduced angular dispersion of the photon propagation direction. Additionally, we use the technique of Russian Roulette to reduce the variance and runtime of the simulations,³¹ where photons with a weight less than 0.0001 have a 10% chance of survival.

Table 1

Levels of angular discretization used in cFAMC simulations.

Simulations that employ a source-detector separation ${\rho }_{\mathrm{sds}}={\ell }^{\star }$ are performed with a voxel dimension of $\mathrm{\Delta }l=0.02{\ell }^{\star }$, and the physical ROI has dimensions of $1.8{\ell }^{\star }\times 1.38{\ell }^{\star }\times 1.8{\ell }^{\star }$ for a total volume of $4.47{\ell }^{{\star }^3}$. This results in a computational ROI with a total number of bins $N_{\mathrm{tot}}=n_{\mu }{n}_{\varphi }\times 90\times 35\times 90=n_{\widehat{\mathit{s}}}\times 2.835\times {10}^5$. Simulations using the larger source-detector separation ${\rho }_{\mathrm{sds}}=3{\ell }^{\star }$ employ a larger voxel size of $\mathrm{\Delta }l=0.03{\ell }^{\star }$ to accommodate the larger ROI of $5.4{\ell }^{\star }\times 4.11{\ell }^{\star }\times 5.4{\ell }^{\star }$ with volume $120{\ell }^{{\star }^3}$. This provides a computational ROI with $N_{\mathrm{tot}}=n_{\mu }{n}_{\varphi }\times 180\times 69\times 180=n_{\widehat{\mathit{s}}}\times 2.2356\times {10}^6$. Given the voxel dimensions and the ROI, the memory required for a single radiance estimate, as a function of the solid angle discretization, is $n_{\widehat{\mathit{s}}}\times 2.268\mathrm{MB}$ for ${\rho }_{\mathrm{sds}}={\ell }^{\star }$ and $n_{\widehat{\mathit{s}}}\times 17.9\mathrm{MB}$ for ${\rho }_{\mathrm{sds}}=3{\ell }^{\star }$, using double precision floating point numbers.

5.1.

Highly Scattering Media ${\mathsf{\mu }}_{\mathsf{s}}^{\prime }/{\mathsf{\mu }}_{\mathsf{a}}=$100

We first consider a highly scattering medium probed using a source-detector separation, ${\rho }_{\mathrm{sds}}=3{\ell }^{\star }$. Figures 5(a) to 5(c) display the spatially resolved absorption sensitivity ${\sigma }_{{\mu }_a}(r)$ and depth-resolved cumulative sensitivity ${\mathrm{\Sigma }}_{{\mu }_a}(z)$ distributions for the normally oriented detector, while Figs. 5(d) to 5(f) display these same results for an oblique detector oriented at $\alpha =30\mathrm{deg}$. Purely diffusive coupled forward-adjoint light transport occurs in locations where the computation of the sensitivity distribution is relatively unaffected by the degree of angular discretization used. In Figs. 5(a), 5(b), 5(d), and 5(e), this occurs at depths $z>2{\ell }^{\star }$ and lateral locations $y>{\ell }^{\star }$ for both normal and angled detectors. However, for depths $z<2{\ell }^{\star }$ and lateral locations $y<{\ell }^{\star }$, the use of no or coarse angular discretization under predicts the sensitivity distribution. The use of coarse angular discretization also under predicts the sensitivity distribution in the superficial region between the source and detector as seen in Figs. 5(a), 5(b), 5(d), and 5(e). However, the results of Figs. 5(c) and 5(f) show that the depth-resolved absorption sensitivity provided by coarse discretization agrees well with the calculations obtained using both moderate and fine angular discretizations. Figures 5(c) and 5(f) also display the photon hitting density,¹⁸ which substantially underpredicts the sensitivity for $z\lesssim 3{\ell }^{\star }$ for both normal and oblique detections. These predictions are noticeably poorer than those obtained using no angular discretization. This is the case as the photon hitting density not only ignores any angular dependence in the radiance but also utilizes the standard diffusion approximation to compute the forward and adjoint light fields, which behave poorly near boundaries and collimated sources/detectors. Comparison of Figs. 5(c) and 5(f) reveals that while the angled detector increases the sensitivity to superficial tissue structures, the differences observed between various angular discretization levels are similar. This is due to the high level of scattering, which effectively randomizes the photon propagation directions over small distances.

Fig. 5

Spatially resolved absorption sensitivity distributions are displayed on a ${\mathrm{Log}}_{10}$ scale for highly scattering media ${\mu }_s^{\prime }/{\mu }_a=100$ with ${\rho }_{\mathrm{sds}}=3{\ell }^{\star }$ for $\alpha =0\mathrm{deg}$ (a–c) and $\alpha =30\mathrm{deg}$ (d–f). Distributions are displayed for the $x\text{-}z$ plane at $y=0$ (a, d) and the $y\text{-}z$ plane at $x={\rho }_{\mathrm{sds}}/2$ (b, e) where the colormap provides the results obtained using fine angular discretization. The corresponding depth-dependent cumulative absorption sensitivity is plotted for the four angular discretization cases and the photon hitting density (PHD) (c, f).

Figure 6 displays the sensitivity distributions for source-detector separation, ${\rho }_{\mathrm{sds}}={\ell }^{\star }$. In this case, the forward and adjoint radiances possess much greater directional dependence as compared to the ${\rho }_{\mathrm{sds}}=3{\ell }^{\star }$ case considered above. An examination of Figs. 6(a), 6(b), 6(d), and 6(e) reveal that the computations made with no angular discretization result in significant overestimations of the penetration and underestimation of the magnitude of the absorption sensitivity distribution. The underestimation of the sensitivity distribution at proximal depths is due to strongly direction-dependent radiance coupling at locations in the immediate vicinity of the source/detector. By contrast, sensitivity overestimation at further distances from either the source or detector arises from coupling of one strongly directed light field with the other that has become more directionally randomized. This latter scenario reduces the sensitivity, while the coupling of two isotropic fields produces an overestimation. This is consistent with recent computations that have shown the localized persistence of a strong directional dependence in the radiance distribution at distances up to $3{\ell }^{\star }$ from collimated sources at a depth of ${\ell }^{\star }$ in strongly scattering media.³⁵

Fig. 6

Spatially resolved absorption sensitivity distributions are displayed on a ${\mathrm{Log}}_{10}$ scale for highly scattering media ${\mu }_s^{\prime }/{\mu }_a=100$ with ${\rho }_{\mathrm{sds}}={\ell }^{\star }$ for $\alpha =0\mathrm{deg}$ (a–c) and $\alpha =30\mathrm{deg}$ (d–f). Distributions are displayed for the $x\text{-}z$ plane at $y=0$ (a, d) and the $y\text{-}z$ plane at $x={\rho }_{\mathrm{sds}}/2$ (b, e), where the colormap provides the results obtained using fine angular discretization. The corresponding depth-dependent cumulative absorption sensitivity is plotted for the four angular discretization cases and the PHD (c, f).

The use of coarse discretization provides marked improvements in prediction accuracy that begins to resemble the results obtained using fine angular discretization. This is displayed in the spatially resolved cross sections at depths $z>{\ell }^{\star }$ and lateral positions beyond the source and/or detector. However, Figs. 6(c) and 6(f) show that the coarse angular discretization still overestimates the cumulative absorption sensitivity distribution as a function of depth. Figure 6(f) shows that this overestimation is more prominent in the case of oblique detection. In the case of normal detection, Fig. 6(c) shows that the depth-dependent absorption sensitivity calculated using moderate angular discretization is in strong agreement with results computed using fine discretization. However, Fig. 6(f) shows the use of moderate angular discretization for the case of oblique detection overpredicts the results provided using fine angular discretization for depths $z\lesssim {\ell }^{\star }/2$. The lack of agreement between the absorption sensitivity functions computed using different degrees of angular discretization demonstrates the high degree of angular anisotropy in both the forward and adjoint light fields at depths up to $2{\ell }^{\star }$, when using source-detector separations on the order of ${\ell }^{\star }$. Of particular significance is the misleading picture of the absorption sensitivity obtained when using no angular discretization. These results consistently underestimate the sensitivity to superficial tissue regions while overestimating the sensitivity to deeper ($z>{\ell }^{\star }/2$) tissue volumes. Moreover, the results provided by the photon hitting density for the depth-dependent cumulative absorption sensitivity dramatically underpredict those of cFAMC, and agreement between the diffuse predictions with MC occurs at depths beyond the ROI. These photon hitting density results are consistent with the poor fluence estimates obtained using low-order approximations of light transport resulting from narrow beam irradiation of highly scattering media near the source and boundaries.^36,37

5.2.

Moderately Scattering Media ${\mu }_s^{\prime }/{\mu }_a=$1

When considering moderately scattering media (${\mu }_s^{\prime }/{\mu }_a=1$), we see that the reduction in scattering relative to absorption serves to reduce the angular dispersion of the photon propagation direction and results in increased anisotropy of the forward/adjoint radiances relative to highly scattering media. Figure 7 displays the absorption sensitivity distributions in moderately scattering media for a source-detector separation of ${\rho }_{\mathrm{sds}}=3{\ell }^{\star }$. Examination of Figs. 7(a), 7(b), 7(d), and 7(e) reveals that for both normal and oblique detections, the use of no angular discretization leads to an overprediction of the absorption sensitivity at large depths and lateral positions and an underestimation of the sensitivity in regions proximal to the source/detector as well as in the proximal region between the source and detector. This underestimation at proximal locations and overestimation at distal locations occur via the same mechanism observed in Fig. 6. Specifically, both the forward and adjoint light fields are highly directionally dependent, which leads to an underestimation of the sensitivity distribution in the proximal region of the source/detector and overestimations at distal locations. Due to the reduction in scattering relative to absorption in this case, the overestimations extend over a larger region since directional dependence of the radiance persists over an extended region.³⁵

Fig. 7

Spatially resolved absorption sensitivity distributions are displayed on a ${\mathrm{Log}}_{10}$ scale for moderately scattering media ${\mu }_s^{\prime }/{\mu }_a=1$ with ${\rho }_{\mathrm{sds}}=3{\ell }^{\star }$ for $\alpha =0\mathrm{deg}$ (a–c) and $\alpha =30\mathrm{deg}$ (d–f). Distributions are displayed for the $x\text{-}z$ plane at $y=0$ (a, d) and the $y\text{-}z$ plane at $x={\rho }_{\mathrm{sds}}/2$ (b, e), where the colormap provides the results obtained using fine angular discretization. The corresponding depth-dependent cumulative absorption sensitivity is plotted for the four angular discretization cases and the PHD (c, f).

Figures 7(c) and 7(f) display the depth-dependent absorption sensitivity in the cases of normal and oblique detections, respectively. These results emphasize the need for fine angular discretization to obtain accurate results due to the persistent directional nature of the forward-adjoint light fields in this case. While the use of coarse discretization provides improved results for proximal depths $z\le 1.5{\ell }^{\star }$, the predictions for depths $z>1.5{\ell }^{\star }$ are still overpredicted and are more pronounced when using oblique detection [Fig. 7(f)]. We also show results using the photon hitting density that dramatically underpredicts those of cFAMC at superficial depths. This is not surprising since the standard diffusion approximation, upon which the photon hitting density is based, is not valid for the case of moderately scattering media. Although it appears that the photon hitting density predictions converge with the cFAMC predictions at depths $z>3{\ell }^{\star }$, this is in fact not the case, and the predictions diverge at larger depths. This divergence leads to an overprediction relative to the cFAMC results and is consistent with fluence predictions in moderately scattering media.³⁷ These results also demonstrate the effectiveness of using oblique detection to enhance the sensitivity to superficial tissues regions. In fact, the comparison of Fig. 7(f) with Fig. 7(c) shows that the peak absorption sensitivity in the superficial tissue region predicted using fine angular discretization is enhanced nearly twofold even when using ${\rho }_{\mathrm{sds}}=3{\ell }^{\star }$.

We conclude by examining the absorption sensitivity results for moderately scattering media at a source-detector separation ${\rho }_{\mathrm{sds}}={\ell }^{\star }$ in Fig. 8. Similar to the ${\rho }_{\mathrm{sds}}=3{\ell }^{\star }$ case, Fig. 8 shows that the use of no angular discretization when coupling the forward-adjoint light fields consistently provides incorrect sensitivity distributions. Likewise, the use of coarse angular discretization improves the calculations but remains inaccurate for all locations. Moreover, the resulting depth-dependent sensitivity distributions are inaccurate regardless of detector orientation. Figures 8(c) and 8(f) provide results that display improvements when using a moderate level of angular discretization $(\mu ,\varphi )=(8,16)$ but still overpredicts the results relative to those computed using fine angular discretization at superficial depths $z\lesssim 1.5{\ell }^{\star }$. Moreover, the inaccuracies in using the photon hitting density are exacerbated for this shorter source-detector separation and again did not converge with cFAMC computations even at distal locations beyond the ROI. The differences in the results obtained using moderate and fine angular discretizations in Fig. 8(f) are amplified when using oblique detection as compared to normal detection shown in Fig. 8(c). Furthermore, similar to the ${\rho }_{\mathrm{sds}}=3{\ell }^{\star }$ case, the oblique detector is nearly twice as sensitive to the superficial region of the moderately scattering turbid medium, as seen by comparison of the peak values of the depth-dependent absorption sensitivity distributions.

Fig. 8

Spatially resolved absorption sensitivity distributions are displayed on a ${\mathrm{Log}}_{10}$ scale for moderately scattering media ${\mu }_s^{\prime }/{\mu }_a=1$ with ${\rho }_{\mathrm{sds}}={\ell }^{\star }$ for $\alpha =0\mathrm{deg}$ (a–c) and $\alpha =30\mathrm{deg}$ (d–f). Distributions are displayed for the $x\text{-}z$ plane at $y=0$ (a, d) and the $y\text{-}z$ plane at $x={\rho }_{\mathrm{sds}}/2$ (b, e), where the colormap provides the results obtained using fine angular discretization. The corresponding depth-dependent cumulative absorption sensitivity is plotted for the four angular discretization cases and the PHD (c, f).

5.3.

Relevance of Coupled Forward-Adjoint Monte Carlo and Applications

Using our cFAMC framework, we have provided absorption sensitivity distributions computed using various levels of angular discretization for a homogeneous tissue sample. However, it is important to note that the proposed cFAMC framework can be applied to heterogeneous tissues with arbitrary boundaries. These results display the importance of angular discretization to obtain accurate predictions of coupled forward-adjoint light fields. Thus, accounting for the detailed directional dependence of light transport is critical when using optical diagnostic methods to probe mesoscopic tissue volumes. Such considerations are particularly important for moderately scattering media, that is, when probing tissues using visible wavelengths, as the tissue is less effective in randomizing the direction of light transport. This has the effect of extending the spatial region over which detailed consideration of directional light transport is required. Furthermore, the use of fine angular discretization may be insufficient in certain spatial locations. Future work will examine the use of spatially varying angular discretization schemes and/or adopt other means of angular discretization³⁸ as a means to improve the computational efficiency and accuracy. Moreover, the underlying computational efficiency of the cFAMC methods may be enhanced further by adopting advanced methods for photon propagation and tallying such as next event estimation.^39,40

Our cFAMC framework is also of value to analyze light transport over macroscopic spatial scales, even though the need for angular discretization may be limited to locations near boundaries and collimated sources/detectors. For example, the use of the ROI in which to compute and couple forward-adjoint light fields, as well as our novel method for the tracking and weight deposition of the photon paths within the ROI, provides computational efficiencies regardless of the transport domain under consideration.

The accurate computation of the absorption sensitivity distributions is not only important in the resolution of inverse problems but can also be used to assess the efficacy of different measurement configurations and/or wavelengths to probe localized tissue regions. For example, the determination of the spectral variation of sensitivity functions can identify those wavelengths that give the best contrast in their sensitivity distribution and in turn aid the a priori design of measurements for functional optical imaging. For example, such a tool could provide enhanced optical specificity in high-density optical tomography,⁴¹ while reducing the number of source-detector pairs needed for the tomographic inverse problem.

While consideration of the absorption sensitivity has dominated the field of diffuse optical tomography,¹⁴ the cFAMC approach can also provide a means to efficiently compute the scattering sensitivity and enable an evaluation of a given source-detector configuration to detect structural modifications of the cellular and extra cellular compositions.¹⁴ Furthermore, higher order light transport approximations that provide more accurate assessments than the diffusion approximation have been shown to enhance the tomographic reconstruction of absorptions and scattering⁴² and bioluminescence.⁴³ Thus, the transport rigorous approach that cFAMC offers can provide important improvements for these optical imaging modalities.