2.1.

Monte Carlo Software Package

The developed Monte Carlo software (MCtet) uses a tetrahedral mesh to represent the objects, which has several benefits compared with voxels or layered arrangement. Arbitrary shapes can be taken into account adaptively by refining the mesh in regions of high curvature on the one hand and taking large elements for less structured regions on the other hand. This leads to a well-balanced state of accuracy and memory demand. In this program, each tetrahedron is expressed by four integers, linking the four nodes to the matrix, where the coordinates of all nodes are saved. Additionally for each of a tetrahedron’s four faces, the connected neighbor tetrahedra, the normal vectors, and the total reflectance angles are stored to accelerate the simulation speed. Different optical properties can be assigned to each tetrahedron individually by addressing a material number to each element. A flow chart of the MCtet is shown in Fig. 1.

Fig. 1

Flow chart of the Monte Carlo software package. The parts inside the red box handle photons propagating inside the turbid medium and are executed on the GPU. The parts outside the red box treat the photons traveling through air and are executed on the CPU based on the ray-tracing kernel Embree. Both parts are running in parallel. (PF, phase function)

Some functions of the Monte Carlo simulation are similar to those described by Wang et al.,²⁶ like the calculation of the step length, the survival roulette, the decision if a photon is reflected at a refractive index mismatch, and the calculation of a direction for the HG phase function. The functions needed due to the tetrahedral representation of the objects, most importantly the intersection test of a photon with the tetrahedron borders, are comparable with the functions described by Shen and Wang,^27,28 except for the following improvement. Shen and Wang checked all four faces if the ray intersects the plane of the face. By temporarily storing the face number, where a photon entered a tetrahedron, only three faces need to be checked for the subsequent intersection test. All parts of the program concerning the light propagation of a photon inside the diffuse medium by the Monte Carlo algorithm can be found inside the red box in Fig. 1.

The parts inside the red box are implemented with OpenCl.⁴¹ Therefore, the program can be executed in parallel on GPUs and CPUs of different vendors compared with a CUDA-implementation, for which the software would be limited to NVIDIA-GPUs. Nested OpenMP is used to spread the work load across multiple GPUs and to use several CPU cores for the ray tracing per GPU. As the program is implemented with OpenCl, it is also possible to spread the work load of the Monte Carlo algorithm among CPU cores and GPUs at the same time. In most cases, this is not done for two reasons. First, the performance of Monte Carlo simulation for light propagation on a GPU is up to thousandfold faster compared with a CPU.³³ Second, the CPU cores are used in MCtet to handle the light source as described in the following chapter.

2.2.

Simulating Arbitrary Light Sources

In this paper, a concept is presented where neither a preprocessing step is needed nor additional tetrahedra to realize almost any desired light source. Additionally the method enables to track photons leaving a concave object, traveling through air and potentially reenter again the scattering medium at a different position. The elements in the red box of Fig. 1 are the basic Monte Carlo functions running on the GPU, which were expanded for more challenging tasks, such as anisotropic light propagation. Light sources as well as photons leaving the object are dealt by the elements outside the red box based on the fast ray-tracing kernel Embree on the CPU.⁴² Embree is a collection of high-performance ray-tracing kernels developed at Intel. The target users of Embree are graphics application engineers, who want to improve the performance of their applications by leveraging the optimized ray-tracing kernels of Embree.⁴³ The basic idea is to use a highly optimized ray-tracing kernel for photons traveling in air where they do not change their direction and use a parallelized Monte Carlo algorithm for photons propagating in turbid media. Thus, the whole computational power of a system, the GPU as well as the CPU, is leveraged at the same time.

In a first step, the CPU is handling the light source in a way that it computes the first intersection point (IP) of the photon with the object. Therefore, the Embree function rtcIntersect is used to calculate the IP of a ray with an object. A ray is defined by an origin and direction.

There are two situations to define the starting position and direction of a ray. They are either determined by the desired light source or are set to the position and direction of a photon leaving the object. Next the IP with the surface of the object is calculated, if the ray hits the object.

