1. Introduction

The ever-increasing usefulness of lasers in medical applications requires the knowledge of optical properties of biological tissues. Many applications require values for the absorption coefficient, $\mu _a$, and the reduced scattering coefficient, $\mu _s^\prime$. These coefficients cannot be measured directly, but instead are derived from experimental measurements of reflectance and transmittance paired with a predictive model such as diffusion theory, adding-doubling, and Monte Carlo photon transport [1–4].

A common technique for measuring the reflectance and transmittance of a sample uses integrating spheres. A photodetector attached to an integrating sphere records signals from a series of measurements, and simple relations between the signals calculate the reflectance and transmittance [4–6]. Figure 1 shows how integrating spheres are used to make these measurements.

 

Fig. 1. Depictions of how integrating spheres would be used in optical experiments measuring (a) reflectance and (b) transmittance as well as (c) a dual-sphere setup measuring both reflectance and transmittance.

Download Full Size | PDF

Monte Carlo photon transport is a common predictive model to calculate reflectance and transmittance of tissue samples based on their optical coefficients. Determining the optical coefficients from reflectance and transmittance measurements is an inverse problem, solved by a lookup table or multivariable minimization.

The reflectance and transmittance depend on the optical properties of the tissue sample as well as the geometry of the measurement configuration [7–9]. Therefore, integrating spheres affect the measured reflectance and transmittance. This is because the presence of an integrating sphere introduces the possibility that a photon exits the sample, reflects within the integrating sphere, and is then re-incident on the sample. These interactions cause some photons which should have been measured as reflectance to interact with the sample again and instead be measured as transmittance—and similarly, initially transmitted photons may instead be measured as reflected. Traditional Monte Carlo photon transport methods estimate the total reflectance and transmittance but do not account for the geometry of the measurement configuration, specifically that of integrating spheres. For the accurate prediction of the measured reflectance and transmittance, and therefore accurate solutions to the inverse problem of determining the optical properties of the sample from these measurements, the simulation geometry must correspond to the measurement geometry. Figure 2 shows the modeled measured total reflectance (specular plus diffuse) and total transmittance (unscattered plus diffuse) for various configurations of integrating spheres. This figure demonstrates that values can often differ by $20\%$ or more depending on the measurement configuration and sample optical properties, and in general the presence of an integrating sphere reduces the measured transmittance.

 

Fig. 2. Simulated effect of integrating spheres—with properties listed in Table 4—on the measured reflectance and transmittance of samples (a) with $\mu _a=1\, \textrm{cm}^{-1}$ and varying $\mu _s^\prime$ and (b) with varying $\mu _a$ and $\mu _s^\prime=10\, \textrm{cm}^{-1}$. All samples had refractive index $n=1.4$, anisotropy factor $g=0.5$, and thickness $t=0.135\, \textrm{cm}$.

Download Full Size | PDF

In this work, we introduce an extension to existing Monte Carlo photon transport methods to simulate the propagation of photons in an integrating sphere—discussed in detail later. This approach requires significantly less computation than ray-tracing photon paths inside a sphere as well as existing theory to correct for the change in measured reflectance and transmittance when integrating spheres are present [10]. This existing theory requires the calculation of the reflectance and transmittance for multiple illuminations, whereas our method only requires a single simulation with the desired illumination to find the reflectance and transmittance when integrating spheres are present. We have implemented this approach, titled Integrating Sphere Monte Carlo (ISMC), in standard Python $3$ for ease of use and CUDA C++ for quick computation times, which can be found in Code 1 [11,12]. We have implemented our method in CUDA because—while not strictly required—Monte Carlo photon transport codes benefit greatly from parallel implementations and CUDA offers significantly more parallelization than traditional multithreading or multiprocessing techniques [11,13].

Propagating the experimental uncertainties from integrating sphere measurements and statistical uncertainties from the Monte Carlo method requires additional Monte Carlo simulations on the inverse solvers. We are able to propagate all uncertainties through our method to calculate a final uncertainty.

This work is primarily concerned with baffled integrating spheres and samples illuminated directionally at normal incidence by well-collimated lasers. However, integrating spheres need not have baffles and samples can be illuminated diffusely. Our model can be generalized to simulate such cases and modifications to accommodate those cases are discussed later.

2. Integrating sphere measurements

