Ray Tracing

Yet another rendering method

  • Rasterization is an object-based approach to scene rendering.
    • Loop over triangles and see which pixels are filled
  • Ray tracing, by contrast, colors the pixels first, then identifies them with objects later.
    • Loop over pixels and see which triangle fill the pixel

Ray tracing is:

  • usually more costly,
  • but far more flexible than rasterization
  • Basic operation is to find the closest intersection between a ray and objects

Pseudocode of Ray Tracing

1
2
3
4
5
6
7
8
for all pixels {
  ray = generate_camera_ray( pixel ) 
  for all objects {
  	hit = intersect( ray, object )
  	if (“hit” is closer than “first_hit”) {first_hit = hit} 
  }
  pixel = shade( first_hit ) 
}

Ray-Sphere Intersection

…

Ray-Box intersection

Slab Test

Efficient ray-axis aligned box test

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
bool intersect(HitInfo& minHit, const Ray& ray) const {
  // set minHit.t as the distance to the intersection point
  // return true/false if the ray hits or not
  float tx1 = (minp.x - ray.o.x) / ray.d.x;
  float ty1 = (minp.y - ray.o.y) / ray.d.y;
  float tz1 = (minp.z - ray.o.z) / ray.d.z;

  float tx2 = (maxp.x - ray.o.x) / ray.d.x;
  float ty2 = (maxp.y - ray.o.y) / ray.d.y;
  float tz2 = (maxp.z - ray.o.z) / ray.d.z;

  if (tx1 > tx2) {
    const float temp = tx1;
    tx1 = tx2;
    tx2 = temp;
  }

  if (ty1 > ty2) {
    const float temp = ty1;
    ty1 = ty2;
    ty2 = temp;
  }

  if (tz1 > tz2) {
    const float temp = tz1;
    tz1 = tz2;
    tz2 = temp;
  }

  float t1 = tx1; if (t1 < ty1) t1 = ty1; if (t1 < tz1) t1 = tz1;
  float t2 = tx2; if (t2 > ty2) t2 = ty2; if (t2 > tz2) t2 = tz2;

  if (t1 > t2) return false;
  if ((t1 < 0.0) && (t2 < 0.0)) return false;

  minHit.t = t1;
  return true;
}

Ray Tracing in Triangle

Ray tracing and transforms

What if your object is transformed by affine transformation?

  • Two options:
    • Inverse transform a ray, compute the intersection, then transform the intersection
    • Transform the object and compute the intersection

Option 1:

  • Bring a ray into the object space (transform the origin and the direction)
  • Compute intersection in the object space
  • Transform the intersection (position and normal)
  • Intersection can be simplified, but it can be slower
  • Option 1 can avoid storing transformed objects

Option 2:

  • Transform the object into the world space
  • Compute intersection in the world space
  • No transformation is needed for intersection
  • Sometimes difficult to do intersection (e.g., sphere becomes an ellipsoid), but can be faster

Ray Marching

Ray marching is a class of rendering methods for 3D computer graphics where rays are traversed iteratively, effectively dividing each ray into smaller ray segments, sampling some function at each step.

Ray-Triangle Intersection (What we do in practice)

Given a ray and a triangle, the objective is to compute (t,α,β,γ)(t, \alpha, \beta, \gamma)

Way 1: Cramer’s Rule (Better)

Using Cramer’s Rule

https://www.geeksforgeeks.org/system-linear-equations-three-variables-using-cramers-rule/

p⃗=αa⃗+βb⃗+γc⃗\vec{p} =\alpha\vec{a}+\beta \vec{b}+\gamma \vec{c} is the Parametric description of a point in a triangle

