Volume Rendering by Scratchapixel

This page should be merged with the “Volume Rendering Draft” one. This is content that is missing from the “Volume Rendering Draft” page which should be the final one.

When traveling through a volume, light is partially absorbed and scattered by it

Ray-Marching Algorithm

  1. Find intersections [math]t_0[/math] and [math]t_1[/math] (camera ray enters/leaves volume)
  2. Divide [math]t_0[/math]-[math]t_1[/math] segment into [math]X[/math] subsegments of identical size
  3. March along the camera ray, for each subsegment:
    • Shoot a light ray from the middle of the subsegment to the light source, compute the distance it travels in the volume and use Beer’s law to get amount of light at that point:
      • [math]L_i(x) = e^{- \text{distance} \cdot \sigma_a} \cdot \text{light colour} \cdot dx[/math]
    • Light is attenuated as it passes through the sample itself, compute the transmission value of this small volume element too
      • [math]T_x = e^{-dx \cdot \sigma_a}[/math]
    • Total contribution of this sample is
      • [math]T_x \cdot L_i(x)[/math]
  4. Combine each sample for its respective contribution to the overall opacity and colour of the volumetric object but
    • Farthest volume element is occluded by the one before along the ray, and so on -> backward (t1 to t0) or forward (t0 to t1) ray-marching

Backward Ray-Marching

[math]\text{final colour} = \text{background colour} \cdot \text{transmission} + \text{result}[/math]

Forward Ray-Marching

[math]\text{final colour} = \text{background colour} \cdot \text{transmission} + \text{result}[/math]

Forward ray-marching allows for early stopping when transparency of volume gets close to 0 (large volume or scattering coefficient)

Better Ray Marching

Interactions between light beams and particles in a medium:

  • Light beams traveling to the eye through the medium lose energy due to:
    • Absorption: light energy absorbed by particles of the medium
    • Out-scattering: light that’s traveling towards the eye is scattered out on its way to the eye
  • Light beams traveling to the eye through the medium gain energy due to:
    • Emission: flame can emit incandescent light
    • In-scattering: some of the light that’s not initially traveling toward the eye is being redirected toward the eye due to scattering

Scattering coeff

Need to account for both absorption and out-scattering in computation for how much light we loose -> need to account for scattering [math]\sigma_s[/math] in Beer’s law alongside absorption coeff [math]\sigma_a[/math]

[math]T = e^{- \text{distance} \cdot (\sigma_a + \sigma_s)}[/math]

Sometimes note extinction coeff as [math]\sigma_t = \sigma_a + \sigma_s[/math]\

We need to use the above when computing transmission value and also

How much light is being scattered towards the eye due to in-scattering is also proportional to the scattering term

So we need to account for it in a bunch of stuff, figure out the formula

Density term

We have assumed that scattering and absorption coefficients are uniform across the media -> homogeneous participating medium. Need to know how to handle a heterogeneous participating medium too, more frequent IRL.

We want some vairable to scale our scattering and absorption coefficients: density

Phase function

In-scattering contribution for a sample computed using (this is the running “result” in the ray marching algorithm):