The “scene” is the description of the simulation object by its surface triangles. If a photon is changing from air to turbid medium or vice versa its position and direction needs to be transferred between CPU and GPU. A detailed explanation about the parallelization of the GPU and the CPU can be found in Table 1, where steps 3 + 4 are repeated consecutively until the desired number of photons is completed. A single IP is defined by its position, direction, and the index of the first intersected tetrahedron, resulting in $7\times 4$-bytes of memory and a total memory consumption of only 28 MB, if $N=1\times {10}^6$ IPs are used and the data are stored in single precision. The positions and directions of photons leaving the object can be saved in the same memory as their starting positions and directions, as shown in Table 1, and no additional memory is needed because at most the number of started photons can be reflected, if no photon is absorbed. Furthermore, the communication is run asynchronously, in parallel with the GPU code, at the cost of a higher needed memory $2\times \mathrm{IPs}$ (56 MB). In general, the memory consumption on the GPU as well as the time needed to transfer the data between CPU and GPU cause only minor restrictions for modern graphics cards as they typically have up to 6 GB of memory and a bandwidth of a few GB/s.

Table 1

Parallelization of the CPU and GPU. The brackets […] symbolize the GPU-memory used at the current step (see “memory” in Fig. 1). An IP is defined by the position, the direction, and the index of the tetrahedron where a ray hits the object. A reflected photon (RP) is defined by the current position and direction a photon is leaving the turbid medium. Tasks at the same step are running in parallel and the synchronization is done at the end of each task.

2.3.

Simulation Setup

A human tooth morphology can mainly be separated in three different regions, the enamel, the dentin, and the pulp, as shown in Fig. 2(a).⁴⁴ Enamel is the outer most regions containing $\sim 2.3\%$ of water by weight and the rest mainly consists of the mineral hydroxylapatite.⁴⁵ For scattering in enamel, the HG phase function was applied. The biggest part of the tooth by volume is the dentin located under the enamel. Dentin contains $\sim 10\%$ of water by weight, 70% of hydroxylapatite by weight and the remaining parts are organic materials, mainly collagen and elastin.⁴⁵ The assumed reduced scattering coefficient ${\mu }_s^{\prime }$ for enamel and the absorption coefficient ${\mu }_a$ for enamel and dentin were chosen depending on the wavelength, comparable with ones reported previously.^46–48 The refractive index $n$ was set constant to 1.633 and to 1.54 for enamel and dentin, respectively. The developed software MCtet is able to consider the important anisotropic light-propagation behavior of dentin caused by the tubules. Tubules are small pipes, embedded in the dentin material, with an assumed radius of $1\mu \mathrm{m}$, which are filled with water, and a typical concentration of $c_A=4.5\times {10}^4{\mathrm{mm}}^{-2}$ was chosen.⁴⁹ The tubules form an aligned structure that extends from the dental enamel junction (DEJ) to the pulp as shown in Fig. 2(a) and serves as supply channel for the tooth.⁵⁰ In Fig. 2(b), the defined orientations of the tubules in the simulation object for every 20th dentin-tetrahedron are illustrated. Tubules are the main scatterers in dentin and cause a high-scattering coefficient, such that the scattering coefficient of dentin is significantly higher than the scattering coefficient of enamel, whereas the absorption coefficient is moderately higher. For both materials, the scattering is much higher than the absorption, causing the typical white color of a tooth. The pulp consists of living connective tissue located at the center of a tooth, which is supplied by blood via the tooth root. Typical values for the optical properties of connective tissue were used for the pulp.^51,52 The absorption coefficient of the pulp is varied in different simulations, based on a defined concentration of oxygenated hemoglobin and deoxygenated hemoglobin as the main absorber of blood in this spectral range.⁵³

Fig. 2

(a) Human tooth morphology with its different components and the schematic orientation of the tubules extending from the DEJ to the pulp.¹⁵ (b) Illustration of the tubules’ orientation (blue lines) for every 20th tetrahedron in dentin of the used simulation object. The solid orange area marks the pulp. (c) Cross section of the used tetrahedral tooth mesh and the optical fibers to demonstrate the Monte Carlo simulation setup. The different colors indicate the different materials (dark blue: gingiva, light blue: dentin, green: enamel, orange: pulp, red: bone).