p⃗=o⃗+td⃗o⃗+td⃗=αa⃗+βb⃗+γc⃗α=(1−β−γ)o⃗+td⃗=(1−β−γ)a⃗+βb⃗+γc⃗ox+tdx=(1−β−γ)ax+βbx+γcxoy+tdy=(1−β−γ)ay+βby+γcyoz+tdz=(1−β−γ)az+βbz+γcz\begin{aligned} \vec{p} & = \vec{o}+t \vec{d} \\ \vec{o}+t \vec{d} & =\alpha\vec{a}+\beta \vec{b}+\gamma \vec{c} \\ \alpha &=(1-\beta-\gamma)\\ \vec{o}+t \vec{d} & =(1-\beta-\gamma) \vec{a}+\beta \vec{b}+\gamma \vec{c} \\ \\ o_x+t d_x & =(1-\beta-\gamma) a_x+\beta b_x+\gamma c_x \\ o_y+t d_y & =(1-\beta-\gamma) a_y+\beta b_y+\gamma c_y \\ o_z+t d_z & =(1-\beta-\gamma) a_z+\beta b_z+\gamma c_z \end{aligned}

[ax−bxax−cxdxay−byay−cydyaz−bzaz−czdz][βγt]=[ax−oxay−oyaZ−oz]\left[\begin{array}{lll} a_x-b_x & a_x-c_x & d_x \\ a_y-b_y & a_y-c_y & d_y \\ a_z-b_z & a_z-c_z & d_z \end{array}\right]\left[\begin{array}{l} \beta \\ \gamma \\ t \end{array}\right]=\left[\begin{array}{l} a_x-o_x \\ a_y-o_y \\ a_Z-o_z \end{array}\right]

Solve the equation with Cramer’s Rule to find β,γ,t\beta, \gamma, t.

Cramer’s Rule https://rosettacode.org/wiki/Cramer's_rule:

Given