Integrating sphere measurements determine two important values: the measured reflectance, $M_R$, and the measured transmittance, $M_T$. Calculating these values with integrating spheres has been discussed in great detail before, so we merely present an overview of single-sphere measurements of samples irradiated directionally, normal to their surface, with relevant equations [4–6]. The arguments to each measurement value of $R$ and $T$ denote the exact type of illumination and sample being measured. This nomenclature follows the notation set out in Prahl (2011) [4]. Here, Prahl defines $R(p,q)$ to be the power measured by the detector on the sphere for a sample with reflectance, $p$, for normally incident light, and reflectance, $q$, for diffuse incident light. The notation is the same for $T(p,q)$; $p$ and $q$ still represent values of reflectance. This nomenclature was intended to be as close to independent of incident angle as possible.

Calculating $M_R$ requires one value and three measurements:

 

Fig. 3. Configurations necessary to make the measurements required for the calculation of $M_R$.

Download Full Size | PDF

Calculating $M_T$ is much the same, except without a standard. The required measurements are:

 

Fig. 4. Configurations necessary to make the measurements required for the calculation of $M_T$.

Download Full Size | PDF

In addition to $M_R$ and $M_T$, the regular, or unscattered, transmittance, $M_U$, can also be measured. The measurements necessary to calculate $M_U$ do not require any integrating spheres, only an aperture and a standalone photodiode—or some other form of photon-measuring tool such as a power meter. Three measurements determine $M_U$:

In practice, $M_U$ is difficult to determine accurately [4]. This is due in part to the massive difference in $U_{100}$ and $U_0$, which makes it possible to saturate the detector or have a poor signal-to-noise ratio during each measurement, respectively. The use of an optical filter in the measurement of $U_{100}$ can help remedy this; however, uncertainty in the optical density of the filter can introduce error in the final measurement of $M_U$ as well.

 

Fig. 5. Configurations necessary to make the measurements required for the calculation of $M_U$.

Download Full Size | PDF

At normal incidence, the regular, or specular, reflectance cannot be measured since it coincides with the incoming beam. Instead, the beam may be incident $8^\circ$ relative to the sample normal in order to include the specular component of the reflection in the sphere, in which case $M_R$ would measure the regular plus diffuse reflectance. This incident angle is the standard for such experiments [14].

Simulation of integrating spheres requires a few additional values: the inner wall reflectance, diameter of the sphere, and diameter of each port. The supplier of the sphere or additional experiments provide these values [4,6]. The simulation of the regular transmittance experiment also requires the maximum permissible polar angle measured from the central axis of the incident laser. This is the angle threshold for regular transmittance: ISMC scores any photon that transmits through the sample with a polar angle less than or equal to this angle as regular transmission. This angle, denoted as $\theta _U$, is calculated from the aperture diameter, $A$, and the distance of the aperture from the sample, $l$, using

3. Method

In the typical Monte Carlo photon transport method, the sample properties, beam properties, and number of photons are defined prior to simulation. In our methodology, we additionally require the inclusion or exclusion of each sphere and the relevant properties of any included sphere. These properties, which can vary between the two spheres, are:

As mentioned before, regular reflectance is often included in the reflectance measurements using integrating spheres. To simulate this case, we define our beam to be incident on the sample at $8^\circ$ relative to the sample normal—which is the standard for such experiments [14]. Otherwise, we define our beam to be incident normal to the sample.

With these properties, photon propagation begins by using the Monte Carlo photon transport method with non-weighted tissue boundary interactions [15]. We use a non-weighted approach for these interactions because in weighted models, fractions of photon packets can leave the sample while the rest of the packet continues to propagate in the sample [16]. With our integrating sphere model, that would require independently propagating an increasing number of different fractions of a single photon packet if part of the packet entered the integrating sphere and was re-incident on the sample. However, weighted absorption can still be used to increase the computational efficiency of propagation within the sample [16]. The methods of propagation in the sample, scoring of absorbed photons, and internal reflection of photons are unchanged from the original Monte Carlo photon transport method [15,16]. Our integrating sphere method is employed when a photon is either reflected, backscattered, or transmitted from the sample. When this happens, we first update the photon’s direction as it leaves the sample by Snell’s Law, ignoring any surface roughness. Then, we check if there is an integrating sphere placed on that side of the sample. If not, ISMC immediately scores the photon as either reflectance or transmittance and then terminates it. If there is a sphere present, we first determine if the photon will leave directly through the opposite port. To do this, we check to see if the direction cosine of the photon in the $\hat {z}$ direction, defined as $\mu _z=\cos {\theta }$ where $\theta$ is the polar angle of the photon’s direction, is within a given range. This amounts to

