2.1.

Normal Vector Estimation and Smoothing

In an effort to better estimate the direction of the mismatched surface, an edge detection algorithm, the Sobel filter,²³ is applied to create a gradient map. By applying such filters to a 3D map with each voxel containing its respective index of refraction, not only can the edges be calculated but also the direction of change is specified. In addition, by applying a Sobel–Feldman type operator, an additional smoothing effect is observed. More sophisticated algorithms to better approximate the normal vectors to the surface certainly exist, but this report is limited to investigating the simple $3\times 3\times 3$ Sobel filter to illustrate the usefulness of this approach and highlight the low computing cost of this additional step. In the 2D case, the Sobel operator can reduce from a maximum angle error of 45 deg using the facet normal of a voxel down to an estimated maximum error value of 1.36 deg.²⁴ The $3\times 3\times 3$ Sobel filter is readily available in MATLAB through the function imgradientxyz.

Figure 1 illustrates the method for specifying the boundaries and the unit vectors that are normal to that surface in each voxel containing the surface. An expanded explanation containing the implementation details is available in Sec. 5. The steps are (following the order of Fig. 1) as follows.

Fig. 1

Proposed algorithm to handle curvature and oblique angles in MC light simulations.

The justification for the smoothing is that when calculating direction vectors, the Sobel filter is local in scope and suffers from the staircase pattern that forms when discretizing into a regular grid. Thus, the smoothing helps the algorithm to estimate curvature more realistically as it gives a dependence of each direction vector on its neighbors.

These normal vectors that are calculated near the edges, where a refractive index mismatch occurs, are then saved into a storage matrix that can be consulted by the MC software. In Fig. 1, to illustrate the concept, the $x$, $y$, and $z$ directions of the normal vectors for each voxel of the simulation domain were saved into three 3D cubes. The implementation of the gradient calculations and smoothing steps was made in MATLAB and the smoothed gradient map was then saved in a binary file that was subsequently read by mcxyzn. It is important to note that this map is only consulted when a photon is about to cross a boundary with a refractive index mismatch, thus, it only alters the performance of the MC software when this normal vector needs to be read from memory. If memory is an issue, it is trivial to create a 3D hash table, a data structure that maps keys to values, that would be associated with a matrix that would only store the direction vectors for the nonzero vectors. Each entry of the hash table would give the index (or key) corresponding to the direction vector in the look-up table. As a photon crosses a boundary into a new voxel with a different refractive index, the hash-table pointer of the new voxel would point to the hash table to find the surface normal vector. The benefit being that the 3D storage of one key value at each voxel uses less memory than having to store a triplet of vector values at each location.

2.2.

Level of Details in Calculating Fresnel Equations

In Fig. 2, we show two methods: the facet normal approach, where only the normal vectors to the facets of voxels are considered, and the surface normal approach, where the normal estimate to the mismatched boundary is being used in the propagation of photons. It is important to note that the facet normal approach can lead to significant errors in the angles of transmittance and reflectance but also can create large discrepancies in the proportion of light that transmits versus reflects at the mismatched boundaries. The surface normal approach reduces these errors significantly.

Fig. 2

Using virtual surface that lies at a 45 deg angle for simplicity, various transmission and reflectance examples are shown for (a) the facet normal approach and (b) the surface normal approach. In the facet normal approach, small variations in the angle of incidence at which a photon hits a mismatched boundary can lead to dramatic changes in the reflected and transmitted angles. Using very similar scenarios under the surface normal approach corrects for these errors in angles.

When a photon is at the mismatched boundary, simply using the normal gradient calculated from the Sobel filter already represents a significant improvement over the facet normal approach. A downside of this approach is that when modeling a curvature in a low-scattering medium (e.g., light “bouncing off” the surface of a spherical object in water), the change in angle direction is still discrete in nature and not gradual as one might expect. To further enhance the proposed algorithm, the interpolated surface normal approach is introduced in the next section to deal with curvature.

2.3.

Intermediate Angle Interpolation Approach

To model gradual changes in normal angles, the use of the trilinear interpolation^27,28 is introduced. The basic method interpolates the values of gradient at the eight voxels surrounding the photon, applying separately to the $x$, $y$, and $z$ gradients. In the modification used in this work, the interpolation only considers the voxel gradients on the side of the surface from which the photon is propagating. The interpolated normal vector is then normalized. The details on how to apply this interpolation are provided in Sec. 5.2.

