Collimated Whole Volume Light Scattering in Homogeneous Finite Media

Crepuscular rays form when light encounters an optically thick or opaque medium which masks out portions of the visible scene. Real-time applications commonly estimate this phenomena by connecting paths between light sources and the camera after a single scattering event. We provide a set of algorithms for solving integration and sampling of single-scattered collimated light in a box-shaped medium and show how they extend to multiple scattering and convex media. First, a method for exactly integrating the unoccluded single scattering in rectilinear box-shaped medium is proposed and paired with a ratio estimator and moment-based approximation. Compared to previous methods, it requires only a single sample in unoccluded areas to compute the whole integral solution and provides greater convergence in the rest of the scene. Second, we derive an importance sampling scheme accounting for the entire geometry of the medium. This sampling strategy is then incorporated in an optimized Monte Carlo integration. The resulting integration scheme yields visible noise reduction and it is directly applicable to indoor scene rendering in room-scale interactive experiences. Furthermore, it extends to multiple light sources and achieves superior converge compared to independent sampling with existing algorithms. We validate our techniques against previous methods based on ray marching and distance sampling to prove their superior noise reduction capability.


INTRODUCTION
I NDOOR scenes in video games, AR, VR and other interactive experiences are constrained by the room geometry. Destruction, fire, dust and smoke are gameplay elements that are often incorporated in those experiences. The microscopic particles created in these events can fill up the entire room and lead to light scattering. This phenomenon is observed as light shafts stretching from windows, doorways or other elements of the geometry.
Most previous solutions targeting real-time rendering focused on scattering in outdoor scenarios and did not account for the room geometry. The overall effect of these solutions is that the outdoor environment is greatly obscured by dense fog which might be undesired look for a certain experience. Furthermore, the rendering algorithm should incur low computational cost to deliver noise-and artifact-free images within the allotted rendering time.
Integration of the light scattered towards a camera can be expensive since an unbiased solution can potentially involve solving an integral with infinite number of dimensions [1], expressed as a Neumann series. Most real-time applications limit the evaluation to a single scattering event. Integration is performed by taking a set of samples along the camera ray defined by pinhole perspective projection. Visibility and transmittance at each scattering event along the camera ray is determined by performing multiple intersection tests against the scene geometry and evaluating the distance traveled through the medium to reach a light source.
Many techniques were proposed in the real-time rendering community to efficiently estimate the single-scattering integral of radiative transfer [2]. They range from numerical integration by a Monte Carlo integrator or ray marching for estimation based on projection on Fourier basis functions or moments [3], [4] which can be further filtered. These basic frameworks are derived by separating geometry from radiance evaluation.
By considering the many facets of integration of single scattering in indoor scenes, we propose a set of techniques for sampling and integration of radiance. We explain most elements of our techniques with respect to rectilinear boxshaped medium. However, these ideas are general enough to extend to other geometry and they already seamlessly integrate within existing real-time game engines. Our central idea is to use the remarkably simple solution of unoccluded radiance in a medium constrained by a polygonal shape to derive 1) an algorithm for analytic integration of the unoccluded single scattering in box-shaped medium, 2) importance sampling algorithm matching exactly the distribution of the geometry-dependent transmittance along a camera ray, 3) an optimized Monte Carlo integration method exploiting variable caching for computing the radiance scattered inside the medium from a collimated light source towards the camera.
We are now going to clarify the nature of those three main contributions. Each of them can be considered separately since they are applicable within very different rendering frameworks.

Analytic Integration of Unoccluded Single Scattering and Ratio Estimators
The types of lighting conditions discussed in this work involve directional lights, used to approximate illumination by a very distant light source such as a star illuminating a planet, and collimated light sources masked by a gobo. We exploit the assumption of collimated light to construct a skewed coordinate system with axes corresponding to a pair of light and camera rays and origin at the camera position. We intersect the box-shaped medium against the plane of this coordinate system and form linear segments. We will call them trapezoid segments for brevity. Radiance is analytically integrable, assuming homogeneous medium within those segments. We apply this integration in estimators that split the computation into two terms: unoccluded singlescattering term and an ratio term accounting for visibility. Our integration method solves the first term, while we use existing techniques to estimate the ratio between the occluded and unoccluded radiance in the second term. We investigate integration with single-scattering ratio estimators and approximate estimators based on moments [4]. This framework forms our first integration approach. The application of a ratio estimator enables standalone use of this technique since it can be easily combined with existing sampling schemes, such as equidistant and distance sampling algorithms. Previous approaches used it for shadowing BRDFs [5] or considered it as an approximation when in fact it can asymptotically converge to the correct solution of the integral. We furthermore show that it generalizes to multiple scattering and progressive rendering.

Whole Volume Importance Sampling
Our second contribution is an importance sampling scheme applicable in Monte Carlo integrators. Sampling is carried according to the geometry of the entire medium, thus incorporating both camera and light rays in the computation. This differs from classic techniques that limit sampling only according to transmittance along the current light or camera ray, depending on the integrator. The resulting sampling strategy can deliver images with lower variance compared to existing methods. This technique is applicable on its own in path tracing frameworks where existing approaches enable combining multiple sampling strategies [6].

Optimized Monte Carlo Integrator
The third contribution regards the implementation of the integrator using our sampling scheme. Our method is split into two steps: geometry traversal and segment selection. The geometry of the medium is successfully broken down into trapezoid-shaped segments by traversing the edges of the box from corner to corner to connect the incident and outgoing exit points from the medium. The repeated traversal for each sample can be avoided by caching variables describing each segment. Based on this data a segment can be selected and samples can be drawn proportional to transmittance accounting for both the light and camera ray in paired fashion. This is the main optimized algorithm which we developed with real-time rendering application in mind. The final computation is significantly simplified and delivers performance improvements over ray marching. Mixed reality room-scale experiences are one potential target for our technique. The additive nature of light shafts allowed us to successfully deploy the technique in a non blocking optical see-through AR setup (Fig. 2). We further see our techniques as general enough to incorporate them in video games and even as a sampling strategy in path tracing frameworks.