Once ISMC scores and terminates the photon, another photon launches. This process repeats for as many photons as desired. The uncertainty in each measurement of photon counts follows ${u(M) \propto \frac {1}{\sqrt {N}}}$, where $N$ is the total number of photons in the simulation. This gives a good estimate for the number of photons required for a desired uncertainty. Figure 6 shows a flowchart overview of the method.

 

Fig. 6. A flowchart of our integrating sphere method.

Download Full Size | PDF

It is important to note that this method accurately accounts for crosstalk in dual-sphere experiments. Photons are free to interact with one of the spheres, be re-incident on the sample, propagate through the sample, and finally interact with the opposite sphere.

3.1 Derivation

While integrating sphere theory is well-established [17,18], one detail critical to this method is missing: the probability that a photon in the integrating sphere is re-incident on the sample.

A photon that enters an integrating sphere must end its interaction with the sphere—after any number of reflections from the sphere wall—by being absorbed into the sphere wall, re-incident on the sample, incident on the detector, or by escaping through any open ports on the sphere. Each of these events necessarily has finite probability and these probabilities must sum to unity.

If we denote the reflectance of the sphere wall as $\rho _w$, and enforce that the wall transmits no light, then the probability that a photon is absorbed into the sphere on any given strike is $1-\rho _w$.

For an integrating sphere with a nearly ideal Lambertian coating, the flux between any two differential surface elements is independent of their locations on the sphere surface [19]. It follows that the flux incident on a differential surface element originates uniformly from all other differential surface elements of the integrating sphere. Therefore, the fraction of reflected flux exiting through a port is equal to the fractional area of that port relative to the entire integrating sphere. In terms of individual photons, the probability of exiting the port is this fractional area.

With these interaction probabilities established, we find the probability of re-incidence after a photon enters the integrating sphere. First, the photon must not escape through the source/optional port, which has a probability of $(1-f_p)$. Then, the photon must not be absorbed into the sphere wall, which has a probability of $\rho _w$. Finally, the photon must enter the sample port, which has a probability of $f_s$. Therefore, the probability that the photon is re-incident on the sample after just one reflection is $\rho _w(1-f_p)f_s$. This is the case for a photon entering the sphere in an unknown random direction. In Monte Carlo photon transport, however, the exit angle of the photon is known, so the simulation of the sphere can be omitted if it is known that the photon will immediately exit through the source/optional port. Thus, in practice, the probability that the photon will be re-incident on the sample after one reflection is actually $\rho _w f_s$. The lack of the dependence on $f_d$ of the first reflection comes from the presence of the baffle, which prevents photons from reaching the detector immediately after exiting the sample. It directly follows that if we denote the total fractional port and detector area as $f$ then the probability of the photon being re-incident on the sample after $n$ reflections is $\rho _w^n(1-f)^{n-1}f_s$. Since we want the total probability that the photon is re-incident on the sample, we sum these probabilities over all possible values for the number of reflections to find the total probability, $P$, that the photon is re-incident:

As stated before, we know the direction of photons leaving the sample and want to know if they are on a trajectory that takes them immediately out through the source/optional port. We start with the equation for the surface area of a spherical cap, $A_{cap}$. Given a sphere of radius $r$ and a cap with polar angle $\alpha$, which is the angle from the line segment connecting the peak of the cap to the center of the sphere to the line segment connecting the rim of the cap to the center of the sphere, the surface area is [20]

To derive the angle of re-incidence, we assume that the integrating sphere’s inner wall is an ideal Lambertian surface. One of the properties of integrating spheres is that the flux incident on the wall is independent of viewing angle [19]. This means that the probability of a photon having polar angle $\theta$ and azimuthal angle $\phi$ is

