Page updated: April 14, 2020
Author: Curtis Mobley
View PDF

Photon Tracking

The previous page showed how to determine photon path lengths and scattering angles using the beam attenuation c, the scattering phase function β̃, and a uniform U[0,1] random number generator. This page shows how to combine those two processes to create random photon paths through an absorbing and scattering medium.

Photon Tracking

There is more than one way to simulate photon paths, and each will give the same answer. However, some techniques can be numerically much more efficient than others. Indeed, a reasonable approach to developing a Monte Carlo algorithm for a particular problem is to

first figure out how to numerically simulate a process as it occurs in nature, and
then figure out how to simulate another, perhaps artificial, process that will give the same answer as the ”natural” process, but with less computational time.

This page illustrates this two-step development process.

Consider first how photons propagate through a medium. Loosely speaking, a photon travels until it interacts with a particle, e.g. a molecule of water or chlorophyll. It is then either absorbed by the particle and disappears, or it is scattered into a new direction and continues on its way until in interacts with another particle.

Recall the albedo of single scattering, ωo = bc. If there is no absorption, b = c and ωo = 1. If there is no scattering, b = 0 and ωo = 0. ωo = bc thus can be interpreted as the probability of photon survival in any particular interaction. When a photon encounters a particle, we can randomly decide if the photon is to be absorbed or scattered as follows:

Draw a random number 𝔯 from a U[0,1] distribution.
Compare 𝔯 with ωo.
  • If 𝔯 ωo, then the photon is scattered.
  • If 𝔯 > ωo, then the photon is absorbed.

If the photon is absorbed, tracing stops and a new photon is emitted from the source and tracing begins anew. If the photon is scattered, two new U[0,1] random numbers are drawn and used to determine new polar and azimuthal scattering directions ψ and χ as shown on the previous page. Another random number is used along with c to determine the distance traveled before another interaction.

Figure 1 illustrates this process for two photons, which also introduces the geometry to be used in numerical simulations below. A source emits a collimated beam of photons, which are then recorded by an annular, ring, or ”bullseye” detector some distance away. The red photon is emitted by the source, undergoes one scattering, and is then absorbed by a particle. The green photon is emitted, undergoes two scatterings, and is recorded by a detector.


Figure 1: 1. Illustration of type 1 photon tracing. The red photon is absorbed and the green one reaches the target.

This process mimics what happens in nature. We will call this ”type 1” photon tracing (there are no standard names for ways of tracing photons). Note that all of the computations used to trace the absorbed photon are wasted because the photon never reached the detector. Nature can afford to trace innumerable photons and waste some by absorption, but that is not advisable for most numerical simulations. We therefore seek other ways to trace photons.

The previous page showed that the mean free path or average distance between interactions with the medium is 1c. These interactions can result in either absorption or scattering of the photon, as just described. Rather than tracing one photon at a time as nature does, consider a source emitting ”packets” of many photons. Then view each interaction as having a fraction 1 ωo of the photons in the packet being absorbed, and the remaining fraction ωo being scattered, all in the same direction. Let the packet be emitted with an initial weight of w = 1, which can represent one unit of energy, power, or some number of photons. At each interaction, the current weight w is multiplied by ωo to account for the loss of energy or number of photons by absorption (that is, a fraction ωo continues onward). The scattered packet then carries a reduced weight. If the photon packet reaches the target, the current weight w is tallied. Another packet is then emitted from the source and traced. This tracing process, which we’ll call ”type 2,” is illustrated by the green photon track in Fig. 2. After two scatterings, as illustrated, the detected photon packet has weight w = ωo2.


Figure 2: 2. Illustration of type 2 photon tracing in green and type 3 in red.

A third ray-tracing process can be envisioned. The mean distance traveled between scattering events is 1b. We can thus use 1b and a random number to determine the distance between scattering interactions, and the initial weight of w = 1 is not changed at each interaction because all photons in the packet are viewed as being scattered. Then, if the photon packet reaches the target, absorption is treated as a continuous process occurring along the entire photon path. Assuming homogeneous water, the final weight tallied is then e(a), where is the total path length in meters and a is the absorption coefficient. The red track in Fig. 2 illustrates this ”type 3” ray tracing. The red track shows a total path length of 1 + 2 + 3, so the final weight is e[(1+2+3)a].

The two tracks in Fig. 2 are drawn as though each track were generated by exactly the same sequence of random numbers. Because 1b > 1c, the individual type 3 photon paths will be greater than the type 2 paths. The scattering angles are the same. Thus these two tracing types clearly lead to different results, photon packet by photon packet. However, numerical simulation of many photon packets shows all three of these photon tracing types yields the same distribution of energy at the detector.

To summarize, the three types of ray tracing considered here are

Type 1:
Individual photons are tracked, and photons can be absorbed.
Type 2:
Photon packets are tracked, with a packet weight being multipled by ωo at each interaction.
Type 3:
Photon packets are tracked, with track lengths determined by the mean free path for scattering, no weighting at scattering events, and absorption treated as a continuous process based on total path length.

Numerical Comparison of Tracking Types