OVERVIEW
The main problem being solved is how to efficiently integrate radiance within a well-defined medium's bounds Fig. 1. We present a set of algorithms for efficiently computing single scattering by a box-constrained medium which are applicable to rendering realtime volumetric lighting effects in indoor scenarios. Our techniques correctly separates the indoor environment from the rest of scene while rendering at real-time frame rates with improved performance and variance reduction (a). Furthermore, they generalize to multiple scattering and can improve convergence in scenes with multiple collimated light sources (b). Fig. 2. AR application. Our technique deployed to an optical seethrough AR headset experience as an additive lighting effect where virtual light shines through virtual windows onto the real world. The finite extent of the medium allows the distant buildings in the environment to be still visible and not obstructed by a thick fog. specified as a rectilinear box-shaped medium (Fig. 3). We are going to use this medium representation to clearly separate indoor and outdoor geometry allowing the distant landscape to be still visible from a room in our environment ( Fig. 2), but retain the compelling effect of light scattering in a room filled with microscopic dust and smoke particles.
The main concern when performing this procedure is how to avoid corner cases and reuse computations. Scattering within a participating medium is governed by the equation of radiative transfer [2]. We will consider paths representing single-scattering events, which connect each point along the camera ray to the light source. Furthermore, we determine the bounds of the integration domain w.r.t. the medium geometry. In the case of a directional light source and a pinhole camera, the entire integration domain lies on a single plane for each individual pixel on the imaging plane of the camera sensor.
The plane is defined by the direction of the camera ray and light source. We will call each of these planes a lightview plane. We outline the section of the box formed by intersecting the (yellow) light-view plane against the (blue) medium geometry in Fig. 3. The shape of the section can be a point or a polygon with three to six edges. Moreover, scattering happens only in the (green) section between the camera and closest occluder or medium exit point. Combining those factors makes it non-trivial to construct an algorithm that handles all corner cases and performs all operations analytically. We chose to solve this problem by iteratively breaking down the section into trapezoid-shaped domains that allow analytic integration. The iterative loop is executed by selecting the starting point p L which is the exit point of the ray starting at the camera origin and pointing towards the light source. The algorithm then traverses the outline of the section until it reaches the exit point of the camera ray p C .
The vectors parallel to each segment of the boundary can be determined in advance in the coordinate system of the box where each of the three parallel pairs of faces is orthogonal to one of the axes. This is accomplished by finding the vectors that are both orthogonal to the specific axis and lie on the light-view plane 1 . We use each q in the process of traversing the outline to determine the next direction that lies on the adjacent face. Their direction is selected to progress the computation towards the exit point of the camera ray. Each of these vectors is parallel to the orange arrows in Fig. 3 and they outline the direction of traversal.
The specified geometry allows analytic integration for unoccluded geometry by following the integral form of the radiative transfer equation (RTE) [2]. However, shadowing needs to be further applied to recreate the appearance of crepuscular rays cast from doorways, windows and other geometry masking out light. We are going to discuss three techniques to achieve it. The first approach takes the analytic integration of the single-scattered radiance in a section of the box determined by the camera ray and directional light source and pairs it with a ratio estimator or approximate estimator based on projection on moments. Both approaches use the same basic framework of analytically computing the unoccluded radiance and then numerically integrating the ratio between unoccluded and occluded radiance. In the case of momentbased estimation it uses a heuristic formula which leads to an approximate solution, while ratio estimators with a sampling-based numerical integrator are capable of  Table 1 for symbols explanation.  1. We name these face vectors as they lie on a specific face without being explicitly assigned to a certain position.
computing an asymptotically converging correct solution. Even in multiple scattering cases.
The second approach exploits the fact that the integral over all segments can be analytically inverted. Thus an efficient sampling strategy can be derived that accounts for scattering within the entire section of the box. The derived sampling function and probability density function (PDF) can be used in general path tracing frameworks.
Finally, an efficient integrator is developed that exploits the sampling method to analytically integrate radiance using caching to avoid executing at every iteration the geometry traversal algorithm outlined in Fig. 3. Its optimized equation allows cancellation of the exponential part of the integrand since the probability density function in use matches that component perfectly. We used this integrator for real-time rendering applications on the GPU.
The following section will proceed with explanation of how our methods fit within ratio estimators and Monte Carlo integrators (Section 3) used to solve the equation of radiative transfer [2]. Following those general concepts, the discussion proceeds with geometry traversal and integration of individual segments of the box section (Section 4). Sampling is derived afterwards by inverting the integral analytically (Section 5). Finally, a Monte Carlo integrator is defined using the sampling algorithm (Section 6).