The estimate for the number of photons needed to achieve a given uncertainty follows from counting statistics. The “Square-Root Rule for Counting Experiments” gives us that the uncertainty in any discrete random counting variable is just the square root of the number counted [21]. In the Monte Carlo method, the number of photons being absorbed, reflected, or transmitted are the counted events. For example, the number of photons reflected, $r$, has an uncertainty of $u(r)=\sqrt {r}$. However, we are interested in the reflectance, $R=r/N$ where $N$ is the number of photons. Propagating this uncertainty in $r$ gives us $\frac {r \pm \sqrt {r}}{N} = R \pm \sqrt {\frac {R}{N}}$. Since the absorptance, reflectance, and transmittance are all inclusively bounded between $0$ and $1$, we know the uncertainty in each of these must be bounded above by $\frac {1}{\sqrt {N}}$.

3.2 Assumptions and modifications

We have made a few assumptions in the derivation of our model. We provide an overview of those assumptions and possible modifications to the model here.

We have ignored the effect of surface roughness during any refraction the photon may experience. It is not known exactly how critical this simplification is; however, this is a common assumption in Monte Carlo photon propagation [1,2,15,16,22]. Moreover, during integrating sphere measurements, the sample is fitted between two glass slides, making its surface somewhat smooth.

In the derivation of Eqs. (7), (9), and (10), we have assumed that the baffle is large enough to block all photons leaving the sample port from immediately reaching the detector. Additionally, in our derivation of Eq. (6), we have assumed that the baffle is small enough not to block photons leaving the sample port from escaping through the source/optional port. Commercially available baffled integrating spheres appear to follow this convention, allowing for this assumption. Within these limits, the size of the baffle does not affect any of the aforementioned equations, as the baffle is treated the same as the walls in integrating sphere theory [17,18].

However, an integrating sphere without a baffle allows some of the light leaving the sample to be incident on the detector without ever interacting with the sphere. This light is lost and should immediately be scored as diffuse reflectance or transmittance, depending on which sphere it is in. This contribution can be treated in exactly the same way as photons which leave directly through the source/optional port without interacting with the sphere, as discussed above. Both the polar and azimuthal angles of the photon should be checked as it exits the sample—and enters the sphere—to see if the photon is on a trajectory which will take it directly to the detector. The probability of re-incidence, Eqs. (7) and 10, need not be modified, since this angle threshold check eliminates photons which do not interact with the sphere wall.

Much of integrating sphere theory, including our derivations above, assume that the inner walls of the sphere are an ideal Lambertian surface [17–19]. Many commercial coatings exist today which exhibit nearly ideal Lambertian reflectance, validating this assumption.

In our derivation of Eq. (6) we have assumed that photons exit the sample, and enter the sphere, at the center of the sample port. More realistically, this derivation could have included the photon’s radial distance from the center, which would be known for every photon in the Monte Carlo photon propagation step. Furthermore, photons that exit the sample far enough from the center of the port such that they never enter the sphere could have been accounted for as well. However, for sufficiently thin samples, or samples with a sufficiently low reduced scattering coefficient, this distance would be negligible. Samples used are often very thin, or have low reduced scattering coefficients, so that enough light transmits to be measurable.

Similarly, we have placed all re-incident photons at the center of the sample—which is also the center of the sample port. Equation (8) could be modified to sample $x$ and $y$ uniformly over the surface of the sample that is exposed to the port. This simplification may have the effect of slightly increasing the predicted reflectance and transmittance due to photons that would usually travel far enough in the radial direction to escape the sides of the sample instead having their position reset at the origin, allowing them to be incorrectly scored as reflected, transmitted, or absorbed.

A related effect is that of finite samples. As in the original Monte Carlo photon transport method, we have assumed that the sample is effectively infinite in the $x$ and $y$ directions [15]. However, in reality, if the reduced scattering coefficient is large enough, photons may escape from the sides of the sample and be lost—rather than be counted as reflected or transmitted. This may cause any inverse solver using this method to overestimate the absorption coefficient [23]. Modifying the photon propagator to account for finite sample size may help remedy this issue.

As stated previously, our model can simulate samples illuminated both normally and at $8^\circ$ relative to normal—to include regular reflectance. However, there is another method of illumination commonly used with integrating spheres: diffuse illumination. This entails the laser being incident on the sphere wall, instead of the sample [5,6]. To simulate this case, each photon would need to interact with the integrating sphere on the reflection side of the sample first, meaning that photon propagation would instead start at the “Backscattered or Reflected” path in the flowchart shown in Fig. 6.