{a1x+b1y+c1z=d1a2x+b2y+c2z=d2a3x+b3y+c3z=d3\left\{\begin{array}{l} a_1 x+b_1 y+c_1 z=d_1 \\ a_2 x+b_2 y+c_2 z=d_2 \\ a_3 x+b_3 y+c_3 z=d_3 \end{array}\right.

which in matrix format is

[a1b1c1a2b2c2a3b3c3][xyz]=[d1d2d3]\left[\begin{array}{lll} a_1 & b_1 & c_1 \\ a_2 & b_2 & c_2 \\ a_3 & b_3 & c_3 \end{array}\right]\left[\begin{array}{l} x \\ y \\ z \end{array}\right]=\left[\begin{array}{c} d_1 \\ d_2 \\ d_3 \end{array}\right]

Then the values of x,yx, y and zz can be found as follows:

x=∣d1b1c1d2b2c2d3b3c3∣∣a1b1c1a2b2c2a3b3c3∣,y=∣a1d1c1a2d2c2a3d3c3∣∣a1b1c1a2b2c2a3b3c3∣, and z=∣a1b1d1a2b2d2a3b3d3∣∣a1b1c1a2b2c2a3b3c3∣.x=\frac{\left|\begin{array}{lll} d_1 & b_1 & c_1 \\ d_2 & b_2 & c_2 \\ d_3 & b_3 & c_3 \end{array}\right|}{\left|\begin{array}{lll} a_1 & b_1 & c_1 \\ a_2 & b_2 & c_2 \\ a_3 & b_3 & c_3 \end{array}\right|}, \quad y=\frac{\left|\begin{array}{lll} a_1 & d_1 & c_1 \\ a_2 & d_2 & c_2 \\ a_3 & d_3 & c_3 \end{array}\right|}{\left|\begin{array}{lll} a_1 & b_1 & c_1 \\ a_2 & b_2 & c_2 \\ a_3 & b_3 & c_3 \end{array}\right|}, \text { and } z=\frac{\left|\begin{array}{lll} a_1 & b_1 & d_1 \\ a_2 & b_2 & d_2 \\ a_3 & b_3 & d_3 \end{array}\right|}{\left|\begin{array}{ccc} a_1 & b_1 & c_1 \\ a_2 & b_2 & c_2 \\ a_3 & b_3 & c_3 \end{array}\right|} .

For efficient computation:

https://math.stackexchange.com/questions/3734743/connection-between-cross-product-and-determinant

det⁡(a⃗,b⃗,c⃗)=(a⃗×b⃗)⋅c⃗\det(\vec{a}, \vec{b}, \vec{c}) = (\vec{a} \times \vec{b}) \cdot \vec{c}

where ⋅\cdot is dot product, ×\times is cross product.

We accept the solution only if

tclosest>t>01>β>01>γ>01>1−β−γ>0\begin{aligned} t_{\text{closest}}> & t > 0 \\ \\ 1 > & \beta > 0 \\ 1 > & \gamma > 0 \\ 1 > 1 - & \beta - \gamma > 0 \end{aligned}

where tclosestt_{\text{closest}} is the smallest positive tt values so far

Way 2: Signed volume

Using the Triple Product, we can get the Signed volume of parallelepiped

V=A⃗⋅(B⃗×C⃗)V = \vec{A} \cdot (\vec{B} \times \vec{C})

where ⋅\cdot is dot product, ×\times is cross product.

1
2
3
4
float getTripleProduct(float3 A, float3 B, float3 C) const {
    float V = dot(A, cross(B, C));
    return V;
}

If V<0V < 0, it is inside the triangle, else if V>0V > 0 it is outside the triangle.

q⃗=o⃗+td⃗=o⃗+d⃗\vec{q} = \vec{o} + t\vec{d} = \vec{o} + \vec{d}

  • t=1t = 1 when define q⃗\vec{q}

If Va,Vb,VcV_a, V_b, V_c share same sign, it has Intersection.

At the same time, we can get the barycentric coordinates for free because:

α=VaVa+Vb+Vcβ=VbVa+Vb+Vcγ=VcVa+Vb+Vc\begin{aligned} \alpha & =\frac{V_a}{V_a+V_b+V_c} \\ \beta & =\frac{V_b}{V_a+V_b+V_c} \\ \gamma & =\frac{V_c}{V_a+V_b+V_c} \end{aligned}

Recall to interpolate any attribute value:

1
2
3
float getInterpolatedAttribute(float3 barycentricCoords, float3 attribute) const{
    return barycentricCoords.x*attribute.x + barycentricCoords.y*attribute.y + barycentricCoords.z*attribute.z;
}

Multiple ways to compute t:

  • Ray-plane intersection
  • Interpolated position + dot product
  • Another signed volume

Shading

Supp materials:

Reflection types:

  • Diffuse: matte paint, paper
  • Specular: mirror
  • Glossy: plastic, rough metallic surface

BRDF

  • BRDF (Bidirectional Reflectance Distribution Function)

https://boksajak.github.io/blog/BRDF

Reflected Radiance is the Sum of contributions from all directions.

Properties of BRDF:

  • Range: [0, infinity]
    • Max value is not 1 because BRDF is not reflectance, it is a concentration of reflected energy
  • Dimensions: 4 (assuming uniform material)
    • Two for incoming directions
    • Two for outgoing directions
  • Reciprocal
    • BRDF stays the same for swapped viewer/light
    • BRDF(ωo,ωi)=BRDF(ωi,ωo)BRDF(\omega_o, \omega_i) = BRDF(\omega_i, \omega_o)
  • Energy Conservation
    • Reflections do not produce extra energy

Extra: BSDF

  • BSDF (Bidirectional Scattering Distribution Function)
  • Generalization of BRDF to transparency

Lambertian

  • Uniformly reflects light over all directions
  • Constant BRDF = Lambertain BRDF
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
static float3 shadeLambertian(const HitInfo& hit, const float3& viewDir, const int level) {
        float3 L = float3(0.0f);
        float3 brdf, irradiance;

        // loop over all of the point light sources
        for (int i = 0; i < globalScene.pointLightSources.size(); i++) {
            float3 l = globalScene.pointLightSources[i]->position - hit.P;

            // the inverse-squared falloff
            const float falloff = length2(l);

            // normalize the light direction
            l /= sqrtf(falloff);

            // get the irradiance
            irradiance = float(std::max(0.0f, dot(hit.N, l)) / (4.0 * PI * falloff)) * globalScene.pointLightSources[i]->wattage;
            brdf = hit.material->Kd / PI;

            if (hit.material->isTextured) {
                brdf *= hit.material->fetchTexture(hit.T);
            }
            L += irradiance * brdf;
        }
        return L;
}

Shadow ray tracing

  • Return irradiance only if the light is visible
  • ignore intersections that are too close (set a threshold distance)

https://web.cse.ohio-state.edu/~shen.94/681/Site/Slides_files/shadow.pdf

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
static float3 shadeLambertian(const HitInfo& hit, const float3& viewDir, const int level) {
        float3 L = float3(0.0f);
        float3 brdf, irradiance;

        // loop over all of the point light sources
        for (int i = 0; i < globalScene.pointLightSources.size(); i++) {
            float3 l = globalScene.pointLightSources[i]->position - hit.P;

            // the inverse-squared falloff
            const float falloff = length2(l);

            // normalize the light direction
            l /= sqrtf(falloff);

            // get the irradiance
            irradiance = float(std::max(0.0f, dot(hit.N, l)) / (4.0 * PI * falloff)) * globalScene.pointLightSources[i]->wattage;
            brdf = hit.material->BRDF(l, viewDir, hit.N);

            if (hit.material->isTextured) {
                brdf *= hit.material->fetchTexture(hit.T);
            }
            L += irradiance * brdf;
        }
        return L;
}

static float3 shadeLambertianShadow(const HitInfo& hit, const float3& viewDir, const int level) {
    float3 L = float3(0.0f);
    float3 brdf, irradiance;

    // loop over all of the point light sources
    for (int i = 0; i < globalScene.pointLightSources.size(); i++) {
        float3 l = globalScene.pointLightSources[i]->position - hit.P;

        // the inverse-squared falloff
        const float falloff = length2(l);

        // normalize the light direction
        l /= sqrtf(falloff);

        // get the irradiance
        irradiance = float(std::max(0.0f, dot(hit.N, l)) / (4.0 * PI * falloff)) * globalScene.pointLightSources[i]->wattage;
        brdf = hit.material->BRDF(l, viewDir, hit.N);

        if (hit.material->isTextured) {
            brdf *= hit.material->fetchTexture(hit.T);
        }

        //shadow ray
        Ray shadowRay = Ray(hit.P,l*sqrtf(falloff));
        HitInfo shadowHitInfo;
        //Set tMin to a tiny value to avoid self-intersection
        //Light length was sqrtf(falloff) so it is within the length of the light vector
        if(globalScene.intersect(shadowHitInfo, shadowRay, 0.000001f, sqrtf(falloff)))
        {
            continue;
        }
        L += irradiance * brdf;
    }
    return L;
}

Specular Reflection

https://raytracing.github.io/books/RayTracingInOneWeekend.html#metal/mirroredlightreflection

For Mirror

  • Recursive tracing is used
  • Trace a new ray toward the reflected direction
    • ω⃗r=−2(ω⃗i⋅n⃗)n⃗+ω⃗i=ω⃗i−2(ω⃗i⋅n⃗)n⃗\vec{\omega}_r=-2\left(\vec{\omega}_i \cdot \vec{n}\right) \vec{n}+\vec{\omega}_i = \vec{\omega}_i-2\left(\vec{\omega}_i \cdot \vec{n}\right) \vec{n}
1
2
3
float3 reflect(const float3& incident, const float3& normal) {
    return incident - 2.0f * dot(incident, normal) * normal;
}
  • You need to set the recursion
    • level of recursion = number of bounces
    • use level to trade off of how accurate you want it to be
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
static float3 shade(const HitInfo& hit, const float3& viewDir, const int level) {
  ...
  } else if (hit.material->type == MAT_METAL) {
        if (level >= threshold) {return float3(0.0f);}
        float3 reflectedDir = reflect(viewDir, hit.N);
        // Create a reflection ray
        // Invert the direction such that -reflectedDir
        Ray reflectionRay = Ray(hit.P, -reflectedDir);
        HitInfo reflectionHitInfo;
        if(globalScene.intersect(reflectionHitInfo, reflectionRay, 0.000001f))
        {
            // Invert the direction again such that -reflectionRay.d
            //return shade(reflectionHitInfo, -reflectionRay.d, level+1);
            return hit.material->Ks * shade(reflectionHitInfo, -reflectionRay.d, level+1);
        }
        else {return float3(0.0f);}
	}
...

Specular Refraction

For Glass

  • Recursive tracing is used
  • Just now the direction become a refracted direction
    • ω⃗t=η1η2(ω⃗−(ω⃗⋅n⃗)n⃗)−(1−(η1η2)2(1−(ω⃗⋅n⃗)2))n⃗\vec{\omega}_t=\frac{\eta_1}{\eta_2}(\vec{\omega}-(\vec{\omega} \cdot \vec{n}) \vec{n})-\left(\sqrt{1-\left(\frac{\eta_1}{\eta_2}\right)^2\left(1-(\vec{\omega} \cdot \vec{n})^2\right)}\right) \vec{n}
    • Note that the thing (called as kk) in sqrt must be positive
      • k determines whether total internal reflection occurs. If k is negative, indicating total internal reflection, the function returns false.
  • Be aware of what object ray is entering and exiting
    • Entering (from η1\eta_1 to η2\eta_2): ω⃗⋅n⃗<0\vec{\omega} \cdot \vec{n} < 0
    • Exiting (from η2\eta_2 to η1\eta_1): ω⃗⋅n⃗>0\vec{\omega} \cdot \vec{n} > 0

Some terminology:

  • ior : Index of Refraction
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
static bool refract(const float3& incident, const float3& normal, float ior, float3& refractedDir) {
    float etaI, etaO, eta;
    float3 vecN;
    float cosThetaI, k;

    bool exiting; // is it exiting from etaO?

    if (dot(incident, normal) > 0.0f) {exiting = true;} else { exiting = false;}

  	// If bug happens, check whether if etaI/etaO is reversed (which one is the material (ior) ?)
    if (exiting)
    {
        vecN = -normal;
        eta = ior;
    }
    else
    {
        vecN = normal;
        eta = 1.0f / ior;
    }
    cosThetaI = dot(incident, vecN);
    k = 1 - (eta * eta) * (1 - cosThetaI * cosThetaI);
    
    if(k < 0.0f) 
    {
        return false; // Do total internal reflection
    }

    refractedDir = eta * (incident - cosThetaI * vecN) - (sqrt(k) * vecN);

    return true;
}

You need to set the recursion

  • level of recursion = number of bounces
  • use level to trade off of how accurate you want it to be

Fresnel Reflection

Environment mapping

Note that Environment mapping is not surface texture.

  • If ray hits nothing, return a value from image.

Light Probe Images

  • HDR : high dynamic range images

https://www.pauldebevec.com/Probes/

"If for a direction vector in the world (Dx, Dy, Dz), the corresponding (u,v) coordinate in the light probe image is (Dx*r,Dy*r) where r=(1/pi)*acos(Dz)/sqrt(Dx^2 + Dy^2)"

If you apply that formula you can convert a ray’s direction to uv coordinates where uv is in the range [-1.0, 1.0]

So basically:

1
2
3
4
5
float2 hdrRay2texcoords(const float3& raydir) {
  float PI = 3.14159265358979f;
  float r = (1/PI) * acos(raydir.z) / sqrt(raydir.x * raydir.x + raydir.y * raydir.y);
  return {raydir.x*r, raydir.y*r};
}

Image-based lighting (IBL)

  • environment mapping usually refers the use of images for rendering specular reflections/refractions
  • image-based lighting refers to the use of images as lighting in general (including specular and others such as Lambertian)