For all of the following simulations, the same incisor tooth model (left maxillary central incisor, tooth 21) was used and a cross section is shown in Fig. 2(c). The tetrahedral tooth mesh was generated at the Department of Orthodontics at Ulm University based on cone beam-computed tomography of the specimen.^54,55 The model consists of $1.9\times {10}^5$ nodes and $9.0\times {10}^5$ tetrahedra. Five different materials were separated for the simulation. The tooth itself is consisting of three materials: enamel, dentin, and the pulp as described before. The whole tooth is grounded in bone, which was assumed to be a highly scattering, low absorbing medium, and has no significant influence on the presented results. On top of the bone, a 1.5-mm-thick gingiva tissue layer is present. For the gingiva, the same optical properties as for the pulp (connective tissue) were assumed except of the absorption, which was fixed to oxygenated hemoglobin with a concentration of $3\mathrm{g}{\mathrm{HbO}}_2/\mathrm{L}$. An overview of the used optical properties is given in Table 2. MCtet was used to calculate transmission spectra through a tooth, considering the anisotropic light propagation caused by the tubules. For each of the investigated cases (different absorptions in the pulp), an analogue simulation has been made to demonstrate the influence of the tubules on the light propagation, where all parameters were kept constant but the scattering in dentin was changed. Each analogue simulation calculates a second transmission spectrum assuming only rotationally symmetrical scattering in dentin, represented by the HG phase function. These simulations are called “HG case,” whereas the simulations considering the anisotropic light propagation caused by the tubules are called “cylinder case.” To estimate the scattering coefficient and anisotropy factor of dentin for the HG case, it is assumed that the tubules would be randomly oriented. In Sec. 2.4, the needed equations are derived to calculate the scattering coefficient and anisotropy factor of dentin. The results are shown in Fig. 3. For the illumination and detection, an optical fiber with an NA of 0.22 and diameter of $365\mu \mathrm{m}$ was modeled. This source type can easily be realized by the concept presented in this paper. The detection fiber was placed lingually and the illumination fiber was placed labially. Both fibers were placed on the tooth as close to the gingiva as possible without overlapping.

Table 2

Assumed optical properties for the different materials of the tooth model from 400 to 700 nm.46–53

Fig. 3

Optical properties of dentin assuming that the tubules are randomly oriented calculated with Eqs. (3) and (4). These properties were used for the HG case.

2.4.

Anisotropic Light Propagation

The anisotropic light propagation behavior caused by aligned, e.g., fibrous, media has been reported for different materials, such as wood or hair.^56,57 In the dental field also anisotropic light propagation appears, caused by tubules in the dentin.^15,49 A detailed description of the human tooth morphology can be found in Sec. 2.3.

To correctly describe the light propagation in a tooth, a scattering phase function obtained by scattering light at a cylinder is used in addition to a rotationally symmetric phase function, in this case the HG phase function, for the Monte Carlo simulation. If a photon in the simulation is scattered by a cylinder at a certain incident angle $\zeta$, its direction is restricted to be on the so-called scattering cone, as shown in Fig. 4(a). For an extreme angle $\zeta \to 90\mathrm{deg}$, all scattered photons directions will be in the same plane, defined by the incident direction and a vector perpendicular to the orientation of the cylinder. For another extreme angle $\zeta \to 0\mathrm{deg}$, the photons effectively do not change their direction as the scattering cone is extremely narrow and is pointing in forward direction. The effect, for $\zeta \to 0\mathrm{deg}$, that photons are scattered in the direction of the cylinders if they travel parallelly to the elongation of the cylinders is the reason for the so-called light guiding effect described in Sec. 3.2. For the Monte Carlo simulation, the probabilities for the scattering angles $p(\lambda ,\zeta ,\theta )$ on the scattering cone in Fig. 4(a), as well as the scattering efficiency $Q_{\mathrm{sca}}(\lambda ,\zeta )$, are needed for all possible incident angles $\zeta$ for the desired wavelength $\lambda$. The scattering efficiency is defined as $Q_{\mathrm{sca}(\lambda ,\zeta )}=\frac{C_{\mathrm{sca}}(\lambda ,\zeta )}{G}$, where $C_{\mathrm{sca}(\lambda ,\zeta )}$ is the scattering cross section and $G$ is the particle cross-sectional area projected onto a plane perpendicular to the incident beam, in this case $G=2r$ where $r$ is the radius of the cylinder.⁵⁸ These values are obtained by solving Maxwell’s equations for an incident electromagnetic plane wave scattered by an infinitely extended cylinder.⁵⁸

Fig. 4

(a) Scattering cone for one incident angle $\zeta$. (b) Cylindrical scattering phase functions $\mathrm{p}(\lambda ,\zeta ,\theta )$ on the scattering cone for $\lambda =500$, 600, and 700 nm, for $\zeta =10\mathrm{deg}$, 40 deg, and 70 deg, as well as the scattering efficiency $Q_{\mathrm{sca}}(\lambda ,\zeta )$ for these wavelengths.