4. Verification

4.1 Verification without integrating spheres

To verify our implementation of the base Monte Carlo photon transport method, we compared results with those previously obtained from analytical models and other Monte Carlo simulations [2,16,22,24–26].

The first verification was for a sample with a refractive-index-matched boundary. The total reflectance, $R_t$, and total transmittance, $T_t$, were calculated for a sample with the following optical properties: refractive index $n=1$, absorption coefficient $\mu _a=10.0\, \textrm{cm}^{-1}$, reduced scattering coefficient $\mu _s^\prime=22.5\, \textrm{cm}^{-1}$, anisotropy factor $g=0.75$, and thickness $t=0.02\, \textrm{cm}$. One run of $500,000$ photons was done without any integrating spheres present. The “Square-Root Rule for Counting Experiments” calculated the uncertainty [21]. These results and those from previous analytical models and Monte Carlo simulations are presented in Table 1 [2,16,22,24,26].

Table 1. Verification for a sample with a refractive-index-matched boundary.

View Table | View all tables in this article

The second verification of the base Monte Carlo photon transport implementation was for a sample with a refractive-index-mismatched boundary. The sample had the following optical properties: refractive index $n=1.5$, absorption coefficient $\mu _a=10.0\, \textrm{cm}^{-1}$, reduced scattering coefficient $\mu _s^\prime=90\, \textrm{cm}^{-1}$, anisotropy factor $g=0.0$, and infinite thickness. The diffuse reflectance, $R_d$, was found from one run of $10,000$ photons without any integrating spheres present. The same method as before calculated the uncertainty of the diffuse reflectance. Our results are presented alongside results from previous analytical models and Monte Carlo simulations in Table 2 [2,16,22,25,26].

Table 2. Verification for a sample with a refractive-index-mismatched boundary.

View Table | View all tables in this article

Our implementation of the base Monte Carlo photon transport method is consistent with previous works as demonstrated by Tables 1 and 2. Moreover, our uncertainty calculation method produces uncertainties comparable to those of previous works—while using the same total number of photons and without the need to run multiple simulations.

4.2 Verification with integrating spheres

Reflectance and transmittance measurements of a sample with known optical properties verified the integrating sphere method. We used a $2.5 \, \textrm{cm} \times 2.5 \, \textrm{cm}$ sample of a synthetic epidermis tissue phantom ($2\,\textrm{N}$ caucasian skin, SynDaver, USA) and a collimated $532\,\textrm{nm}$ laser with a beam radius of $3\, \textrm{mm}$. The sample was illuminated normal to its surface. Throughout the experiment, no damage to the sample was observed. The optical properties of the sample are listed in Table 3.

Table 3. Optical Properties of our tissue phantom at $\lambda =532\,\textrm{nm}$ according to a National Institute of Standards and Technology (NIST) report [27].

View Table | View all tables in this article

A single integrating sphere ($4$P-GPS-$033$-SL, Labsphere, USA) with a diameter of $3.3$ inches, three $1$ inch ports, and one $1.5$ inch port measured the diffuse reflectance, $R_d$, of the sample. A plug made of the same diffuse reflecting material as the inner sphere wall capped one of the three $1$ inch ports. The three open ports were the sample port, source/optional port, and detector port. The detector port was always a $1$ inch port, while the sample port and source/optional port alternated between $1$ inch and $1.5$ inch depending on the orientation of the sphere. The same integrating sphere measured the diffuse transmittance, $T_d$. The sphere was reoriented when measuring the diffuse transmittance to ensure the baffle was in the correct position. During the measurement of $T(r_s^{direct},r_s)$, the optional port was left open. The properties of the sphere relevant to our method in each orientation are shown in Table 4. The sphere was removed and an aperture with radius $3\, \textrm{mm}$ was placed $78.3\, \textrm{cm}$ behind the sample to measure the regular transmittance, $T_u$. This corresponds to measuring transmitted photons with a polar exit angle of at most $0.22^\circ$. The regular reflectance was not measured since it coincided with the incident beam. The raw measurement values are presented in Table 5. Our results alongside simulation results from one run of $100,000$ photons are presented in Table 6.

