2.1.

Experimental Setup

The experimental setup used in the measurements and previously reported³ is based on a diode laser operating around $760\mathrm{nm}$ where the molecular oxygen A-band is located. The laser was tuned across a sharp individual line [typically a half width at half maximum (HWHM) of a few gigahertz] by applying a saw-tooth ramp to the driving current of the diode laser. To achieve sensitive lock-in detection, a sine-wave was superimposed to the saw-tooth ramp in a wavelength-modulation scheme.

The light was guided to the sample by a fiber. Two detection geometries were used, as shown in Fig. 2 . In the backscattering geometry, a right-angled prism located in the center of the detector was used to launch the light into the scattering medium. An annular aperture with a central block was used to collect photons backscattered from the sample. In the transmission geometry, the fiber was placed on the opposite side of the sample with respect to the detector. The detection was made by a photomultiplier to provide efficient photon collection. A blocking colored glass filter was placed in front of the detector to reduce ambient light.

Fig. 2

Phantom models used in experiments and simulations representing (a) measurements on the human frontal sinuses with a backscattering detection geometry, and (b) measurements on the human maxillary sinuses with a transmission detection geometry.

The signal from the detector was split into two parts: one directly connected to an oscilloscope, and an other processed in a lock-in amplifier. The peak-to-peak value of the oxygen absorptive signal, recorded as the second derivative, was estimated and normalized to the amount of light reaching the detector. By applying a method called standard addition, a so-called equivalent mean path length, $⟨L_{eq}{⟩}_{exp}$ , could be estimated corresponding to the distance light has to travel in ambient air to obtain the same signal. This is described in Eq. 1, where $c_{air}$ corresponds to the gas concentration in ambient air, $⟨L_{sample}⟩$ to the mean path length in the sample, and $c_{sample}$ to the gas concentration in the sample.

2.2.

Computer Simulations

The Advanced Systems Analysis Program, ${\mathrm{ASAP}}^{\mathrm{TM}}$ , is a ray-tracing software that can be used to simulate light propagation in optical systems. Its implementation is based on the Monte Carlo method and the Henyey-Greenstein approach, which makes it well suitable for simulations in biological tissue and other light scattering materials.⁴ The basic idea of Monte Carlo is to apply a probability view of light propagation, where an individual photon can be absorbed or laterally scattered by a certain angle according to relevant probability laws applying to the situation. After launching a large number of individual photons, the probability distributions are transformed into intensity distributions that can be compared to measured data. Basically, the physical experiment is repeated in the computer, considering a multitude of identical incident photons that later experience different fates.⁵

All simulated models considered in the present work incorporate scattering materials representing the human tissue and an air space representing the sinuses. By studying the rays reaching the detector in the simulations, the equivalent mean path length could be calculated for all models. For every ray, ${\mathrm{ASAP}}^{\mathrm{TM}}$ computes information such as the flux at every surface it enters or is reflected off, as well as the distance it has traveled in all media. We then calculate the equivalent mean path length $⟨L_{eq}{⟩}_{sim}$ by adding up, for all rays, the product of the distance $L$ each individual ray has traveled in the air space and its corresponding flux $\varphi$ when reaching the detector, and dividing it by the total flux reaching the detector. This is described in Eq. 2, where $n$ is the total numbers of rays reaching the detector.

It should be noted that rays reaching the detector that have not traveled in the air gap will have $L=0$ , and only contribute to the simulated equivalent mean path length with a diluting factor. These rays are referred to as short-cut photons.

2.3.

Investigated Sinus Models

The phantom model representing the human sinuses used in the experiments and the simulations contains two scattering slabs made of white Delrin-type plastic, with optical properties similar to human tissue (see Fig. 2). The scatterer closest to the light source is referred to as the primary scatterer $S_1$ , and the other one is referred to as the secondary scatterer $S_2$ . In the case of the transmission geometry, $S_2$ is the scatterer closest to the detector, and in the backscattering geometry it is the one located above $S_1$ . The thicknesses of the scatterers are denoted $l_1$ and $l_2$ . In-between the scatterers, a variable air gap was placed. The annular aperture determining the detection area was simulated by a variable absorptive mask.

To describe imaging of the frontal sinuses, a new model was created for simulations. In the case of imaging of the frontal sinuses, the sinuses were represented by two ellipsoids embedded in a scattering medium. The same backscattering detection geometry shown in Fig. 2a was used.