THEORETIC FOUNDATION AND RELATED WORK
Light scattering in uncorrelated participating medium is traditionally modeled using the radiative transfer equation (RTE) [1], [2]. Graphics rendering systems estimate the solution of its integral form [1], [7], [8]. Volume rendering systems based on path tracing [8] solve light transport by alternating between generating path segments representing a scattering event and then connecting them to a light source, and selecting a new direction of light propagation with probability proportional to a phase function. Photon tracing systems [9] solve radiative transfer by accumulating the contribution of a set of points spread throughout the medium, and a system working at interactive rates was shown by Jimenez et al. [10].
One common simplification is to treat the medium as locally or completely homogeneous. Then path segments are generated by drawing distance samples from distribution with exponential probability density function proportional to transmittance along the current path segment, also known as distance sampling [1], [8]. However, this simplification also gives rise to more sophisticated algorithms tracing well-defined primitives, such as points, beams [11], planes and volumes [12], even surfaces [13]. Unconverged images using those custom integrators show structured artifacts in the shape of the primitive which is less desirable than noise introduced by path tracing approaches.
In real-time graphics most methods are constrained to single scattering. We first assume homogeneous medium and ignore light emitted by the medium, leading to the equation The total radiance Lðp o ; v o Þ arriving at the camera is a sum of the in-scattered radiance Lðp o À tv o ; v i Þ from all directions S at each point along the current ray and radiance arriving after scattering off a surface L s ðv i ; v o Þ -both exponentially attenuated by the medium according to Beer's law based on the extinction s t and distance. The ray is selected to end either at the surface t S or exit point from the medium t C depending on what is closest ðt SC ¼ minðt S ; t C ÞÞ. Radiance at each scattering event is weighted by the scattering coefficient s s and phase function fðv i ; v o Þ, depending on the incident v i and outgoing v o light direction. Single scattering means that only direct light illumination is computed at each scattering event  Fig. 3) v o ;ṽ o outgoing light vector in world and box coordinates (cf. Fig. 3) p o ;p o camera origin in world and box coordinates (cf. Fig. 3) d L distance to exit point from origin along light ray (cf. Eq. 4) p L ;p L light ray medium exit point (cf. Eq. (4)) p 1 ; p 2 ; p 3 positions of corners in world coordinates (cf. Fig. 3) t C distance to exit point along camera ray (cf. Eq. (5)) p C ;p C camera ray medium exit point (cf. Fig. 3) t C;1 ; t C;2 distance to corner projected on camera ray (cf.  1)) L s radiance scattered by a surface (cf. Eq. (1)) L i radiance emitted by a light source (cf. Section 3) number of samples (cf. Eqs. (2) and (3)) K cached variables (cf. Section 4.2) ðxÞ generator of uniformly distributed values (cf. Eqs. (2) and (3)) T integral over transmittance (cf. Section 4.2) q face vector (cf. Section 2) d 1 . . . distance between corner and camera ray (cf. Section 4.4) d E;S ; d E;E distance to corner at start and end of trapezoid (cf. Section 4.4) t, t 0 distance along camera ray (variable and randomcf. Eqs. (2) and (3)) d m variable distance to medium exit point (cf. Eqs. (2)) P probability (cf. Section 5.2) r random value (cf. Eqs. (2) and (3)) The distance to the medium boundary d m ðtÞ leads to spatially-varying transmittance dependent on its geometry. Another common assumption that we will follow is uniform illumination ðL i ¼ constÞ. Multiple light sources can be linearly combined. See Appendix A, which can be found on the Computer Society Digital Library at http://doi. ieeecomputersociety.org/10.1109/TVCG.2021.3135764 for a detailed derivation. The traditional method for solving this equation is using Monte Carlo integration employing importance sampling. The integrand is evaluated at a fixed number of samples and weighted by the probability density of each sample ðPDFðt 0 ÞÞ, The visibility function V ðt; v i Þ accounts for occlusion by opaque geometry. Ray marching [14] is considered as the most basic method of integrating single-scattering along a ray by performing equidistant sampling using the truncated RTE. The offline rendering community has extensively studied sample generation in participating media. The most straightforward sampling approach is distance sampling based on inversion of the transmittance term which is performed in closed form for homogeneous media [1]. It works with bidirectional path tracing [15] and Metropolis light transport [16]. Kulla and Fajardo [17] proposed equi-angular sampling which distributes samples closer to the light source in participating media. More advanced techniques were developed in the case of connecting multiple paths throughout participating media [18]. We propose another approach to integrate the radiative transfer by using importance sampling. Our sampling strategy, which we call whole volume distance sampling, accounts for the complete shape of the medium when distributing samples along the camera ray. In comparison, traditional distance sampling generates samples with probability proportional only to transmittance along the current (camera) ray. Other approaches consider different integration methods of this equation with improved performance. Analytic solutions and approximations exist for a camera and point light source immersed in unoccluded participating media [19], [20]. Very early works rely on shadow volumes for integration [21]. In the real-time community different implementations exist of this basic approach [22], such as rendering shadow planes [23] and interleaved sampling with the programmable shading pipeline [24]. Moro et al. [25] further optimized shadow plane placement and used 3D textures to store cloud densities. Wyman and Ramsey [26] proposed a hybrid between ray marching and shadow volumes employing bilateral upsampling in the shadowed regions. Munoz [27] extensively studied integration strategies to derive a set of algorithms with improved convergence properties. Sampling in a set of points and gradient-based interpolation was studied for rendering heterogeneous medium [28]. Radial basis functions (RBF) were used for low-frequency approximation of density fields and illumination by low-frequency environment light sources [29]. Other approaches use voxelized shadow volumes for scattering [30], [31]. More specialized methods were proposed to reparameterize the integral in epipolar coordinates [32], [33] or linearly rectifying the integration domain [3]. Epipolar coordinates were further used with a min-max tree to speed-up traversal [34]. We outline an algorithm in our work based on Monte Carlo integration that exploits the analytic solution of unoccluded radiance under the abovestated assumptions in a finite medium to cancel out the transmittance component and estimate only visibility.
Ratio estimators are another approach to estimate radiance based on unoccluded radiance estimate and a ratio of the visibility divided by the sum of all weights. Previously, Heitz et al. [5] applied Monte Carlo ratio term on unoccluded area light approximations. Predating this work, techniques based on filtering shadow maps computed volumetric single scattering [3], [4] using an approximate ratio term and exact unoccluded term. In our work we derive and study variations of ratio estimators that asymptotically converge to the correct solution in the form of a product between unoccluded analytic expression or Monte Carlo estimate, and the ratio of two Monte Carlo estimates (occluded and unoccluded). The more widely discussed and related to previous works is its single-scattering variant, where the distance ðt 0 ¼ Sððk=ðN sample À 1ÞÞÞÞ is generated by an importance sampling function SðxÞ which transforms samples from uniform distribution ðxÞ to the target distribution. Estimation is separated in two terms. The first term is the unoccluded integral T box ðv i ; v o ; t S Þ. It is essentially an integral over transmittance at each point along the camera ray where light source radiance L i is separated out. Our analytic computation of the single-scattering radiance in a box section corresponds to that term (Section 4). The second term is a ratio between the current occluded and unoccluded estimate. Evaluation is performed numerically at a finite number of samples N sample . In the case of completely unoccluded section, the estimator converges in a single sample, since the ratio term, which is dependent on visibil-ityV ðv i ; v o ; t 0 Þ, equals 1. Ratio estimators can be further generalized to multiple scattering by evaluating both the unoccluded and occluded light transport component. A generic formula for the analytic unoccluded term in that particular case does not exist for complex medium, therefore it can be analytically integrated after each scattering event and combined with next scattering events using stochastic methods. Other surfaces can be handled either by folding them into the ratio term or splitting the evaluation using the linearity of light transport. More details are explained in the supplemental material, available online, since this step does not significantly alter the ratio estimator framework. It has to be noted that progressive rendering requires accumulation of three surfaces containing both unoccluded terms (estimated and analytic) and the occluded term. Otherwise, the denominator in the estimator will not be properly canceled and it will lead to bias. Further information about potential bias and convergence behavior is explained in the supplemental material, available online.
The previously proposed estimators using Fourier coefficients [3] and moments [4] compute the ratio termV based on heuristics. These approaches were inspired by earlier works on filtering shadow maps [35]. In those previous works computation is targeted towards outdoor scene rendering. Estimators of radiance scattered in the atmosphere were extensively studied in the graphics literature [36], [37], [38], [39]. Our techniques, on the other hand, are focused on indoor scenes constrained by the room geometry. We use the integral separation according to a ratio estimator using both numerical and approximate moment-based estimation [4] and combine it with our unoccluded single-scattering solution in rectilinear box medium.