Table 4. Properties of the integrating sphere used to verify the model in each orientation.

View Table | View all tables in this article

Table 5. Raw measurement values for the integrating sphere measurements of our tissue phantom.

View Table | View all tables in this article

Table 6. Reflectance and transmittance values found by experiment and simulation.

View Table | View all tables in this article

Table 6 shows that our integrating sphere method produces results similar to experimental results for diffuse reflectance and diffuse transmittance. We suspect better calibration of experimental materials, as well as proper uncertainty estimation for the nominal reflectance of the reflectance standard used would bring results closer to agreement. However, our method does not agree with the measured regular transmittance. This discrepancy is likely due to the fact that the regular transmittance is the most difficult quantity to measure experimentally, for reasons discussed earlier [4].

5. Application to measuring optical properties

We present two methods to solve the inverse problem of predicting the absorption coefficient and reduced scattering coefficient of a sample: lookup tables and multivariable minimization. Each method has its own weaknesses and advantages. This section provides a general outline of both methods.

For both methods, it is important to pay close attention to how the laboratory measurements were made. Single-sphere setups require two separate simulations by ISMC: one with only a reflection sphere present and one with only a transmission sphere present. Moreover, the regular transmittance of a sample is usually measured without any spheres present, so this measurement must also be simulated separately.

5.1 Lookup table method

To create a lookup table, we first must know the range of allowed values for both the absorption coefficient and reduced scattering coefficient. Within these ranges, we can simulate many samples to obtain values for the reflectance and transmittance. Interpolating these reflectance and transmittance values effectively represents all samples in the given ranges of $\mu _a$ and $\mu _s^\prime$. Finally, inverting the data creates a lookup table that matches a given pair of reflectance and transmittance values to a unique pair of $\mu _a$ and $\mu _s^\prime$.

We devised a fictional sample to demonstrate the lookup table method. The sample had a refractive index $n=1.4$, anisotropy factor $g=0.9$, thickness $t=1\, \textrm{cm}$, and was infinite in the transverse directions. The range of possible absorption coefficients for the sample was from $0.01\, \textrm{cm}^{-1}$ to $2\, \textrm{cm}^{-1}$ inclusive. The range of possible reduced scattering coefficients was from $0.01\, \textrm{cm}^{-1}$ to $44\, \textrm{cm}^{-1}$ inclusive. A zero-sphere setup was simulated. These optical properties, sphere configuration, and in particular large thickness were chosen to create a fictional setup that almost entirely fills the lookup table. The absence of any spheres and infinite transverse size of this fictional sample allow the assumptions made in the model to be valid. Figures 7(a) and 7(b) show the lookup tables for the absorption coefficient and reduced scattering coefficient of this fictional sample. These example tables required $11 {\min }\, 33.337\,\textrm{s }$ of real time to generate using the CUDA implementation of the method on a single desktop computer with an Intel Xeon W-2123 CPU running at $3.6\, \textrm{GHz }$ with $32~\mathrm {GB}$ of RAM and two Nvidia Titan Xp GPUs running at $1.582\, \textrm{GHz }$. This is roughly $72$ times faster than the Python implementation of the method running on the same desktop, which required $834\,{\min }\,40.840\, \textrm{s }$ of real time to generate these tables.

 

Fig. 7. Lookup tables for the (a) absorption coefficient and (b) reduced scattering coefficient for a sample with $n=1.4$, $g=0.9$, and $t=1\, \textrm{cm}$ measured with a zero-sphere setup. Marked points represent values obtained directly from simulation that were interpolated from.

Download Full Size | PDF

This method is particularly useful for visualizing results. Moreover, the evaluation of several samples can be done at once, so long as they share a refractive index, anisotropy, and thickness. However, many simulations are required before any evaluations; in the above example, $900$ simulations. This many simulations can take a considerable amount of time if the implementation is not well-parallelized.

5.2 Minimization method

An alternative to creating a lookup table is the multivariable minimization of an objective function. Since the output of ISMC is approximately smooth and continuous, as evident by Fig. 7, several common minimization algorithms can be used. We define an objective function, denoted Obj $(\mu _a,\mu _s^\prime)$, which uses ISMC to calculate a theoretical diffuse reflectance, $R_d$; diffuse transmittance, $T_d$; and regular transmittance, $T_u$, and then returns the Euclidean distance between these and their corresponding experimentally measured values—denoted by $M_R$, $M_T$, and $M_U$. By minimizing Obj $(\mu _a,\mu _s^\prime)$, we arrive at an absorption coefficient and reduced scattering coefficient that closely match the coefficients of the experimental sample.