Resulting values are exemplary shown in Fig. 4(b) for three wavelengths for the case of tubules, where a radius $r=1\mu \mathrm{m}$, a refractive index $n_i=1.33$ for the cylinders, and $n_o=1.54$ for the surrounding medium was assumed. Calculations are not done during the simulation but in a preprocessing step, where $p(\lambda ,\zeta ,\theta )$ is tabulated in a cumulative way to efficiently obtain a direction. An example how to draw a scattering angle $\theta$ on the scattering cone with a random number $\chi$ is shown in Fig. 5. In MCtet, the orientation of the cylindrical scatterers as well as the concentration $c_A$ can be defined for each tetrahedron individually. Hence, the light propagation during the simulation is highly depending on a photon’s direction relative to the orientation of the cylinders inside the current tetrahedron. In addition to the scattering phase function, also the scattering coefficient ${\mu }_s(\lambda )$ has to be adjusted.

Fig. 5

Cumulative scattering phase function, calculated for the scattering of a plane wave at an infinitely large cylinder, for $\lambda =700\mathrm{nm}$ and an incident angle $\zeta =40\mathrm{deg}$, to directly draw a scattering angle $\theta$ based on a random number.

The scattering coefficient

At each scattering event, either HG scattering or cylinder scattering is happening. Based on a random number $\chi \in [\mathrm{0,1}]$ scattering at a cylinder happens if

In Sec. 3, the influence of a cylindrical phase function is shown for a tooth model by comparing simulations where the cylindrical scattering is considered with simulations where only HG scattering is assumed. On comparison, the corresponding optical properties for the HG phase function are needed. Therefore, it is assumed that all cylinders are oriented randomly in space, instead of being aligned. Then, a scattering coefficient

3.1.

Validation

To check the correct implementation of the developed Monte Carlo software, comparisons with the software packages MCML, CUDAMCML, MCOnline, MCxyz, and timOS were done.^{26,27,33,59,60} All comparisons were performed on an Intel Core i7-950 (3.07 GHz) with eight CPU cores using a Microsoft Windows 7 system. The used graphics card was an NVIDIA GeForce GTX 480. The simulations for the MCML, MCxyz, and timOS software packages were performed on one CPU core and were supposed to scale linearly with the number of additionally used cores. The calculations done by CUDAMCML were performed exclusively on the GPU. Results with MCOnline were generated using the online user interface at Ref. 61. Simulations with MCtet were performed on two CPU cores and the GPU simultaneously. MCtet is designed to spread the work across multiple GPUs, decreasing the simulation time linearly with the number of used GPUs. For the following tests, the relative errors were computed as the fraction $\frac{x_0-x}{x}$, where $x_0$ is the respective results computed by MCTet and $x$ is the corresponding results computed by CUDAMCML, MCML, MCOnline, MCxyz, and TimOS. Comparisons with these five software have been performed for three different problem sets.

The first case proofs the correct implementation of the feature allowing photons to leave and reenter a concave object. Therefore, the spatially resolved reflectance for a three-layered slab, with a thickness of 4 mm for each layer, was calculated. The optical properties of the first and the last layer are set to ${\mu }_a=0.01{\mathrm{mm}}^{-1}$, ${\mu }_s=5{\mathrm{mm}}^{-1}$, $g=0.8$, and $n=1.4$, whereas the middle layer as well as the surrounding medium were set to be air. Please note that the second layer needs no representation for the MCtet software but must be explicitly defined as transparent for all other programs. To simulate $4\times {10}^8$ photons, MCML needed $3.6\times {10}^4\mathrm{s}$, CUDAMCML needed 27.71 s, and MCtet took 103.44 s, respectively. The spatially resolved reflectance calculated by MCtet as well as the relative errors compared with the other programs are shown in Fig. 6(a).

Fig. 6

(a) Spatially resolved reflectance $I(\rho )$ (red curve) and the relative errors compared with CUDAMCML and MCML for a three layered medium. (b) Spatially resolved fluence $F(\rho )$ (red curve) at a depth $z=2\mathrm{mm}$ and the relative errors compared with CUDAMCML, MCML, MCOnline, and MCxyz for a single slab. (c) Histogram of the relative error of the fluence comparison per tetrahedron with timOS for a cube.