3.1.

Comparison between Experiments and Simulations

To simulate the phantom experiments described in Ref. 3, which represent measurements on the human sinuses, we created a model in ${\mathrm{ASAP}}^{\mathrm{TM}}$ consisting of two scattering materials $S_1$ and $S_2$ , respectively, with a variable air gap in-between (see Fig. 2). The thicknesses $l_1$ and $l_2$ were varied, consistent with the experiment. The optical properties (the anisotropy factor $g$ , the scattering coefficient ${\mu }_s$ , and the absorption coefficient ${\mu }_a$ ) were measured for one of the Delrin plates used in the experiments by the integrating-sphere technique^{6, 7} to be $g\approx 0.87$ , ${\mu }_s\approx 14{\mathrm{mm}}^{-1}$ , and ${\mu }_a\approx 0{\mathrm{mm}}^{-1}$ at $760\mathrm{nm}$ . These values were used in the simulations. A small absorption ${\mu }_a=0.0005{\mathrm{mm}}^{-1}$ was added to better fit the experimental data. The index of refraction $n$ for the Delrin plates was set to 1.48 and for the air gap to $n=1$ .

Initial simulations, including all reflection and refraction phenomena according to the Fresnel and Snell formulas, were performed. It was noticed that the influence of the reflection was of minor importance for angles below the critical angle in our case. Thus it was suppressed in subsequent simulations, resulting in a manageable computer time. Above the critical angle, the total internal reflection was fully taken into account. The motivation to this approach has its roots in the fact that the solid angle for impinging photons corresponding to total internal reflection is much larger than that of shared reflection and refraction. Further, for photons impinging toward an optically denser medium, reflection has a small contribution as compared to penetration into the medium for the dominant fraction of the solid angle. For example, with two 10-mm-thick Delrin plates spaced by $10\mathrm{mm}$ , the underestimation of the true $L_{eq}$ is about 8% and will, for larger spacings, become smaller (numerical values are provided in the Appendix.

In Fig. 3 the simulated and the experimental $L_{eq}$ are shown for the backscattering detection geometry, representing measurements on the frontal sinuses [model according to Fig. 2a]. The inner detection radius $r_i$ was about $6\mathrm{mm}$ , and the outer $r_o$ about $10.5\mathrm{mm}$ in both simulations and experiments. Each simulated data point consists of launching ${10}^6$ rays into the model.

Fig. 3

The equivalent mean path length $L_{eq}$ as a function of air-gap distance in phantom experiments and simulations for the backscattering detection model shown in Fig. 2a, representing measurements on the human frontal sinuses. The optical properties of the primary and secondary scatterers are in the simulations set to $g=0.87$ , ${\mu }_s=14{\mathrm{mm}}^{-1}$ , and ${\mu }_a=0.0005{\mathrm{mm}}^{-1}$ . (a) Fixed primary scatterer thickness $l_1=3\mathrm{mm}$ . (b) Fixed secondary scatterer thickness $l_2=30\mathrm{mm}$ .

As during the experimental trial,³ it was observed in the simulation that an $S_2$ thicker than $30\mathrm{mm}$ did not influence the outcome. This can be understood because the photons being scattered deeper than $30\mathrm{mm}$ into $S_2$ have a very low probability to be detected. As can be seen in Fig. 3, the agreement between the simulations and the experiments is very good. It can be seen that introducing a thicker secondary scatterer results in $L_{eq}$ increasing. This is because a thicker secondary scatterer results in a higher probability for the photons to be scattered back to the primary scatterer and be detected after passing through the air gap. On the contrary, a thicker primary scatterer results in a decrease of $L_{eq}$ . The photons have now a lower probability of reaching the secondary scatterer, and even if they do, the backscattered photons from the secondary scatterer, are now more likely to be scattered outside the finite size of the detector. Because of the finite size of the detector, the signals first increase with increasing air gap until a certain distance, where they start to decrease. The maximum of the signal should be located at an air gap distance as large as possible to prevent two air gaps from giving rise to the same $L_{eq}$ value, which could be critical in the case of real human measurements. The signal behavior is explained in more detail in Ref. 3.

The deviation between the experimental data and the simulated data can be due to many reasons. In the experiments, the thicknesses of the Delrin plates, the air gap, and the annular detection aperture could not be determined as precisely as they can be in the simulations. During the simulation it was also observed that a small change in one of the optical properties influenced the outcome. It should also be noted that the same optical properties were used for all the different Delrin plates in the simulations. This might not be the case, since the plates used in the experiments did not come from the same batch.

Experimentally it is impossible to explore the influence of multiple passes through the air gap. In contrast, this can easily be done in the simulations. Table 1 shows the number of times the detected photons pass over the air gap, for $l_1=10\mathrm{mm}$ and $l_2=30\mathrm{mm}$ , and for different air-gap distances. It should be noted that only even numbers of passes are possible, since the photons have to travel back and forth to be detected. Similar results were observed for the other thicknesses discussed in Fig. 3.⁸

Table 1

Distribution of the number of passes through the airgap for l1=10mm and l2=30mm and different air-gap distances in the backscattering detection geometry according to Fig. 2a.

It can be seen that signals due to multiple passes can be significant, in contrast to what was intuitively assumed in Ref. 3. This illustrates the power of detailed Monte Carlo simulations, taking into account the actual regime of values of the optical constants.

It can also be seen from Table 1 that a larger air-gap distance results in fewer multiple passes. This phenomenon can be explained by the finite size of the detector. It was observed from the simulations that a thicker $S_1$ resulted in a increase of multiple passes. This can be understood because a thicker primary scatterer can collect oblique rays to be detected by the finite size detector that have been through more multiple passes. It was also observed that a thicker secondary scatterer increases the number of multiple passes, as a thicker secondary scatterer allows more photons to be scattered back to the primary scatterer.⁸

In Fig. 4 the simulated and the experimental $L_{eq}$ are shown for the transmission detection geometry, representing measurements on the maxillary sinuses [model according to Fig. 2b]. The mask radius was about $2.5\mathrm{mm}$ in both simulations and experiments. Each simulated data point consists of launching ${10}^6$ rays into the model.

Fig. 4

The equivalent mean path length $L_{eq}$ as a function of air-gap distance in phantom experiments and simulations for the transmission detection model shown in Fig. 2b, representing measurements on the human maxillary sinuses. The optical properties of the primary and secondary scatterers in the simulations are in the simulations set to $g=0.87$ , ${\mu }_s=14{\mathrm{mm}}^{-1}$ , and ${\mu }_a=0.0005{\mathrm{mm}}^{-1}$ . (a) Fixed primary scatterer thickness $l_1=10\mathrm{mm}$ . (b) Fixed secondary scatterer thickness $l_2=10\mathrm{mm}$ .

Again, very good agreement between the experimental and simulated data can be seen. The obtained $L_{eq}$ for the different scatterer parameters increases with increasing air-gap distance. This can be understood, since in the case of the transmission detection geometry, all photons detected must have traveled over the air gap. This is in contrast to the backscattering detection geometry, where the short-cut photons (photons that have only been scattered in $S_1$ before being detected) dilute the observed $L_{eq}$ . The signal behavior is explained in more detail in Ref. 3.

It can be seen that the photon statistics are not as good as for the backscattering detection geometry. For the transmission geometry, the percentage of detected photons in relation to the number of incident photons is much lower than for the backscattering case, despite the fact that all photons are carrying information of the sinuses. A much longer simulation time is needed to achieve good photon statistics. Again, the deviation between the experimental and the simulated data can be due to many reasons, as discussed before.

The number of multiple passes was also investigated in the transmission detection geometry for the models representing measurements on the human maxillary sinuses. Table 2 shows the number of times the detected photons pass the air gap, for $l_1=10\mathrm{mm}$ and $l_2=10\mathrm{mm}$ , and for different air distance. In this geometry, only odd number of passes are possible. Similar results were observed for the other thicknesses discussed in Fig. 4.⁸

Table 2

Distribution of the number of passes through the air gap for l1=10mm and l2=10mm and different air-gap distances in the transmission detection geometry according to Fig. 2b.

Again, as in the case of backscattering detection geometry, it can be seen that signals due to multiple passes can be significant. It can also be seen that for a larger air-gap distance, fewer multiple passes occur. As discussed before, this phenomenon can be understood from the finite size of the detector. It was also observed that a thicker $S_2$ results in more detected multiple pass rays. This can be understood, since more oblique rays will be collected by the finite size of the detector. The inverse argument is applicable to explain why a thicker $S_1$ results in less multiple passes.⁸

3.2.

Influence of the Scattering Coefficient

The scattering coefficient and the thickness of the facial bone may vary between different humans, especially with age.⁹ The scattering coefficient of bone is significantly smaller for children than for adults. In our model, a lower scattering coefficient has similar influence as a smaller thickness of the primary scatterer, i.e., a lower scattering coefficient of $S_1$ results in a higher probability for the photons to reach $S_2$ be scattered back and be detected by the finite size of the detector.

Figure 5 shows the influence of the scattering coefficient of $S_1$ and $S_2$ . The simulations were performed for the backscattering detection geometry shown in Fig. 2a, with an annular detection aperture of dimensions $r_i=6\mathrm{mm}$ and $r_o=10.5\mathrm{mm}$ . The optical properties used in the simulations are taken from different sources and presented in Table 3 .^{10, 11, 12} In all simulations, except for Delrin, ${\mu }_a\approx 0{\mathrm{mm}}^{-1}$ was used. For each trial with a set of optical properties, simulations were performed with the same air-gap intervals as in Figs. 3 and 4, and by launching ${10}^6$ rays into the model for each air-gap distance.

Fig. 5

Simulated equivalent mean path length $L_{eq}$ as a function of air-gap distance, using backscattering detection geometry according to Fig. 2a, with a primary scatterer thickness of $10\mathrm{mm}$ and a secondary scatterer thickness of $30\mathrm{mm}$ . An annular detection aperture with an inner radius $r_i$ of $6\mathrm{mm}$ and another radius $r_o$ of $10.5\mathrm{mm}$ , and optical properties of primary and secondary scatterer according to Table 3 were used in the simulations.

Table 3

Optical scattering properties used in the simulations taking from different sources.10, 11, 12

From the simulations it was observed that only the reduced scattering coefficients ${\mu }_s^{\prime }$ and ${\mu }_a$ influence the result. The difference between optical properties corresponding to adult skull in vivo and neonate skull in vivo is only the reduced scattering coefficient. By comparing the curves obtained for these cases in Fig. 5, it can be seen that a lower ${\mu }_s^{\prime }$ results in a larger air-gap distance needed to obtain the maximum signal. This is because a decrease in ${\mu }_s^{\prime }$ can also be seen as a decrease in the thickness of the primary scatterer. This phenomenon can also be seen in Fig. 3a. The influence of the absorption coefficient can be seen by comparing the simulated curves obtained for Delrin and adult skull in vitro. A lower ${\mu }_a$ gives rise to a higher integrated signal.

3.3.

Annular Detection Aperture

In the experiments, the choice of annular detection aperture was done intuitively. Here we investigate the influence of the aperture geometry using the simulations.

In Fig. 6 , $L_{eq}$ is plotted as a function of the air-gap distance for different annular detection apertures in the case of the backscattering detection geometry. The thicknesses of $S_1$ and $S_2$ were set to 5 and $30\mathrm{mm}$ , respectively. The simulations were done with neonatal optical scattering properties (shown in Table 3). Again, the simulations were performed with the same air-gap intervals as in Figs. 3 and 4, and by launching ${10}^6$ rays into the model for each air-gap distance.

Fig. 6

Simulated equivalent mean path length $L_{eq}$ as a function of the air-gap distance for different annular detection apertures in the case of backscattering detection geometry according to Fig. 2a. The simulations were performed for a model with a primary scatterer thickness of $5\mathrm{mm}$ and a secondary scatterer thickness of $30\mathrm{mm}$ , and with optical properties corresponding to neonate (Table 3). (a) The inner radius of the annular detection aperture $r_i$ is kept constant at $4\mathrm{mm}$ , while the outer radius $r_o$ is varied between 5 and $10\mathrm{mm}$ . (b) The outer radius of the annular detection aperture $r_o$ is kept constant at $10.5\mathrm{mm}$ , while the inner radius $r_i$ is varied between 4 and $9\mathrm{mm}$ .

In Fig. 6a, it can be seen that increasing $r_o$ while keeping $r_i$ constant increases the signal. This is because a larger $r_o$ allows more photons traveling over the air gap to reach the detector. This also gives rise to a larger photon flux reaching the detector. In Fig. 6b, $r_o$ is kept constant while $r_i$ is varied. It can clearly be seen that increasing $r_i$ results in a higher $L_{eq}$ . This is equivalent to blocking out more of the short-cut photons that have only been scattered in the first medium, and thus only influence $L_{eq}$ by diluting it. As can be seen in the figure, an increase in $r_i$ also pushes the maximum of the signal toward larger air-gap distances. This might be very important for preventing two different air-gap distances giving rise to the same $L_{eq}$ value.

The conclusion from these different curves obtained with different annular detection apertures is that an outer detection aperture as large as possible should be used. In the current setup, a photomultiplier (Hamamatsu 5070A) with a radius of $10.5\mathrm{mm}$ is used, which limits the outer radius. The simulations also show that an inner radius as large as possible should be used with an upper limit set by the experimental signal-to-noise ratio. Increasing the power of the diode laser used should thus be advantageous. Similar results were obtained and the same conclusion was drawn for scattering properties according to an adult, and are discussed in detail in Ref. 8.

A study of the influence of the detection aperture in the transmission detection geometry was also performed. In this detection geometry, the size of the aperture did not have the same influence on the signal as in the backscattering detection geometry. As mentioned before, all photons reaching the detector have traveled over the air gap distance, and short-cut photons need not to be blocked. This is discussed in more detail in Ref. 8.

3.4.

Frontal Imaging

To investigate the feasibility of imaging the human frontal sinuses, a model according to Fig. 7a was studied. The two frontal sinuses were represented by two air-filled ellipsoids ( $25\mathrm{mm}$ high, $20\mathrm{mm}$ wide, $10\mathrm{mm}$ deep) separated by $5\mathrm{mm}$ and embedded at a depth of $10\mathrm{mm}$ in a scattering material, with optical properties corresponding to adult in vitro (Table 3). The backscattering detection geometry used in the simulations is shown in Fig. 2a. Figure 7b shows the results of the simulations when using an annular aperture with $r_i=5\mathrm{mm}$ and $r_o=10.5\mathrm{mm}$ . The whole image consists of 1000 simulated data points with the light and detection arrangement scanned across the model. For each position, $0.5\times {10}^6$ rays were launched into the sample. The locations of the ellipsoids are also indicated. The figure clearly shows the possibility of spatially resolving the cavities.

Fig. 7

(a) Model used in ${\mathrm{ASAP}}^{\mathrm{TM}}$ to simulate imaging of the frontal sinuses. (b) Simulated equivalent mean path length $L_{eq}$ obtained when scanning over the model using backscattering detection geometry according to Fig. 2a with an annular detection aperture with an inner radius $r_i$ of $5\mathrm{mm}$ and another radius $r_o$ of $10.5\mathrm{mm}$ . Optical properties according to adult in vitro in Table 3 were used.

Many factors might be considered to determine what constitutes a “good” image. It might be important to maintain the largest possible contrast between maximum and minimum signals, e.g., corresponding to the centers of the ellipsoids (the sinuses), and the center gap between them. To be able to resolve the two sinuses, it might also be critical to use a sufficiently small detection aperture to not blur the image.

In Fig. 8 , the cross sections of the simulated backscattered image are shown for different detection apertures. In Fig. 8a, $r_i$ was kept constant at $4\mathrm{mm}$ , while $r_o$ was varied between 5 and $10\mathrm{mm}$ . The real cross section is also included for the sake of clarity. It can be seen that selecting a larger outer detection radius results in a higher $L_{eq}$ value. In Fig. 8b, $r_o$ was kept constant at $10.5\mathrm{mm}$ , while $r_i$ was varied between 4 and $9\mathrm{mm}$ . It can be seen that the spatial resolution is rather independent of the chosen detection radii. Then it is important to obtain a sufficiently large $L_{eq}$ value while maintaining a large photon flux.

Fig. 8

Simulation of frontal sinus imaging using backscattering detection geometry according to Fig. 2a with different apertures. (a) The inner radius of the annular detection aperture $r_i$ is kept constant at $4\mathrm{mm}$ , while the outer radius $r_o$ is varied between 5 and $10\mathrm{mm}$ . (b) The outer radius of the annular detection aperture $r_o$ is kept constant at $10.5\mathrm{mm}$ , while the inner radius $r_i$ is varied between 4 and $9\mathrm{mm}$ .

Table 4

Numerical values for the phantom simulations when employing the backscattering detection geometry (l1=l2=d=10mm) . The different cases are explained in the text. For each case, 0.5∙106 rays were launched into the sample. In case 3, the criterion that each ray was allowed no more than two divisions was set.