Our implementation of this minimization technique uses Powell’s method [28]. This method does not require the input function’s derivative—the function does not even need to be differientiable [29]. However, Powell’s method does require an initial guess, which we denote as $(\mu _a^0, \mu _s^\prime 0)$. In practice, this guess can be very far from the actual solution, but the minimization converges faster the closer this guess is to the solution. A pseudocode outline of our minimization technique is shown in Algorithm 1. ISMC’s Python $3$ implementation uses SciPy’s Powell’s method, while ISMC’s CUDA C++ implementation uses PRAXIS, a refinement of Powell’s method implemented in C++ [29–31].

This technique benefits greatly from the added measurement of the regular transmittance of the sample. Furthermore, for coefficients of just one sample—or several samples that do not share a refractive index, anisotropy, or thickness—this method is much faster. The $cost$ value in Obj can also account for the regular reflectance, if that is measured. The anisotropy, $g$, can also be solved for using this method, since Powell’s method can find the minimum of a function of any number of variables; Obj would just need to be modified to take $g$ as an input.

6. Propagation of uncertainties

All laboratory measurements have uncertainties, including integrating sphere measurements. Monte Carlo methods are by definition stochastic and therefore have some intrinsic uncertainty in their results. Here we present methods to propagate all of these uncertainties through to final uncertainties in the predicted optical coefficients.

Ideally, the uncertainty in all simulated values would be as small as possible. However, since the time and space complexity of the Monte Carlo photon transport method increases much faster with larger numbers of photons than the $\frac {1}{\sqrt {N}}$ decrease of the uncertainty, a compromise between sufficiently precise results, reasonable computation times, and memory requirements must be found. We have chosen values of $N$ which result in $\approx 1\%$ uncertainty in the simulated values of $R$ and $T$.

6.1 Integrating sphere measurements

Each of the values in Eqs. (1), 2, 3, and 4 has some uncertainty which must be propagated through to find the total uncertainty in the measured diffuse reflectance, measured diffuse transmittance, measured regular transmittance, and angle threshold for regular transmittance. The total uncertainty in $M_R$, $M_T$, $M_U$, and $\theta _U$ are found using the “law of propagation of uncertainty” [32].

6.2 Lookup table method

Propagating uncertainties through a lookup table requires the forward propagation of simulation parameter uncertainties, followed by the backward propagation of integrating sphere measurement uncertainties.

The forward propagation must propagate the uncertainties in the refractive index, inner wall reflectance of both spheres, and port fractions of all ports in both spheres as well as statistical uncertainties associated with the Monte Carlo approach. Several simulations accomplish this, each with different input parameters drawn from Gaussian distributions for each parameter. Each distribution has a standard deviation equal to that of the uncertainty in that parameter. Sampling two more Gaussian distributions—one for the interpolated reflectance values and one for the interpolated transmittance values—determines uncertainties in interpolated lookup table values. The standard deviations of these distributions are the uncertainties reported by the simulation. The final lookup table with an associated uncertainty table, from input parameters and Monte Carlo uncertainties, is thus given by the mean and standard deviation respectively of these simulations. During this process, it is possible for a value given by the distributions for reflectance or transmittance to be outside the evaluated range of the lookup table. When this happens, the uncertainty is set to be undefined, rather than be incorrect. The consequence of this is that the table that represents the uncertainty from simulation alone will often appear smaller than the original lookup table, since the missing uncertainty values almost always fall on the edge of the original table.

Backward propagation of integrating sphere measurement uncertainties occurs after the creation of the lookup table and simulation uncertainty lookup table. Sampling two more Gaussian distributions for $M_R$ and $M_T$, with standard deviations equal to the associated uncertainties, creates several lookups for each optical property. The means of these lookups are the final predicted optical coefficients. The standard deviations of these lookups are the uncertainty in each optical coefficient from experimental uncertainties alone. The final uncertainty adds this uncertainty to the uncertainty from simulation alone in quadrature. A flowchart showing both the forward and backward propagation of uncertainties through the lookup table method is shown in Fig. 8.

 