INTEGRATION OF UNOCCLUDED RADIANCE
We are now ready to discuss our first contribution. The first method performs radiance integration in a section of the box enclosed between the camera origin and the view and light ray. Integration is only performed by evaluating transmittance at each scattering event along the camera ray. The radiance integration method is going to be central in our ratio estimator-based solutions and some of its intermediate results are essential to explain the importance sampling scheme. We are first going to discuss the coordinate system as it is essential to understand how the algorithm is being evaluated. Then we are going to discuss the main loop of the algorithm to give a perspective why the different steps are required and then we are going to explain the underlying procedures in detail. We provide pseudo-code for all algorithms in Appendix E, available in the online supplemental material.

Coordinate System
All operations are performed in box coordinates where the x-, y-and zaxis of each point are within range ½À1; 1. If the camera origin lies outside the box, it is moved to the first intersection with the box, or the contribution is assigned to 0 when the camera ray does not intersect the box. The orientation, scale and translation of the box medium is expressed as a matrix B 3Â4 . The origin of the light ray is transformed in box coordinates ðp o ¼ B 3Â4 ½p o 1 T Þ, and the incident light vector is similarly expressed in non-normalized box coordinates ðṽ i ¼ B 3Â4 ½v i 0 T Þ. The intersection of the light ray with the box determines the starting point of the traversal algorithm, The intersection algorithm based on the slab test [40], [41] is explained in Appendix B, available in the online supplemental material. Since that information doesn't change, it can be assigned as a shader constant in the final implementation. The exit point along the camera ray is needed as a destination point of the traversal, as seen in Fig. 3. It can be computed by rasterization of the box or in a similar fashion in box coordinates ðṽ o ¼ B 3Â4 ½v o 0 T Þ by performing an intersection,

Integration Algorithm
The exit points along the light ray and camera ray give us the starting and end point and determine the boundaries of the box section that will be traversed by the algorithm. The main objective is to break down the box section defined by the intersection of the light-view plane against the box into multiple segments (Fig. 3). Since the light source is assumed to emit uniform light, the integral can be split into a product of integral over transmittance and emitted radiance from the light source. Three is the number of iterations performed by the algorithm which corresponds to the maximum number of edges involved in a box section. The main explanation for this fact is that a box is composed of three parallel planes and any set of rays constrained to a line segment inside the box can intersect at maximum three faces. We are going to outline the main steps in each iteration and discuss them in detail in their own subsections. Each iteration starts by finding the next corner. Then the distance to the corner and along camera ray is computed. They are used to integrate radiance analytically for a trapezoid segment, The trapezoid is geometrically defined by the distance of the starting point (t C;S ) and end point (t C;E ) along the camera ray combined with distance from the starting point (d E;S ) and end point (d E;E ) to the edge of the box. Each iteration advances the distances along the camera ray and recomputes the distance to the edge. Transmittance of all segments is summed to derive the final integral (T total ¼ P 3 k¼1 T k ). The sampling methods and optimized Monte Carlo algorithm use this procedure. They further pre-cache a set of variables (K) and the individual integral values for each trapezoid segment (T ).

Selecting the Direction of Traversal to Next Corner
The direction of traversal of the boundary outline is defined as the intersection of the light-view plane with the box which leads to progress towards the exit point along the camera ray. We will call face vectors the direction of each line formed by this computation. They are associated with a pair of parallel faces. We compute them in advance before traversing the outline to determine the position of each light segment. They are computed to be perpendicular to each axis (x-, y-and zaxis) and the light-view plane. The direction of these vectors is selected to follow the direction of progress towards the camera ray exit point or the shortest path to the same point. The direction of progress is determined in skewed coordinates determined by the incident and outgoing light direction vector. The exit point is projected on the outgoing light direction and its sign is taken as a direction of progress. Detailed derivation is provided in the supplemental material and the appendix, available in the online supplemental material. The default case is to multiply by the progress along the camera ray ðv ¼ g progress ða; Àṽ o ÞÞ. When the light plane is parallel to the camera ray direction, traversal is not performed, instead the integral over transmittance is directly computed. Therefore, this case is not considered by the algorithm. Another possibility is that the vector connecting the two points is parallel to the light plane ðjvj < Þ. Then the direction along the camera ray is selected instead ðv ¼ Àða Áṽ o ÞÞ. The final face vector is determined by multiplying the value representing the direction of progress (q ¼ v a).
The resulting face vectors are used when selecting the next direction of traversal. At each corner there is a choice between one of three directions corresponding to each face vector associated with each axis. A face vector is selected only if the current position lies on an edge that connects a face corresponding to it and progress can be made over the remaining coordinates by traversing that face vector. The rest of the algorithm finds the smallest distance that must be traversed based on all coordinates given the selected face vector.

