Page updated: February 15, 2021
Author: Curtis Mobley
View PDF

Ray Tracing

The previous page showed how to determine ray 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 ray paths through an absorbing and scattering medium.

Ray Tracing

There is more than one way to simulate ray 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 rays propagate through a medium. Loosely speaking, a ray 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 ray survival in any particular interaction. When a ray encounters a particle, we can randomly decide if the ray 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 ray is scattered.
  • If 𝔯 > ωo, then the ray is absorbed.

If the ray is absorbed, tracing stops and a new ray is emitted from the source and tracing begins anew. If the ray 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 rays, which also introduces the geometry to be used in numerical simulations below. A source emits a collimated beam of rays, which are then recorded by an annular, ring, or ”bullseye” detector some distance away. The red ray is emitted by the source, undergoes one scattering, and is then absorbed by a particle. The green ray is emitted, undergoes two scatterings, and is recorded by a detector.


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

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

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 ray, as just described. Rather than tracing one ray at a time as nature does, consider a source emitting “bundles” of many rays (often called “photon packets” in the literature). Then view each interaction as having a fraction 1 ωo of the rays in the bundle being absorbed, and the remaining fraction ωo being scattered, all in the same direction. Let the bundle be emitted with an initial weight of w = 1, which can represent one unit of energy, power, or some number of rays. At each interaction, the current weight w is multiplied by ωo to account for the loss of energy or number of rays by absorption (that is, a fraction ωo continues onward). The scattered bundle then carries a reduced weight. If the ray bundle reaches the target, the current weight w is tallied. Another bundle is then emitted from the source and traced. This tracing process, which we’ll call ”Type 2,” is illustrated by the green ray track in Fig. 2. After two scatterings, as illustrated, the detected ray bundle has weight w = ωo2.


Figure 2: 2. Illustration of Type 2 ray 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 rays in the bundle are viewed as being scattered. Then, if the ray bundle reaches the target, absorption is treated as a continuous process occurring along the entire ray 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 ray 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, ray bundle by ray bundle. However, numerical simulation of many ray bundles shows all three of these ray 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 rays are tracked, and rays can be absorbed.
Type 2:
Ray bundles are tracked, with a bundle weight being multipled by ωo at each interaction.
Type 3:
Ray bundles 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 rays, 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 rays being sent from the source and using Type 1 ray tracing. rays were traced until they were absorbed. Figure 3 shows some of the resulting statistics.

The red histogram shows the percent of rays that traveled an optical distance τ1 τ τ2 between interactions, for a bin size of τ2 τ1 = 1. The theoretically expected fraction of rays 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 ray 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 rays were absorbed. Thus the value for the first bar shows that about 18% of the rays 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 ray 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 rays 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 rays 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 rays. Because oceanic phase functions scatter much more light at small scattering angles than at large angles, most rays that are scattered just a few times will hit the detector near its center. To even out the numbers of rays (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 rays crossing the detector plane at some distance r from the detector center are tallied in bins as follows:

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

Simulations were made with Nemit = 106 rays (for Type 1) or ray bundles (Types 2 and 3) emitted from the collimated source. Figure 4 shows the distribution of rays or bundles 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 rays (depending on the tracing type) reach the target plane without being scattered. Most rays are scattered 3 or 4 times, and very few rays are scattered more than 10 times. The right panel shows that for zT = 15, almost no rays reach the target plane unscattered, and most undergo 5 to 25 scatterings, with a peak around 10 or 15, depending on the way the rays are traced. For Type 1 ray tracing, almost no rays 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 ray bundles absorbed, there are broad tails in the number of scatterings, although for Type 2 a bundle being scattered 40 times has an almost negligible weight of w = ωo40 = 0.00013.


Figure 4: Distribution of rays 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 rays reaching a τ = 10 × 10 area of the target plane for the first 104 emitted rays, tracing Type 1, and the target plane at zT = 5. Note that of the 104 emitted rays, only 2973 reached the target plane (of which 2966 are in the area plotted). This is less than 30% of the rays making a contribution to the answer of how much power is detected; i.e., 70% of the calculations were wasted. The ray-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 rays were scattered more than about 30 off of the optical axis. The right-hand panel shows the distributions of the rays reaching the target plane by number of scatterings and total distance traveled. Every ray must of course travel a distance of at least τ = 5 before reaching the target plane, and most rays were scattered several times.


Figure 5: Distribution of rays reaching the target plane at zT = 5 for 104 emitted rays and Type 1 tracing. The left panel shows the spatial distribution of points where the rays 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 ray bundles eventually intersect the target plane. Only 7% of the emitted rays were wasted. Those rays 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 rays 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 rays ever reached the target plane. Over 98% of the ray-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 rays used to computed the statistics. For Types 2 and 3 about 71% and 82%, respectively, of the emitted rays eventually reach the target plane. The statistical noise is now much smaller (but still noticeable) because of the larger number of ray bundles.

The total optical distance distributions for Types 2 and 3 show broad tails. In both cases, fewer than one fourth of the rays made it to the target plane after traveling a total distance of τ = 1516. About 30% of the rays 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 rays 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 rays 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 rays, the second is scattered rays inside the first detector ring, and the last plotted point is for rays 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 rays (open circles) are different for the three tracing types. For Type 1, the distribution of the number of rays is the same as the power distribution because the rays all retain their initial weight of w = 1. Thus power detected is simply the number of rays detected. For Types 2 and 3, more rays 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 rays (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 rays emitted. There are obvious differences in the distributions for the three tracing types. However, if Nemit = 106 rays are traced, as in Fig. 13 these differences almost disappear. This indicates that the three ways to trace rays all give identical predictions of the detected power, to within some amount of statistical noise, which can be reduced by tracing more rays.

However, the computation time required by the three tracing types can vary greatly. Recall from Fig. 5 that about 30% of the rays 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 rays to achieve some desired level of statistical noise, we would have to emit and trace over three times as many rays 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 rays reached the target plane (Fig. 8), whereas about 80% of the rays reached the target plane for Types 2 and 3 (Figs. 9 and 10). Getting the same number of rays on target would thus require emitting over 40 times as many rays (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 rays 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 ray 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 Ray Tracing:

Loading Conversation