For the second test case, the spatially resolved fluence at a depth $z=2\mathrm{mm}$ for a single slab with a thickness of 4 mm was computed and compared with CUDAMCML, MCML, MCOnline, and MCxyz. The optical properties of the slab are set to ${\mu }_a=0.01{\mathrm{mm}}^{-1}$, ${\mu }_s=5{\mathrm{mm}}^{-1}$, $g=0.8$, and $n=1.0$. To simulate $4\times {10}^8$ photons, MCML needed $9.4\times {10}^3\mathrm{s}$, CUDAMCML needed 83.77 s, MCxyz needed $2.3\times {10}^4\mathrm{s}$, and MCtet took 433.76 s, respectively. For MCOnline, it was not possible to set the number of photons or measure the time needed for the simulation. The spatially resolved fluence calculated by MCtet as well as the relative errors compared with the other programs are shown in Fig. 6(b).

For the first two test cases, a very good agreement with CUDAMCML and MCxyz was found and a deviation when comparing with MCML. Similar deviations with the MCML software have been reported previously by Shen and Wang,²⁷ where they explained it by the bias of the used pseudorandom number generator. The comparison with MCOnline showed a good agreement, but a simulation with more photons would be needed for a better discussion.

For the third test case, the fluence inside a cube ($2\times 10\times 10\mathrm{mm}$) was calculated and compared with the software package timOS. The cube is represented by 6172 nodes and 17954 tetrahedra and has the same optical properties as the layers in the slab described in first test case. The fluence was calculated and compared per tetrahedron. A histogram of the relative errors shown in Fig. 6(c) gives a good agreement of the two simulations.

The comparisons show a correct implementation of MCtet, also for the introduced method to consider photons leaving and reentering a concave object. The GPU-based CUDAMCML program has a four up to five times better performance compared with MCTet due to the restriction to layered media and therefore its high specialization. In particular, CUDAMCML does not have to check for the IP with the tetrahedron borders. Comparisons of the simulation times indicate a superior speed improvement of the GPU-implementation compared with the programs running on CPU.

3.2.

Monte Carlo Simulation Results

Using the simulation setup as described in Sec. 2.3, a change of the oxygen saturation, from 100% to 0%, in the pulp was simulated for different concentrations of hemoglobin in the pulp. The concentrations were assumed to be 0.75 g $\mathrm{Hb}({\mathrm{O}}_2)/\mathrm{L}$, 3 g $\mathrm{Hb}({\mathrm{O}}_2)/\mathrm{L}$, and 7.5 g $\mathrm{Hb}({\mathrm{O}}_2)/\mathrm{L}$, respectively. For a concentration of 150 g $\mathrm{Hb}({\mathrm{O}}_2)/\mathrm{L}$ for whole blood, the used values relate to a volume concentration of 0.5%, 2%, and 5% of blood in the pulp.⁶² It was assumed that the change in concentration only affects the absorption coefficient and not the other optical properties. The transmission spectra for the cylinder and HG case as well as for varying hemoglobin concentrations are shown in Fig. 7.

Fig. 7

Transmission spectra for different scattering types in the dentin as well as different concentrations of oxygenated and deoxygenated hemoglobin in the pulp. The left column shows the transmission spectra for the cylinder case, the right column shows the result for the HG case. For the first row, the hemoglobin concentration was set to 0.75 g $\mathrm{Hb}({\mathrm{O}}_2)/\mathrm{L}$, for the second row to 3 g $\mathrm{Hb}({\mathrm{O}}_2)/\mathrm{L}$, and for the third row to 7.5 g $\mathrm{Hb}({\mathrm{O}}_2)/\mathrm{L}$. The black dotted line for the cylinder case indicates a case where 0 g $\mathrm{Hb}({\mathrm{O}}_2)/\mathrm{L}$ was assumed for the pulp. This spectrum not shown for the HG case as it is identical with the other lines. The total illuminated power was set to 1 W per wavelength.