This interpolated normal approach has a higher computational cost than the simple surface normal approach in the previous section. Instead of having to read the values of one direction vector, the interpolated surface normal algorithm is now required to read the gradient direction values of the eight closest neighbors to this intersection point as well as their respective index of refractivity, and then apply the interpolation scheme to these values.

3.1.

Laser Beam Hitting the Surface of a Sphere

In this first example, Fig. 3(a) shows the expected reflected angles from hitting the side surface of a sphere. The facet normal approach in Fig. 3(b) has reflectance directions going back against the laser beam. Here, the error in angle orientation is as high as 80 deg, but the error in reflected angles can be up to 180 deg. The deficiencies of this approach become evident as one deviates from simulating layered models.

Fig. 3

Example of a beam hitting the side surface of a sphere of radius 1 cm. (a) Photon trajectories calculated by Fresnel law, (b) the facet normal method, (c) the surface normal method, and (d) the interpolated surface normal method. The refractivity indices are $n=1.00$ outside the sphere and $n=1.33$ inside the sphere. The values of ${\mu }_s$, ${\mu }_a$, and $g$ are picked such that there is no scattering and negligible absorption along the photon paths. The voxel length is 0.01 cm.

On the other hand, the reflected angles in Fig. 3(c) (surface normal method) and Fig. 3(d) (interpolated surface normal method) mimic almost identically the analytical trajectories shown in Fig. 3(a). The differences in reflected and transmitted angles between the analytical solution and the facet normal method are drastic. In contrast, the difference between the analytical solution and that of the surface normal method, interpolated or not, is small, with less than 5 deg error when juxtaposing the edges of the reflected beams. When comparing the interpolated surface normal approach to the simpler surface normal approach, we show that the use of interpolation alleviates almost completely the problem of bunching observed in Fig. 3(c). There are very small gaps observed and that are believed to be due to the approximations made by the proposed approach (as opposed to the detailed approach mentioned in Sec. 2). While the angles are accounted appropriately, the intersections of the photons with mismatched boundary surfaces remain on the faces of the boundary voxels. The observed discontinuities are likely resulting from not performing additional calculations to move the intersections onto the “virtual surface” (illustrated in Fig. 2), a small trade-off in accuracy to preserve fast performances during the simulation.

3.2.

Lenses

In this second example, a laser beam hitting a lens in a water medium is considered. Figure 4(b) shows that the algorithm is able to specify the focusing of the incoming light by the lens. The algorithm shows excellent agreement with the analytical solution presented in Fig. 4(a). The surface normal approach, either with or without interpolation, adds the ability to properly focus light, which does not occur when using the facet normal approach. Note that the fluence rate inside the lens is increased due to reflectance from the lens surface when using the interpolated surface normal method. The difference in using the facet normal approach as opposed to the interpolated surface normal approach is quite significant. Differences of over 90% in fluence rate are observed near the bottom and the sides of the lens. Large sections of the fluence map also show differences of about 20%.

Fig. 4

Example of a beam hitting a lens ($n=1.52$) with the lower half being either water ($n=1.33$) or aqueous scattering solution ($n=1.33$, ${\mu }_s=100{\mathrm{cm}}^{-1}$, $g=0.90$) using the interpolated surface normal approach (ISN). (a) Problem setup; (b) lens and water (no scattering) using ISN approach; (c) lens and scattering aqueous solution using ISN approach; (d) lens and scattering aqueous solution using the facet normal (FN) approach; (e) percent change caused by mismatched voxel refractive indices, $\frac{{\phi }_{\mathrm{ISN}}-{\phi }_{\text{matched}}}{{\phi }_{\text{matched}}}\times 100\%$, with scattering present; and (f) percent change caused by using the FN method, as opposed to ISN approach, $\frac{{\phi }_{\mathrm{FN}}-{\phi }_{\mathrm{ISN}}}{{\phi }_{\mathrm{ISN}}}\times 100\%$. The voxel length is 0.01 cm and absorption is assumed to be negligible.