Computing the Next Corner
The next algorithm forms trapezoid segments by simultaneously computing in the skewed light-view coordinate system the length of the two line segments connecting its projection on the camera ray with the corner and camera origin. First, overstepping past the exit point along the camera ray (p C ) is prevented by finding the progress along the edge and clamping it to fall within the line segment t cut ¼ min 1; g progress Àp C Àp;w Á À Á À Á . Then the end point along the edge is computed (p edge ¼p þ t cutw ). Next, the projected distance along the camera ray is computed À t C;E ¼ g progress Àp edge Àp o ; Àv o ÁÁ in a similar manner to form the distance vector between the line and the edge. The result of the computation is the distance along the camera ray (t C;E ). Combining it with the origin (p o ) and outgoing light direction (v o ) yields the final point (p cam p o À t C;Eṽo ). The distance from the camera ray and the edge is computed as the difference between the point on the camera ray and edge which is projected on the scaled incident light ray vector The incident light vectorṽ i;s 2 is corrected following the procedure in Appendix C, available in the online supplemental material.

Radiance Computation in a Trapezoid Segment
We now outline how radiance is integrated for each segment. Each segment is a trapezoid section defined by two parallel bases along the incident light vector (ṽ i ) and two legs containing the line segments: along the camera ray and medium boundary. Essentially, a linear segment in a skewed coordinate system. The domain has a slope, The final expression is a product of transmittance to the start of the segment and integral over transmittance in the linear segment, which has remarkably simple solution and continues to be compatible with distance sampling, e Às t c t dt ¼ e Àðt C;S þd edge;S Þ s t À 1 À e Às t ðt C;B Àt C;S Þ c Á ¼ E l ð1 À E u Þ: (6Þ The distance along the camera ray is clamped by the distance to the surface to account for occlusion along the camera ray ðt C;B ¼ clampðt S ; t C;S ; t C;E ÞÞ. In the limit when the slope approaches zero (c ! 0), the computation turns into a linear equation ðT ¼ ðt C;B À t C;S Þ E l Þ. The same expression is used to avoid numerical issues when the size of the segment approaches zero (t C;E À t C;S < ). Note, that it is possible to reduce the number of exponent evaluations by extracting the product of the lower and upper bound part (E l Á E u ) outside of the evaluation and assigning it as lower bound part at each consecutive evaluation (E l ). There are two cases when the incidentṽ i and outgoing lightṽ o vectors are parallel and require to be solved separately. The same procedure applies without traversing the outline. When they are facing opposite directions (ṽ i Áṽ o < 0), the end distance goes to zero, and when they are facing the same direction it is the sum of the distance to both exit points along the incident d L and outgoing light vector t C , The integration algorithm as explained requires an estimator capable to account for shadowing. We chose to implement both a ratio estimator and an approximate momentbased estimator in our work. Another approach is to use Monte Carlo integration. We are going to explore it in the following and compare it against moment-based estimation and ratio estimator in the Results section.
Integration can be trivially extended to a collimated light source with a gobo (convex binary map or stencil) by offsetting the starting point of the integration to lie on the gobo edge and limiting it to the far end of the gobo. The result is exponentially attenuated by the distance from the camera origin to the start of the gobo. Multiple lights can be also combined thanks to the linearity of light transport. An example of both effects is illustrated in Fig. 1. We use our unoccluded radiance computation to compute weights for the multiple lights sampling following Monte Carlo integration similar to Shirley and Wang [42]. We demonstrate this novel generalization to rendering light scattered in finite medium in Fig. 1b. More details regarding gobo masking and multiple light integration are available in the supplemental material, available online.

WHOLE VOLUME DISTANCE SAMPLING
The two important elements required to incorporate a sampling strategy in path tracing or single-scattering integrators are a sampling function and probability density function (PDF). The sampling function distributes samples according to the probability defined by the PDF. In the following we discuss both functions separately.

Importance Sampling Function
Sampling according to transmittance dependent on medium geometry is achieved in two steps. The first step is to build a histogram. The integration of transmittance of each trapezoid segment is according to the procedure explained above. Trapezoid is selected using comparison against the cumulative density function (CDF) of all segments. Considering the geometry along the entire camera line segment is what separates our technique from existing distance sampling techniques. CDF is computed by summing the integral of transmittance up to that segment ( P k T k ) and dividing it by the unoccluded integral of transmittance ðT total . . .Þ ¼ T box ðṽ i ;ṽ o ; t S Þ. Then the appropriate sampling function for each trapezoid segment is computed, where the random value r 2 ½0; 1 is transformed after selecting the segment. Its value is multiplied by the unoccluded integral T total , then the integral over previous segments ( P i k¼1 T kÀ1 ) is subtracted and the whole result is normalized by the integral value of the current segment T i , The outcome is an intermediate random value in range ½0; 1. The next step is to determine the distance along the camera ray within the segment. Position of corners (p 1 ,p 2 , p 3 ) and distance vectors between corners (w 1 ,w 2 ,w 3 ) are required to determine the distance between camera ray and box boundary (Eq. (9)). The computation of these components (Section 4.4) is performed similarly to the analytic radiance integration algorithm (Section 4.2). The final sampling method accounts for distance to the edges (d E;S , d E;E ) and distance along the camera ray (t C;S , t C;E ) in a paired fashion. It differs from distance sampling that accounts only for transmittance along the current path segment [8]. The proposed importance sampling scheme is equivalent to distance sampling only when evaluating a single segment where the distance to the boundary is constant (d E;S ¼ d E;E ). More details are available in the supplemental material, available online. Our importance sampling considers the connection from each path segment to box boundary which is bound by trapezoids. The bases are the path segments of each pair of vertices along the camera ray to the exit point along the light ray. The legs are the box boundary and camera ray. Computation requires similarly the slope (c), bounding of the distance along the camera ray (t C;B ) and upper bound component of the trapezoid segment radiance integral (E u ) computed like previously explained. The inversion method is used to derive the sampling component. Random values r are selected by transforming from uniform distribution to distance by a scaled logarithm function ðÀlog 1 À ð1À ð E u Þ rÞ=ðs t cÞÞ. Note, that the values differ for each trapezoid segment and they are dependent on geometry. In the case when the slope approaches zero (c ! 0), the sampling simplifies to multiplication by the distance between starting and end point. The result of the algorithm is distance offset by the starting distance t C;S , representing distance from the origin along the camera ray. Visibility can be determined by sampling towards the light source from the point at that distance. Successive accumulation is used to solve the Monte Carlo integral (Eq. (2)).

Probability Density Function
The PDF matches perfectly the exponential component in the integral equation and cancels it (Eq. (2)). The probability of selecting each trapezoid segment is equal to the ratio of radiance scattered within the segment and the maximum (unoccluded) radiance ðP ðx 2 ½t C;iÀ1 ; t C;i Þ ¼ T i =T total Þ. Taken on its own this function is defined over continuous domain (x 2 ½t C;iÀ1 ; t C;i ) and fits the definition of probability density function (PDF). The conditional probability of selecting a given distance, provided that a segment is selected, is equal to the transmittance at the particular point divided by the integral over transmittance within the segment ðP ðxjx 2 ½t C;iÀ1 ; t C;i Þ ¼ T 0 i =T i ¼ e Às t ðxþdmðxÞÞ =T i Þ. We can assemble the final probability density function by considering where the point x is with respect to the three projected corners of the box on the camera ray (t C;1 ; t C;2 ; t C;3 ), We provide a straightforward validation of the equation in Appendix D, available in the online supplemental material. Further validation is performed by comparing the sample distribution against the probability density function. Examples are shown in Figs. 4 and 5. The probability density is transmittance over the unoccluded integral in both cases. The probability density of the sampling function is computed numerically by using stratified random samples from uniform distribution. These samples are collected in 256 bins and normalized against the sample count to produce the final probability density function. It is visible that the numerical and analytically computed PDF match in both cases. The function is also affected by geometry and it is not necessarily monotonically decreasing as seen in the first example (Fig. 4). The presence of a corner leads to a function that is not smooth with a clear difference in its rate of change around the distance of the projected corner (Fig. 5). A distance sampling function will fail in the cases when the exit point along the camera ray is much closer than the distance from the origin towards the light source (Fig. 4). It will also deviate in all cases where the distance to the box boundary is not constant (d E;S 6 ¼ d E;E ).

OPTIMIZED MONTE CARLO INTEGRATION
Finally, we provide an optimized integrator incorporating our sampling strategy with variable caching to reduce the number of redundant computations. Monte Carlo integration according to our implementation is performed by connecting a set of light rays along the camera ray. Each vertex is selected by performing importance sampling derived by choosing a trapezoid segment and inverting the analytic integration algorithm. Repeated traversal of the outline of the box is avoided in the final algorithm by caching values. Unoccluded integration follows the algorithm from Section 4 and caches values in a matrix K. The cumulative distribution histogram is built by first computing the probability of selecting a segment ðG k ¼ T k =T total Þ. Then probabilities are summed up to the given segment (H k P k i¼1 G i ). When performing sampling, the random values must be rescaled to fit the range of the selected trapezoid segment. The multiplicative weight is computed as the inverse probability (D k ¼ 1=G k ). While the summand is the previous CDF histogram value scaled by the multiplicative term (M k ÀH kÀ1 D k ). Combined with the cached values (K), it is enough to compute sample location based on random values selected from uniform distribution ðr ðk=N sample ÞÞ. Monte Carlo integration is performed over a number of samples (N sample ). Each sample from random distribution (r) is compared against the step functions of the histogram ðs stepðH 1 ; rÞ þ stepðH 2 ; rÞÞ. The sample is then transformed in segment coordinates (r s r D s þ M s ). Distance is sampled according to the importance sampling (Section 5.1). However, the upper bound integral part (K s;1 ¼ 1 À E u ), slope (K s;2 ¼ s t c), size of segment along camera ray (K s;3 ¼ t C;B À t C;S ), starting position along the camera ray (K s;4 ¼ t C;S ) are all cached in advance. The exponential function is cancelled by the probability density function. Therefore, the integration loop accumulates only visibility The final integral over visibility V total is re-weighted by the number of samples N sample , scattering coefficient s s and phase function fðv i ; v o Þ to compute the solution of the radiative transfer integral (L ¼ 1 N sample s s fðv i ; v o Þ V total T total L i ). Only in the single-scattering case a ratio estimator with the same importance sampling scheme yields the same integration algorithm.

IMPLEMENTATION AND RESULTS
Now that we have a complete picture of how our techniques are built from theoretic perspective, we can discuss the practical implementation details and how they impact execution performance and quality. Furthermore, we compare our estimators and samplers against existing techniques and provide details about how they integrate and perform in more practical approximate frameworks targeted at highperformance graphics.
Multiple integration algorithms were outlined within this work. Our aim was to use them in a real-time application built on Unity. Thus we have implemented everything as Unity scripts using the shading language in this engine. Our implementation does not rely on special features from this engine, so it can be ported to any rendering framework. Additionally, we provide a supplemental HTML page using JavaScript and WebGL for rendering that will assist reproduction of our work. We did not re-implement our technique in any offline rendering framework, but as we stated, we see it as a viable sampling strategy for that kind of application.
The stochastic integration algorithms using our techniques are: Monte Carlo with whole volume distance sampling (MCWVDS), ratio estimator with distance sampling (REDS) and ratio estimator with equidistant sampling (REES). MCWVDS is implemented following the optimized algorithm explained in Section 6. Ratio estimators use the basic framework defined in Section 3 paired with the analytic unoccluded radiance computation outlined in Section 4. REDS incorporates traditional distance sampling to generate samples proportional to transmittance along the current path segment, while REES distributes samples at equal distance which is jittered by noise. Those stochastic integrators converge to the same solution of singlescattering. Additionally, we implemented an approximate estimator based on moments following the method described by Peters et al. [4]. We defer to the code bundled with this previous publication for more details. This method relies on heuristics, so regardless of the sample count, it will deviate from the correct solution when dealing with complex occluders. Finally, we implemented ray marching (RM) which distributes samples at equal distance and jitters them with random offsets and distance sampling (MCDS) distributing samples according to transmittance along the current path segment and ignoring the overall shape of the medium. Both methods directly intersect the box at each step of the evaluation following the algorithm discussed in Appendix B, available in the online supplemental material to determine the distance from the scattering event to the box boundary along the light ray.
Performance Benchmark. All integration techniques are benchmarked in Table 2 at 1080p resolution. We captured times on three different GPUs. Our MCWVDS integration approach depends on a logarithmic function which is not directly supported by the GPU and in the corner cases it fallbacks to linear segment sampling through branching. Regardless, these computations are less demanding than performing an intersection test to compute the distance between the scattering event and box boundary and then evaluating an exponential function used by ray marching. We noticed clear performance improvements on NVIDIA GPUs -possibly related to how the cached values are indexed. Precise comparison between the code emitted for different GPU ISAs is still difficult with the currently available tools. We used the GPU profiler in Unity to measure only passes of our technique. Our AMD setup (*) was on a MacBook Pro laptop where we captured the complete frame times which lead to high timing variance. GPU profiling in Unity was not available on that OS. It is also possible to implement a variant of ray marching which uses our technique to find the corners and replaces the intersection test  [4] provides further performance improvements because it approximates the integral with lower number of dimensions by evaluating a heuristic formula dependent on depth and converting the result into moments which are linearly integrated.
Integrator Comparison. We evaluated how our technique integrates with different estimators on a simple scene with a single box medium and a sphere occluder (Fig. 6). We compared it against jittered ray marching (RM) which uses a slab test to find the medium boundaries at each step (Fig. 6b). This example is considered as previous work since it relies on existing approaches. Overall, it fails to distribute samples proportional to transmittance. Our whole volume distance sampling technique MCWVDS (Fig. 6a) follows closely the transmittance integrand. In the unoccluded sections, it converges in a single sample, while in the rest of the scene, it leads to noticeable variance reduction. Ratio estimator using our unoccluded analytic integral computation (Algorithm 2) and ratio term using ray marching (REES) yields better convergence in the unoccluded parts of the scene, but leads to overall poor distribution of samples in the occluded parts. Distance sampling is another previous work (Fig. 6d). It distributes samples proportional to the distance to the observer, however that does not perfectly match the transmittance integrand component, leading to higher variance. Combining this technique with a ratio estimator and our analytic unoccluded radiance integration (Fig. 6e REDS), leads to improved variance reduction because in certain areas the distribution of transmittance follows more closely exponential decay. Ratio estimators are general enough that it is possible to combine them with more sophisticated integration techniques [43]. We demonstrate multiple scattering with our importance sampling and ratio estimators (Figs. 6f, 6g ,6h, 6i, and 6j) and show their potential as integrators in path tracing frameworks. These example figures were generated with the HTML page bundled in the supplemental material, available online.
Integrator Convergence Rate. Fig. 6 shows results at a fixed number of samples per pixel. We further studied how the algorithms behave when increasing the sample count in Fig. 7a. Overall, our optimized integration algorithm (MCWVDS) achieves the fastest convergence rate. Ratio estimator based techniques similarly achieve lower variance in unoccluded areas which leads to higher convergence rate. Finally, ray marching and distance sampling that are considered as previous techniques perform the worst. Distance sampling performs good at very low sample count, but the fact that it relies on very simple assumptions regarding the geometry leads to lower convergence rate at higher sample count. In the case of multiple scattering distance sampling (Fig. 7b), both ratio estimators and our importance sampling algorithm lead to lower variance. The sample count reflects the number of camera rays and after each scattering event only a single shadow ray connect with the light source. The resulting integral corresponds to a single sample when considering only single-scattering integration spread over multiple path segments. Thus distance sampling retains its advantage over ray marching, even when multiple camera ray samples are taken into account. However, augmenting it with a ratio estimator or using our sampling technique leads to noticeable variance reduction.
Destruction Scene. We tested our MCWVDS integration   Fig. 6. Our main technique using whole volume importance sampling achieves the fastest convergence rate. Ratio estimators combined with our unoccluded radiance computation follow since they converge quickly in unoccluded areas. Distance sampling in that case is closest to the distribution of transmittance, but overall misses important parts of the geometry which places it in second place and results in lower convergence rate when it is used with Monte Carlo integrator. Ray marching performs poorly at lower sample count, but it can outperform distance sampling at higher sample count, when considering single scattering.
algorithm on a dynamic scene involving destruction (Fig. 8) and compared it against ray marching (RM). The relatively simple geometry of the occluders and high extinction contribute to the overall higher variance reduction introduced by our algorithm compared to ray marching. We further evaluated our technique MCWVDS against standard research test scenes. Results for the breakfast room scene are shown in Fig. 9. Samples are preferentially distributed close to the camera and thus resulting in noise reduction. The more heavily occluded scene results in less variance reduction than the destruction scene ( Fig. 8) because it leads to deviation of the integral from the probability density of the importance sampling function. As a rule of thumb our sampling approach has advantage in rooms with higher extinction, while ray marching allows taking more samples in unoccluded sections at the back of the room. Benchmark of the rendering performance is shown in Table 2.
We implemented also a variant of the moment-based single-scattering estimation [4] with our analytic unoccluded radiance integration technique (Algorithm 2). We call it moment-based estimation (ME). Since it relies on our technique, we don't consider it completely as previous work. It was tested in a very simple outdoor scenario where a set of basic shapes lie outside of the medium and thus cast detached shadows (Fig. 10). In that very simple scenario, the moment-based estimation matched really closely the results of the integration.
Fireplace Scene. We further tested the moment-based estimation technique against the fireplace room (Fig. 12) scene where the quality was still acceptable, but overestimation artifacts became more apparent. We noticed significant performance improvements when using this technique (cf. Table 2).
Other Medium Shapes. We demonstrate in Fig. 11 that the principles outlined for a box-shaped medium also extend to other convex shapes. Concave shapes are also possible, but they will require more complex integration procedure. We expand the discussion in our supplemental material, available online.
Failure Case of Moments-Based Estimation. In the case of the breakfast room scene (Fig. 13), the moment-based estimation significantly deviates from the approximated integral and leads to severe artifacts. We do not use noise metrics, such as PSNR and RMSE, for those comparisons since the artifacts have very visible structure caused by the heuristic formulas used by this approach [4]. The original paper prescribes a set of heuristic terms to correct for overestimation issues. However, they fail to remove the issue and might further result in false shadows behind thinner objects which is equally distracting as an artifact. Possible future work is to enhance the heuristic to better match the integral.     Supplemental Video. We provide a supplemental video with further comparisons of our techniques when the camera is in motion. Convergence rate of the optimized Monte Carlo algorithm (Algorithm 8) is demonstrated by using a set of monotonically increasing number of samples. Furthermore, we provide a fully dynamic scene to outline the behavior in interactive scenes of our Monte Carlo integration algorithm and our analytic unoccluded radiance integration algorithm paired with the approach described by Peters et al. [4]. Moment-based estimation might exhibit artifacts when the heuristic deviates from the integral, but provides better performance. While our approach based on Monte Carlo integration leads to noise proportional to the number of samples. Denoising algorithms alleviate this problem. We defer to existing surveys for more details [44].

DISCUSSION
We presented a set of algorithms for integration and sampling that are directly applicable to Monte Carlo integrators and ratio estimators in indoor scenes. We demonstrated them predominately with a rectilinear box-shaped medium and showed how to extend them to generic shapes. Together, they can represent many scenes in video games, AR and VR experiences. Ratio estimators were used to shadow the unoccluded radiance computation by uniform collimated light which is part of our first main contribution. They provide greater flexibility regarding combining different sampling strategies, but lower convergence rate than our final technique. We are unaware of any prior work based on this kind of stochastic estimator for finite participating media which we believe is a valuable contribution on its own. The unoccluded radiance computation was further used in a moments-based approximate estimator that has better performance characteristics than all proposed techniques, but leads to visual artifacts caused by its heuristic nature. We leave as future work any perceptual studies that would evaluate whether those artifacts are objectionable enough to users of systems using our techniques.
We further derived a Monte Carlo sampling strategy which takes into account the entire shape of the medium and we showed the equivalence of estimators using it to a ratio estimator only when considering single-scattering integration. That is an important result showing where both theoretic frameworks lead to equivalent expressions. In the multiple scattering case they significantly deviate since the individual terms of the ratio estimator include complex integrals. The superior convergence rate of our sampling strategy was quantified in a series of experiments (c.f. Figs. 6 and 7). Our application targeted real-time performance and we optimized the integrator for this particular case, however the sampling strategy is compatible with path tracing frameworks, which makes it a contribution on its own. The final optimized Monte Carlo integrator has the potential to outperform previous approaches even at the same sample count (c.f. Table 2) and it is our final contribution.
One major application of our approach is room-scale AR where we surrounded the entire room by the volume to enable compositing of light shafts. The main advantage being that the windows to the outside virtual world are not affected by scattering and allow gameplay elements to be visible even at a great distance outside the room. This configuration was a useful building block to deliver the AR experience shown in Fig. 2.
Limitations. We derived our techniques under the assumption of homogeneous media, collimated light, uniform illumination and sampling decisions based on a single scattering event. We will provide in the next few paragraphs ideas about how to resolve those limitations in the future.
We have shown how to handle convex media, however, concave shapes will require more meticulous splitting to allow the use of our sampling technique (Algorithm 7). Ratio estimators can be directly applied using the signed area to determine the sign for each linear segment when analytically computing the integral over transmittance. The approach can be extended to sampling heterogeneous medium by breaking it down into multiple shapes.
Single-scattering estimators based on Fourier coefficients or moments require the unoccluded term to be computed for all geometry contributing to the radiance scattered towards the camera. Taking the opposite route of starting with semi-infinite and subtracting the contribution of the environment outside the box leads to overall more complex equations which we demonstrate in the supplemental material, available online. Those approximations can consequently lead to objectionable artifacts. Thus, we do not perceive this idea as viable future direction. Furthermore, moment-based estimators in general are more likely to fail with complex geometry since they rely on heuristics to compute the single-scattering integral as seen in Fig. 13. In the future, we are interested in developing heuristics and integration schemes to prevent these artifacts.
We demonstrated our sampling strategies in single-and multiple-scattering scenarios using the main principles of path tracing frameworks. Thus, we are able to leverage existing techniques to incorporate other types of luminaires by exploiting the linearity of light transport.
Finally, we demonstrated throughout our work how to apply finite media estimators to our specific problem. However, we think that those general principles can be further developed and applied to novel estimators and approaches to solving radiative transfer in finite media and they will be a valuable starting point for future research.