To illustrate the results obtained for different ways of tracking photons, a Monte Carlo code was written to simulate the energy received by an annular target as illustrated in Figs. 1 and 2. For the simulations shown here, the IOPs were defined by a = 0.2m1, b = 0.8m1, and a Fournier-Forand phase function with parameter values (n,μ) = (index of refraction, slope of Junge distribution) chosen to give a good fit to the Petzold average particle phase function. Thus ωo = 0.8, and optical distance τ is numerically the same as geometrical distance in meters. A run was made with 106 photons being sent from the source and using type 1 ray tracing. Photons were traced until they were absorbed. Figure 3 shows some of the resulting statistics.

The red histogram shows the percent of photons that traveled an optical distance τ1 τ τ2 between interactions, for a bin size of τ2 τ1 = 1. The theoretically expected fraction of photons traveling an optical distance between τ1 and τ2 between interactions is

τ1τ2 eτdτ = eτ1 eτ2 (1)

The red dots in the figure are the expected values given by this formula. The shortest photon path length between interactions was τ = 1.192 × 107 and the longest was 16.69. The blue histogram shows the distribution of total distances traveled until the photons were absorbed. Thus the value for the first bar shows that about 18% of the photons were absorbed after going a total optical distance between 0 and 1. Note that this distance can represent more than one interaction, i.e., a photon being scattered one or more times before being absorbed. Both the red and the blue histograms sum to 100%. As shown on the previous page, the mean distance traveled between interactions is τ = 1, or 1c in meters. For this particular simulation the actual average was τ = 0.9974 (or 0.9974 m for these IOPs). The small difference is statistical noise resulting from the finite sample size of the numerical simulation. Likewise, the average distance traveled until the photons are absorbed is 1a. For the present case of a = 0.2m1 this gives 5 m. The average for this simulation was 4.9949 m. Since c = 1m1 for this simulation, another way to view this is that the photons were on average scattered four times before being absorbed on the fifth interaction.


Figure 3: Example distributions of the optical distances between interactions and of the total distances traveled before absorption for type 1 ray tracing.

We next compare results for the three different ways of tracking photons. Because oceanic phase functions scatter much more light at small scattering angles than at large angles, most photons that are scattered just a few times will hit the detector near its center. To even out the numbers of photons (or power) detected by each ring, an annular target was defined with a logarithmic spacing for the radii of the detector rings. A logarithmic spacing is often used in instruments so that each detector ring receives roughly equal amounts of power, which reduces the dynamic range needed for the instrument design. The detector simulated here had Nrings = 10 rings with the smallest ring radius being rmin = 0.1 and the largest being rmax = 10. This detector is placed in a target plane some distance zT from the source and centered on the optical axis of the source, as shown in Figs. 1 and 2. The photons crossing the detector plane at some distance r from the detector center are tallied in bins as follows:

Bin 0:
Unscattered photons that hit the detector at r = 0.
Bin 1:
Scattered photons that hit the target plane inside the first detector ring, i.e. at 0 < r < rmin.
Bins 2,...,Nrings + 1:
Photons that hit detector rings 1 to Nrings.
Bin Nrings + 2:
Photons that hit the detector plane outside the outer ring, i.e. at r > rmax.

Simulations were made with Nemit = 106 photons (for type 1) or photon packets (types 2 and 3) emitted from the collimated source. Figure 4 shows the distribution of photons or packets received anywhere in the detector plane as a function of the number of scatterings, for the three types of tracing and for target plane distances of zT = 5 and 15. The left panel shows that for zT = 5 one or two percent of photons (depending on the tracing type) reach the target plane without being scattered. Most photons are scattered 3 or 4 times, and very few photons are scattered more than 10 times. The right panel shows that for zT = 15, almost no photons reach the target plane unscattered, and most undergo 5 to 25 scatterings, with a peak around 10 or 15, depending on the way the photons are traced. For type 1 photon tracing, almost no photons are scattered more than 30 times. Note that the probablity of surviving 30 scatterings is ωo30 = 0.0012. For tracing types 2 and 3, which never have photon packets absorbed, there are broad tails in the number of scatterings, although for type 2 a packet being scattered 40 times has an almost negligible weight of w = ωo40 = 0.00013.


Figure 4: Distribution of photons reaching the target plane as a function of the number of scatterings, for the three scattering types and two target plane distances.

The left panel of Fig. 5 shows the intersection points of the photons reaching a τ = 10 × 10 area of the target plane for the first 104 emitted photons, tracing type 1, and the target plane at zT = 5. Note that of the 104 emitted photons, only 2973 reached the target plane (of which 2966 are in the area plotted). This is less than 30% of the photons making a contribution to the answer of how much power is detected; i.e., 70% of the calculations were wasted. The photon-intersection dots are color-coded to show the number of times each plotted point was scattered. The two black circles show off-axis angles of 30 and 60. Very few photons were scattered more than about 30 off of the optical axis. The right-hand panel shows the distributions of the photons reaching the target plane by number of scatterings and total distance traveled. Every photon must of course travel a distance of at least τ = 5 before reaching the target plane, and most photons were scattered several times.