While it is relatively easy to measure the simulation performance (photons/min), the important alterations that this new, more realistic algorithm has on the photon trajectories make it nontrivial to compare to existing approaches such as the facet normal approach or the simplified approach of ignoring the differences in refractive index. Performance values are reported in Table 1 for the lens examples involving scattering. Speeds (photons/min) were found to be comparable between the surface normal approach and no Fresnel approach. A significant performance drop is not expected going from a facet normal approach to a surface normal approach unless the photon trajectories are significantly altered, which can be the case when dealing with curved and oblique surfaces. This is because the only additional operation that the surface normal approach requires is reading a set of three values indicating the surface normal at the current voxel whenever it encounters a mismatched boundary. The interpolated surface normal approach, on the other hand, drops to about 37% of the speed of its noninterpolated counterpart due to the additional global memory reads and calculations required to perform the interpolation.

Table 1

Reported GPU performance for lens simulations.

3.3.

Biological Tissues: Mouse Head

The mouse example shown in Fig. 5 represents a typical example of biological tissue simulation, where a beam hits the air/tissue surface boundary at an angle. Figure 5(e) compares the error in the simulation results between using the interpolated surface normal method and the simplified approach that ignores the differences in refractive index. There is a noteworthy weaker overall diffusive reflectance using the interpolated surface normal method as photons that have entered the tissue media have a harder time escaping due to the internal reflection at the mismatched boundary. Also, photons that have entered the tissue remain there longer, leading to a higher fluence rate near the air/tissue boundary. When comparing the interpolated surface normal approach to the facet normal approach, Fig. 5(f) shows that the facet normal handles poorly the reflected fluence and the fluence transmitting out of the mouse body. There are also significant differences between more than 30% in the areas near the inside boundary of the mouse.

Fig. 5

Example of a beam hitting a mouse head. The shape is taken from the Digimouse segmentation.²⁹ The inside is simplified to be standard tissue ($n=1.33$, ${\mu }_a=1{\mathrm{cm}}^{-1}$, ${\mu }_s=100{\mathrm{cm}}^{-1}$, $g=0.90$). (a) Problem setup, (b) fluence rate without Fresnel calculations, (c) fluence rate using the facet normal approach, (d) fluence rate using the interpolated surface normal approach, (e) percent change in $\varphi$ caused by mismatched voxel refractive indices, and (f) percent change caused by using the facet normal as opposed to interpolated surface normal approach. The voxel length is 0.01 cm.

5.1.

MATLAB Code to Calculate Direction Vectors

%% T is the tissue volume

%% nv is the vector/list of all the refractive indices contained in the tissue list

%% Removing the repeated refractive index values

unique_nv = unique(nv);

%% Creating empty storage maps for the calculated direction vectors

gradient_map_x = zeros(size(T,1),size(T,2),size(T,3));

gradient_map_y = zeros(size(T,1),size(T,2),size(T,3));

gradient_map_z = zeros(size(T,1),size(T,2),size(T,3));

%% Looping through each unique value of refractive index value

for j = 1:length(unique_nv)

 %% Creating 3D binary map for the current unique refractive index value

 n_map = zeros(size(T,1),size(T,2),size(T,3));

 for i = 1:length(n_map(:))

  %% Setting voxel that matches refractive index of the map to 1

  if nv(T(i)) == unique_nv(j)

   n_map(i) = 1;

  end

 end

 %% Finding the direction vectors for the current refractive index value

 [Gx, Gy, Gz] = imgradientxyz(n_map);

 %% Smooth this direction vector map using smoothn

 smooth_map = smoothn({Gx,Gy,Gz},2);

 %% Selectively save the direction vectors contained within the boundaries of each refractive index value domain

 gradient_map_x(n_map == 1) = smooth_map{1}(n_map==1);

 gradient_map_y(n_map == 1) = smooth_map{2}(n_map==1);

 gradient_map_z(n_map == 1) = smooth_map{3}(n_map==1);

end

%% Normalizing final direction vectors

gradient_magn = sqrt(gradient_map_x.^2+gradient_map_y. ^2+gradient_map_z. ^2);

gradient_magn(gradient_magn==0) = 1;

gradient_map_x = gradient_map_x./gradient_magn;

gradient_map_y = gradient_map_y./gradient_magn;

gradient_map_z = gradient_map_z./gradient_magn;

5.2.

Explanation of the Implementation Details

The numbering in this section follows that of Fig. 1.