[math]L_i(x, \omega) = \sigma_s \int_{S²} p(x, \omega, \omega’) L(x, \omega’) d\omega'[/math]

  • [math]L_i[/math]: in-scattering contribution
  • [math]x[/math]: sample position
  • [math]\omega[/math]: viewing direction (camera ray direction) pointing from the obj to the eye
  • [math]\omega'[/math]: light direction, from the obj to the light
  • [math]L(x, \omega’)[/math]: light contribution, amount of light that’s coming from a particular light direction (here [math]\omega'[/math] at sample point [math]x[/math] after it has traveled a certain distance through the volume
  • [math]p(x, \omega, \omega’)[/math] is the phase function
  • What do we integrate on? All directions over the entire sphere [math]S^2[/math] of directions
  • Isotropic scattering volume: when a photon interacts with a particle, it can be scattered out in any direction within the sphere of possible directions around the particle where every direction is equally likely to be chosen
  • Anisotropic scattering: scatters light in restricted range of directions
  • Phase function tells you how much light is being scattered for a particular combination of directions: view direction [math]\omega[/math] and incoming light direction [math]\omega'[/math].
    • Returns value in [0, 1], models angular distribution of light scattered
    • Integrates to 1 over its domain: the sphere of directions
      • [math]\int_{S^2} f_p(x, \omega, \omega’) d \omega’ = 1[/math]
      • considering all directions from which light cna come around the particle, how much light is being scattered out around that same particle can’t be greateer than the sum of all incoming light -> phase function needs to be normalised over the sphere of directions
    • Reciprocity: result is the same if you swap [math]\omega[/math] and [math]\omega'[/math]
      • [math]f_p(x, \omega, \omega’) = f_p(x, \omega’, \omega)[/math]
    • Only depends on angle between view and incoming light direction, often defined in terms of that angle [math]\theta[/math]
    • So basicallz tells you how much light is likely to be scattered towards the viewer [math]\omega[/math] for an incoming light direction [math]\omega'[/math]
  • Isotropic volumes phase function
    • [math]f_p(x, \theta) = \frac{1}{4 \pi}[/math]
    • Not dependent on incoming light direction and view, and all directions are equally likely to get picked
    • Area of a sphere is [math]4\pi[/math] steradians, so the phase function has to be this value in order to integrate to 1
    • Unit of phase functions is 1/steradian
  • Henyey-Greenstein phase function
    • g: asymmetry factor (ranges from -1 to 1): controls whether light is scattered in the forward or backward direction (when 0 the function is the isotropic one)

Jittering the sample positions

Taking sample always in middle of the segment can lead to banding -> pick random position on each segment instead

Early termination when opaque

breaking out from the ray-marching loop as soon as you detect that the transparency variable is lower than a minimum threshold, can only be done in forward ray marching

But this is statistically wrong: we cut some of the contribution of the volume that’s below a certain threshold

More accurate estimation: russian roulette:

  • when transparency value is lower than some threshold we pick a random number in [0, 1] and test whther it’s greater than 1/d (d is some positive real number greater than 1).
    • If it is, we break from the loop, otherwise we continue but we multiply the current transparency value by d
    • d represents the likelihood that we pass the test (if d=5, the changes of the loop being terminated would be 4 out of 5)

Rendering a 3D Density Field (Heterogeneous Medium)

Heterogeneous volumes (smoke, cloud) have parts more opaque than others: absorption or scattering coeff vary through space

More practical approach: constant value for scattering and absorption coeffs and use density parameter to modulate the volume appearance through space: have a density function that returns the density at a point in space

How to generate the density field?

  • Procedural: use 3D texture (Perlin noise) to procedurally create a space-varying density field
  • Simulation: use fluid simulation program that simulates motion of fluids (smoke) to generate it

Generating a density field using a noise function

  • Procedural noise function: function that procedurally generates a noise pattern through 3D space
    • Take point as argument and returns value of the 3D noise texture at that point (real number [-1, 1])
    • Need to remap it to [0, 1] to modulate density (density is in [0, 1])
  • Just move the definition of the density value inside the ray marching loop and set it to the value of procedural noise function at the sample position -> sample the density as we march along the ray
  • This is how the noise function influences the transmission over distance (plot perlin noise function over distance and transmission over same distance)

We need to make some changes to the way we compute in-scattering though

In-scattering for a heterogeneous participating medium

  • Density varies through the medium as the incoming light ray traverses it to get to the camera ray -> cant use the Li we used above
  • Density also varies along the camera ray

We used raymarching above to estimate the in-scattering along the camera ray

Now we will also use it to estimate the camera and light rays transmission as they travel through a heterogeneous participating medium: break the ray into a series of segments, and estimate the transmission of each one of the segments, assuming that over the segment’s length, the density of the volume element is uniform, and then multiply the total transmission value by the sample’s transmission as we move along the ray

now that we use a procedural texture we may encounter some filtering issues. If the frequency at which you sample the noise function is too low, you may miss some details from the procedural texture and you will eventually get aliasing issues

Volume Rendering based on 3D Voxel Grids

TODO

From Radiative Transfer Equation to Volume Rendering

Derivation of Beer-Lambert Law

Derivative of the radiance for a light beam starting at [math]x[/math] traveling a set [math]s[/math] along the direction [math]\omega[/math]: Rate at which radiance is lost due to absorption:

[math]dL = – \sigma_a L(x, \omega) ds[/math]

Assume loss of radiance due to absorption only to start with, we’ll generalise later.

[math]L(x, \omega)[/math] is the radiance of the light beam at point [math]x[/math] along the direction vector [math]\omega[/math]. It decreases as the distance traveled by the light through the medium increases

Rate at which radiance changes is proportional to the radiance itself at [math]x[/math] where proportionality factor is absorption coeff

[math]L(x, \omega)[/math] is a function, this funciton is the beer-lambert law itself: it decreases as the distance traveled by the light through th medium increases (as [math]x[/math] gets further away from the point where the light beam enters the medium). So we are trying to figure out what [math]L(x, \omega)[/math] is using the derivative equation above.

Let’s first rewrite it as a function of [math]s[/math] (distance that light beam has traveled through the medium from entry point [math]x[/math]) rather than a function of [math]x[/math].

We can substitute [math]x_s = x + s \omega[/math] to get:

[math]\frac{dL(s)}{ds} = – \sigma_a L(s)[/math]

We get an ordinary differential equation (derivative and derivate’s function appear in the same equation so we can apply the following rule:

[math]\frac{dy}{dx} = cy /rightarrow y = e^{cx}[/math]

So:

We get as a result the equation for the Beer-Lambert law, which works if the medium is homogeneous. Can extend that to:

Beer’s Law, Transmittance and Optical Depth

Transmittance: measure of the volume’s opacity, measure of how much light can pass through a volume, fraction of light that is transmitted through the volume:

[math]T = \frac{L_o}{L_i}[/math]

Can be computed using Beer’s law:

[math]T = e^{- \sigma_t s}[/math]

[math]T = e^{- \tau}[/math]

  • [math]s[/math]: distance between two points in the volume
  • [math]\tau[/math]: optical depth/thickness
  • Internal transmittance: accounts for absorption only
  • Total transmittance: accounts for absorption, out-scattering, …

For heterogeneous volume, extinction coeff varies through space -> integrate it along the ray:

[math]\tau = \int_{s=0}^{d} \sigma_t(x_s) ds[/math]

[math]d[/math]: distance traveled by ray through the volume

Final and most general form of Beer’s law:

[math]T(d) = \text{exp} \( – \int_{s=0}^{d} \sigma_t(x_s) ds \)[/math]

In-scattering and the Phase function

Phase function helps in estimating in-scattering: how much of the energy from the light beam passing through the cylinder at some oblique angle is scattered in the [math]- \omega[/math] direction (towards the eye)

Gives fraction of incoming light traveling along [math]\omega'[/math] that’s scattered towards [math]- \omega[/math]

Participating medium can have two types of scattering behaviour:

  • Isotropic: all directions in the sphere of directions are equally likely to be chosen
  • Anisotropic: some directions in the sphere of directions are more likely to be chosen than others (forward, backward, …)

Radiative transfer equation and volume rendering equation

RTE says that these elements contribute to a change of radiance in energy flow direction [math]\omega[/math]: absorption, in-scattering, out-scattering, emission (omitted here). Defines a change of radiance (a derivative) along direction [math]\omega[/math]:

[math] \frac{L(x + s\omega, \omega)}{ds} = – \sigma_t L(x, \omega) + \sigma_s \int_{S^2} f_p(x, \omega, \omega’) L(x, \omega’) d \omega’ [/math]

[math]- \sigma_t L(x, \omega) [/math] accounts for absoprption and out scattering, the other term for in-scattering

[math]s[/math]: positional change in direction [math]\omega[/math]

Note the [math]\sigma_s[/math] in front of the integral: amount of in-scattering proportional to proba that light is being scattered as we defined a dL above

Let’s define

[math]L_s(x, \omega) = \int_{S^2} f_p(x, \omega, \omega’) L(x, \omega’) d \omega'[\math]

RTE is a first-order non homogeneous linear differential equation of the form:

[math]y’ + p(x)y = q(x)[/math]

Need to solve equation in which both funciton y and its derivative y’ are present

Redefine the radiance function as a function of the distance s

[math]y \text{ corresponds to } L(s)[/math]

[math]y’ \text{ corresponds to } \frac{dL(s)}{ds}[/math]

[math]p(x) \text{ corresponds to } \sigma_t[/math]

[math]q(x) \text{ corresponds to } \sigma_s L_s(s)[/math]

[math]y’ + p(x)y = q(x) \text{ corresponds to } L'(s) + \sigma_t L(s) = \sigma_s L_s(s)[/math]

This standard form has the following known solution

[math]y(x) = \int_{t=0}^{x} q(t)e^{- \int_t^x pdt’}dt + C_1e^{-\int_{t=0}^x pdt}[/math]

Replacing, we get:

  • 3D position x is reparameterised by s’ and s” measuring distance along the ray starting from x: x_s=x+s\omega
  • Outer integral (variable s’) runs from entry point x in the volume to the exit point x_s (the point where we want to get a value for the light that has been traveling to that point), so over a distance s, represents the accumulation of light scattered into the ray’s direction along its path
  • Inner integral (variable s”) runs from current scattering point x_{s’} to end point x_s of the ray (so over a distance s’), representing attenuation of light from scattering point to end point
  • the last term accounts for initial radiance L(0) attenuated over entire path length s

This is the Volume Rendering Equation

If we were to consider emission we would add emission term L_e next to in-scattering term:

VRE turns the RTE into an integral which can be solved using Riemann sum

With transmittance of the medium over a distance s:

We can rewrite the VRE as

Stochastic Method for Monte Carlo