For the cylinder case, the transmission spectra of the two hemoglobin types are getting more distinguishable for higher concentrations of hemoglobin in the pulp for a spectral range between 500 and 600 nm and less distinguishable for the spectral range between 400 and 450 nm. This finding is caused by the large difference of the hemoglobin absorption of the two spectral ranges. In contrast, for the HG case, the hemoglobin types cannot be distinguished in their transmission spectra. This behavior can be explained by the light-guiding effect of the tubules described by Kienle et al.¹⁵ It was found that light is guided along the orientation of the tubules due to scattering when illuminating them approximately parallel. Although the illumination is pointed on the enamel, which is scattering and therefore blurs the incident direction, a main direction of the light is still conserved when hitting the dentin. This main direction is pointing toward the pulp, approximately parallel to the preferred direction of the tubules.⁵⁰ Hence, the light-guiding effect is still dominant and major portions of the incident light are guided to the pulp. This can be seen by looking at the fluence, showed for a cross section in the incident plane for 570 nm in Fig. 8. A clearly higher intensity on a straight line toward the pulp can be seen for the cylinder case in Fig. 8(a) (top) compared with the HG case in Fig. 8(a) (bottom). The same effect holds for the parts of the light passing through the pulp and are likely to be guided toward the detection fiber, which can be seen in Fig. 8(b) (top). For the HG case, the direction-independent scattering coefficient is very high and no light guiding effect is presented in dentin. Therefore, less photons get to the pulp and reach the detection fiber, causing that the transmission spectra are insensitive to a change of the oxygen saturation. This direction-independent, high scattering is also the reason why the transmitted intensity can partially be higher for the HG case compared with the cylinder case, in spectral regions where the pulp is strongly absorbing. First, the light is not guided toward this highly absorbing part of the tooth. Second, the scattering in enamel is lower but the refractive index is higher compared with dentin. The combination of the low scattering and relative high refractive index mismatch results in another “light guiding” effect in the enamel. Comparable with an optical fiber, light is guided around the tooth inside the enamel to the detection fiber, although the enamel is not completely covering the dentin in the incident plane at the used tooth model [see Fig. 2(c)]. The effect can be seen by looking at a plane located parallelly to the incident plane $z=1\mathrm{mm}$ toward the incisal edge [Fig. 8(b)], where the enamel is starting to completely surround the dentin. By comparing the two fluence plots in Fig. 8(b), it can be seen that the fluence in the enamel is higher for the HG case, because more light is scattered back to the enamel from the dentin compared with the cylinder case. Light that is scattered back to the enamel is then guided inside the enamel around the dentin due to the refractive index mismatch and low scattering. The same effect is observed for the case where cylindrical scattering is assumed but less pronounced because a smaller fraction of the light is scattered back from the dentin into the enamel.

Fig. 8

Fluence compared for cylinder case (top) and HG case (bottom) for 570 nm and an oxygenated hemoglobin concentration of 3 g ${\mathrm{HbO}}_2/\mathrm{L}$ in the pulp. (a) $z=0\mathrm{mm}$ denotes the cross section through the tooth in the illumination-detection plane. (b) $z=1\mathrm{mm}$ and (c) $z=-1\mathrm{mm}$ denotes a shift of the $z=0\mathrm{mm}$ plane toward the tooth incisal edge and the tooth root, respectively. The outlined material is dentin, in the middle the pulp is located, and the surrounding material is enamel.

The transmission spectra in Fig. 7 for the cylinder case show the typical feature of oxygenated hemoglobin, two dips at 542 and 577 nm and for deoxygenated hemoglobin a single dip at 556 nm, respectively. Interestingly for the HG case only, if at all, a double dip can be observed independently on the type of hemoglobin in the pulp. The double dip is getting more pronounced for an increase of the ${\mathrm{HbO}}_2$ as well as an increase of the Hb concentrations in the pulp. In addition to the pulp, the only tissue containing blood is the gingiva. Thus, this evolving double peak must be caused by the gingiva as the absorption and scattering coefficient of the other materials are smooth for this spectral range. This behavior can be explained by the described light guiding effect in the enamel, which is causing a higher fluence in the enamel. Therefore, more light can interact with the gingiva, which is partially covering the enamel [see Fig. 2(c)], and retrieve to the detection fiber. To clarify this, further simulations have been done where the absorber in the gingiva was changed to deoxygenated hemoglobin with the same concentration 3 g Hb/L. In these simulations, the results for the cylinder case remained almost unchanged, where the oxygen saturation of the pulp could still be distinguished in the transmission spectra in terms of a single and double dip. The transmission spectra for the HG case are still nonsensitive to the change of the hemoglobin type but showed the feature of deoxygenated hemoglobin, a single dip at 556 nm, induced by the gingiva.