Fig. 8. A flowchart showing the propagation of uncertainties through the lookup table method.

Download Full Size | PDF

6.3 Minimization method

Propagating uncertainties through the minimization method involves only the forward propagation of all simulation and experimental uncertainties. This is because the integrating sphere measurements of diffuse reflectance, diffuse transmittance, and regular transmittance are all inputs to this method. Also unlike the lookup table method, the uncertainties in $M_U$ and $\theta _U$ are propagated.

As before, we run several minimizations with different input parameters sampled from Gaussian distributions. The mean of these minimizations is the final reported result and the standard error of the mean is the final uncertainty. Figure 9 shows the propagation of these uncertainties.

 

Fig. 9. A flowchart showing the propagation of uncertainties through the minimization method.

Download Full Size | PDF

7. Example inverse calculation

Using the values for the integrating sphere measurements made on the tissue phantom, listed in Table 6, we present the predicted absorption coefficient and reduced scattering coefficient using both inverse solvers. Table 7 shows the input parameters given to each solver with associated uncertainties. Since tolerances for the integrating sphere used for the experimental measurements could not be found, we estimated the uncertainty of each port fraction to be $5\%$. This is purposely an overestimatation in order to mitigate unreasonably precise results. The results from both methods alongside the exact values for both optical coefficients of the sample are presented in Table 8. The lookup tables, as well as the uncertainties from these tables from the forward propagation only, are shown in Fig. 10. These tables required $11\,{\min }\,6.521\, \textrm{s }$ of real time to generate using the CUDA implementation of the method running on the same desktop mentioned before. This is nearly $600$ times faster than the Python implementation which ran on the same desktop and required $6568\, {\min }\,33.212\, \textrm{s }$ of real time.

Our results from both methods are consistent with one another and with values provided by a NIST calibration report, which was consistent with the Inverse Adding Doubling (IAD) method [4,27]. Both of our methods produced comparable uncertainties to each other and to the NIST calibration report.

Table 7. Input parameters and associated uncertainties used to solve for the optical coefficients of our tissue phantom as well as the method each parameter is used in.

View Table | View all tables in this article

Table 8. Inversely calculated absorption and reduced scattering coefficients of our tissue phantom at $\lambda = 532\,\textrm{nm}$ from integrating sphere measurements using both the lookup table method, the minimization method, and a NIST calibration report [27].

View Table | View all tables in this article

 

Fig. 10. Lookup tables for the (a) absorption coefficient and (c) reduced scattering coefficient for our tissue phantom measured with a single-sphere setup as well as the percent uncertainty from simulation parameters only in the (b) absorption coefficient and (d) reduced scattering coefficient Marked points represent values obtained directly from simulation that were interpolated from.

Download Full Size | PDF

8. Conclusion

We presented a Monte Carlo method for modeling integrating sphere experiments that produces results consistent with laboratory measurements. We additionally demonstrated two methods for solving the inverse problem of determining the optical coefficients of a sample from integrating sphere measurements. Uncertainties from simulation and experimental measurements are propagated through both methods, and we demonstrated that both methods produce accurate and precise results. This model and both inverse solvers are implemented in standard Python $3$ and CUDA C++ and are publicly available on GitHub for use by anyone who may want to determine the optical properties of samples from integrating sphere measurements [11,12].

9. Future work

A logical next step in this work is verifying dual-sphere simulations. Moreover, since this method does not require any modification to the base method of Monte Carlo photon transport in tissue, a more sophisticated model such as Monte Carlo Multi-Layered (MCML) could replace the current propagator [22]. Many of the modifications to the model discussed previously, such as modeling diffuse illumination, finite samples, and more realistic re-incident positions, may help improve both the accuracy and precision of our model. Implementing a weighted re-incidence model, not unlike the weighted absorption model found in many Monte Carlo photon propagation models, could improve the computational efficiency of our method [16,22]. Combining the CUDA and Python implementations of ISMC with PyCUDA in the future would increase both the usability and accessibility of the code. Recent work has presented importance sampling inversion algorithms for goniometric measurements which are significantly more efficient than traditional inversion techniques [33]. Implementing such an approach in ISMC for integrating sphere measurements would drastically decrease the computational time required.