Figure 5: Distribution of photons reaching the target plane at zT = 5 for 104 emitted photons and type 1 tracing. The left panel shows the spatial distribution of points where the photons intersected the target plane. The two black circles are drawn at 30 and 60 angles off of the optical axis. The right panel shows the distributions of number of scatterings and total distance traveled.

Figure 6 shows the same distributions for type 2 scattering. Note than now 93% of the emitted photon packets eventually intersect the target plane. Only 7% of the emitted photons were wasted. Those photons ended up being scattered into directions away from the target plane (either by backscattering or by multiple large-angle forward scatterings). Figure 7 shows the results for type 3 tracing. The distributions are similar to those for type 2, but with over 94% of the photons reaching the target plane.


Figure 6: Same as Fig. 5 but for type 2 tracing.


Figure 7: Same as Fig. 5 but for type 3 tracing.

Figures 8-10 show the corresponding results when the target plane is at zT = 15. Now, for type 1 tracing, only 179 of 104 emitted photons ever reached the target plane. Over 98% of the photon-tracing calculations were wasted! Note also that there is obvious statistical noise in the distributions of the right panel, due to the small number of photons used to computed the statistics. For types 2 and 3 about 71% and 82%, respectively, of the emitted photons eventually reach the target plane. The statistical noise is now much smaller (but still noticeable) because of the larger number of photon packets.

The total optical distance distributions for types 2 and 3 show broad tails. In both cases, fewer than one fourth of the photons made it to the target plane after traveling a total distance of τ = 1516. About 30% of the photons underwent 30 or more scatterings and traveled a distance of τ 30. These broad tails illustrate the phenomenon of pulse stretching in time-dependent problems. If we think of all Nemit photons being emitted simultaneously, then the longer distances traveled correspond directly to later arrival times at the target. Pulse stretching is an important limiting factor in time-dependent applications such as lidar bathymetry or communications with high-frequency light pulses.


Figure 8: Same as Fig. 5 but for a target plane at zT = 15. The black ring in the left panel is 30 off of the optical axis.


Figure 9: Same as Fig. 8 but for type 2 tracing.


Figure 10: Same as Fig. 8 but for type 3 tracing.

Figure 11 shows the distributions of numbers of photons at the target plane and the corresponding power (or energy) for the three tracing types and the detector at zT = 5. The left panel shows the distributions as a function of the detector ring radii, and the right panel is the same information as a function of the bin number defined above. Recall that the first abscissa point is unscattered photons, the second is scattered photons inside the first detector ring, and the last plotted point is for photons outside the last detector ring. The solid-line histogram represents the 10 detector rings for rmin = 0.1 < r < rmax = 10.

There are several things to notice in this figure. First, the distributions of the numbers of photons (open circles) are different for the three tracing types. For type 1, the distribution of the number of photons is the same as the power distribution because the photons all retain their initial weight of w = 1. Thus power detected is simply the number of photons detected. For types 2 and 3, more photons are detected, but each is weighted less to account for absorption along the way. Finally—and most importantly—the power distributions for these three tracing types are identical to within a small amount of statistical noise, which is not visible in these plots.


Figure 11: Distributions of number of photons (open circles) and power (histograms) for the target plane at zT = 5, for the three types of tracing. The left panel shows the distributions as a function of radius from the optical axis. The right panel displays the same information as a function of the bin number.

Figure 12 shows the power distributions for the detector at zT = 15 but only Nemit = 104 photons emitted. There are obvious differences in the distributions for the three tracing types. However, if Nemit = 106 photons are traced, as in Fig. 13 these differences almost disappear. This indicates that the three ways to trace photons all give identical predictions of the detected power, to within some amount of statistical noise, which can be reduced by tracing more photons.

However, the computation time required by the three tracing types can vary greatly. Recall from Fig. 5 that about 30% of the photons reached the target plane at zT = 5 for type 1 tracing, but that over 90% reached the target plane for types 2 and 3 (Figs. 6 and 7). Thus, if we require a certain number of detected photons to achieve some desired level of statistical noise, we would have to emit and trace over three times as many photons for type 1 tracing (hence three times the computer time) as for type 2 or 3. For the detector at zT = 15 and type 1 tracing, fewer than 2% of the emitted photons reached the target plane (Fig. 8), whereas about 80% of the photons reached the target plane for types 2 and 3 (Figs. 9 and 10). Getting the same number of photons on target would thus require emitting over 40 times as many photons (hence 40 times the computation time) for type 1 as for types 2 or 3.


Figure 12: Energy distributions for zT = 15 and Nemit = 104, for the three tracing types. There are obvious differences in the distributions.


Figure 13: Energy distributions for zT = 15 and Nemit = 106, for the three tracing types. The differences seen in Fig. 12 have almost disappeared.

We have now shown that several ways of tracing photons can be devised and that each gives the same distribution of power or energy at a detector some distance away from the source. However, an intelligent choice of the photon tracing algorithm can greatly reduce the needed computations. Moreover, the computation differences depend on the particular problem, e.g., on detector distance from the source as shown here (or on the IOPs, not shown here). However, we can do still more to reduce computation times, which leads us to the next topic of variance reduction techniques.

Comments for Photon Tracking:

Loading Conversation