Image Synthesis
Intersection Shader

• This chapter introduces several intersection shaders for the Vulkan GLSL ray tracing pipeline
• The ray tracing pipeline provides a built-in intersection shader for triangles. For geometric primitives that are not triangles (such as cubes, cylinders, spheres, parametric surfaces, etc.) you have to provide your custom intersection shader
• If the intersection shader determines that a ray-primitive intersection has occurred, it can notify the ray traversal with the function reportIntersectionEXT(...)
• Furthermore, the intersection shader can fill a hitAttributeEXT variable, which can be of user-defined type and can be evaluated in the closest-hit, miss, or any-hit shader

• Parameterization of the ray with:
  $\mathbf{x}(t) = \mathbf{s} + t \,\, \mathbf{r} \quad \forall\, \, t \in [t_{\mathrm{\small min}}, t_{\mathrm{\small max}}]$
  where $\mathbf{s}$ denotes the ray origin and $\mathbf{r}$ the direction
• Parameterization of the sphere using the implicit sphere equation:
  $\begin{aligned}(x - c_x)^2 + (y - c_y)^2 + (z - c_z)^2 - r^2 &= 0\\ (\mathbf{x} - \mathbf{c})^\top \, (\mathbf{x} - \mathbf{c}) - r^2 &= 0\end{aligned}$
  where $\mathbf{c} = (c_x, c_y, c_z)^\top$ denotes the center and $r$ the radius of the sphere
• By inserting the ray equation, the intersection points can be found:
  $\begin{aligned} &(\mathbf{s} + t \,\, \mathbf{r} - \mathbf{c})^\top \, (\mathbf{s} + t \,\, \mathbf{r} - \mathbf{c}) - r^2 = 0\\ \Leftrightarrow &(t \,\, \mathbf{r} + \mathbf{b})^\top \, (t \,\, \mathbf{r} + \mathbf{b}) - r^2 = 0 \\ \Leftrightarrow &t^2 + p \,t + q = 0\end{aligned}$
  with $\mathbf{b} = (\mathbf{s} - \mathbf{c})$, $p = 2 \, (\mathbf{r}^\top \, \mathbf{b}) / (\mathbf{r}^\top\,\mathbf{r})$ und $q = ((\mathbf{b}^\top \, \mathbf{b}) - r^2) / (\mathbf{r}^\top\,\mathbf{r})$
  $\Rightarrow t_{1,2} = - \frac{p}{2} \pm \sqrt{\frac{p^2}{4} - q}$

$\mathbf{s}$

$\mathbf{r}$

$t_{\mathrm{\small in}}$

$t_{\mathrm{\small out}}$

• There can be no, one, or two intersections
  $\begin{align} t_{\mathrm{\small in}} &= - \frac{p}{2} - \sqrt{\frac{p^2}{4} - q}\\ t_{\mathrm{\small out}} &= - \frac{p}{2} + \sqrt{\frac{p^2}{4} - q}\end{align}$
• The "no intersection" case occurs if
  $\frac{p^2}{4} - q < 0$
  This condition should be checked as early as possible.
• Otherwise
  • if $t_{\mathrm{\small in}} < t_{\mathrm{\small out}}$ and $t_{\mathrm{\small in}} \in [t_{\mathrm{\small min}}, t_{\mathrm{\small max}}]$, the closet intersection point is:
    $\mathbf{x}(t_{\mathrm{\small in}}) = \mathbf{s} + t_{\mathrm{\small in}} \,\, \mathbf{r}$
  • otherwise if $t_{\mathrm{\small out}} \in [t_{\mathrm{\small min}}, t_{\mathrm{\small max}}]$, the closet intersection point is:
    $\mathbf{x}(t_{\mathrm{\small out}}) = \mathbf{s} + t_{\mathrm{\small out}} \,\, \mathbf{r}$

$y$

$x$

$z$

$\tilde{\mathbf{n}}$

$\phi$

$\theta$

• The surface normal $\mathbf{n}$ at the intersection point $\mathbf{x}(t_{\mathrm{\small hit}})$ is calculated as follows:
  $\mathbf{n} = \frac{\mathbf{x}(t_{\mathrm{\small hit}}) - \mathbf{c}}{r}$
• To calculate the texture coordinates, the normal can be transformed into the local coordinate system. In this space, the transformed normal $\tilde{\mathbf{n}}=(\tilde{n}_x,\tilde{n}_y,\tilde{n}_z)^\top$ is converted into spherical coordinates with the angles $\phi$ and $\theta$.
  The resulting texture coordinates are:
  $\begin{align} s &= \frac{1}{2\pi} \phi = \frac{1}{2\pi} \operatorname{arctan2}(\tilde{n}_y, \tilde{n}_x)\\ t &= 1.0 - \frac{1}{\pi} \theta \\ &= 1.0 - \frac{1}{\pi} \arccos(\tilde{n}_z) = \frac{1}{\pi} \arccos(-\tilde{n}_z) \end{align}$

• Let $\mathtt{T}$ be a 4 x 4 transformation matrix from object to world space and $\underline{\tilde{\mathbf{p}}}$ a point in object space in homogeneous coordinates
  $\underline{\mathbf{p}} = \begin{pmatrix} x \\ y \\ z \\ 1 \end{pmatrix} = \begin{bmatrix} m_{11} &m_{12}&m_{13} &t_x \\ m_{21} &m_{22}&m_{23} &t_y \\ m_{31} &m_{32}&m_{33} &t_z \\ 0 &0 &0 &1 \\ \end{bmatrix} \begin{pmatrix} \tilde{x} \\ \tilde{y} \\ \tilde{z} \\ 1 \end{pmatrix} = \mathtt{T} \, \underline{\tilde{\mathbf{p}}}$
• The provided matrix gl_ObjectToWorldEXT is a 3 x 4 matrix, i.e. it only contains the upper three rows of $\mathtt{T}$
• The idea for calculating the intersection point in object space is not to transform the 3D object into world space (as before), but instead to transform the ray into object space
• To this end, the inverse transformation $\mathtt{T}^{-1}$ is applied to the ray
• The provided matrix gl_WorldToObjectEXT is a 3 x 4 matrix, i.e. it only contains the upper three rows of $\mathtt{T}^{-1}$

• Parameterization of the ray in world space:
  $\mathbf{p} = \mathbf{s} + t \,\, \mathbf{r} \quad \forall\, \, t \in [t_{\mathrm{\small min}}, t_{\mathrm{\small max}}]$
  where $\mathbf{s}$ denotes the ray origin and $\mathbf{r}$ the normalized direction
• When transforming the ray into the object coordinate system, it should be noted that $\mathbf{r}$ is a direction that must not be translated. Therefore, the last component of the homogeneous coordinates is set to 0.
• This results in the transformed ray in object space:
  $\begin{align} \underline{\tilde{\mathbf{p}}} = \mathtt{T}^{-1} \underline{\mathbf{p}} = \mathtt{T}^{-1} \left(\begin{pmatrix} \mathbf{s} \\ 1 \end{pmatrix} + t \,\, \begin{pmatrix} \mathbf{r} \\ 0 \end{pmatrix}\right) &= \mathtt{T}^{-1} \begin{pmatrix} \mathbf{s} \\ 1 \end{pmatrix} + t \,\, \mathtt{T}^{-1} \begin{pmatrix} \mathbf{r} \\ 0 \end{pmatrix} \\ &= \begin{pmatrix} \tilde{\mathbf{s}} \\ 1 \end{pmatrix} + t \,\, \begin{pmatrix} \tilde{\mathbf{r}} \\ 0 \end{pmatrix}\end{align}$
• The vector $\tilde{\mathbf{s}}$ is provided by the variable gl_ObjectRayOriginEXT and $\tilde{\mathbf{r}}$ by gl_ObjectRayDirectionEXT

• Parameterization of the ray in world space:
  $\mathbf{p} = \mathbf{s} + t \,\, \mathbf{r} \quad \forall\, \, t \in [t_{\mathrm{\small min}}, t_{\mathrm{\small max}}]$
  with gl_WorldRayOriginEXT $\,\mathbf{s}$ and gl_WorldRayDirectionEXT $\,\mathbf{r}$
• Parameterization of the ray in object space:
  $\tilde{\mathbf{p}} = \tilde{\mathbf{s}} + t \,\, \tilde{\mathbf{r}} \quad \forall\, \, t \in [t_{\mathrm{\small min}}, t_{\mathrm{\small max}}]$
  with gl_ObjectRayOriginEXT $\,\tilde{\mathbf{s}}$ and gl_ObjectRayDirectionEXT $\,\tilde{\mathbf{r}}$
• The parameter $t$ is unaffected, i.e. the $t$ that is determined in this way in the object space corresponds exactly to the $t$ in world space
• Important: For transformations that contain scaling, this approach only works because $\,\tilde{\mathbf{r}}$ is not normalized to length 1 by the ray tracing pipeline. Therefore, do not normalize gl_ObjectRayDirectionEXT in the intersection shader!

• Intersection Ray-Sphere (Object Space)

struct HitAttributeType {
  vec3 normal; 
  vec2 texCoord;
}; // type of the "hit" variable

...

void main() { /**** INTERSECTION SHADER ****/

  // ray-sphere intersection in object space
  
  // get bounding box in object space
  vec3 aabbMin;
  vec3 aabbMax;
  gsnGetIntersectionBox(gl_InstanceID, gl_PrimitiveID, aabbMin, aabbMax);
  
  // center and radius of sphere in object space
  vec3 center = (aabbMin + aabbMax) / 2.0;
  float radius = length(aabbMax.x - aabbMin.x) / 2.0;

  vec3 b = gl_ObjectRayOriginEXT - center;
  vec3 r = gl_ObjectRayDirectionEXT; // not normalized
  float rr = dot(r, r);
  float p = 2.0 * dot(r, b) / rr;
  float q = (dot(b, b) - radius * radius) / rr;
  
  // discriminant of the quadratic equation
  float d =(p * p / 4.0) - q;
  
  // if no solution, cancel computation
  if(d < 0.0) return;
  
  // two solutions for quadratic equation
  float tmin = -p / 2.0 - sqrt(d);
  float tmax = -p / 2.0 + sqrt(d);
  
  // intersection at tmin or tmax?
  float t = tmax;
  uint hitKind = 2u;
  if(tmin < tmax && tmin >= gl_RayTminEXT && tmin <= gl_RayTmaxEXT) {
    t = tmin;
    hitKind = 1u;
  }
  if(t >= gl_RayTminEXT && t <= gl_RayTmaxEXT) {
    vec3 hitPoint = gl_ObjectRayOriginEXT + t * gl_ObjectRayDirectionEXT;
    vec3 n = normalize(hitPoint - center);
    hit.normal = n;
    hit.texCoord.s = fract(atan(n.y, n.x) / (2.0 * PI));
    hit.texCoord.t = acos(-n.z)  / (PI);
    reportIntersectionEXT(t, hitKind);
  }
}

• Next, we discuss the intersection of a ray with an Axis-Aligned Bounding Box (AABB)
• In contrast to an oriented bounding box, an AABB is aligned with the coordinate system, which simplifies the computation of the intersection point
• The so-called "slab" method is used (Kay and Kajiya: "Ray tracing complex scenes", SIGGRAPH 1986)
• A "slab" is an interval in one spatial direction that is confined by two parallel planes
• By analyzing the intervals, it can be determined whether and where an intersection occurs
• On the next slide, we will perform some initial experiments in 2D. The blue control points can be moved with the mouse.

• Parameterization of the ray:
  $\mathbf{x}(t) = \mathbf{s} + t \,\, \mathbf{r} \quad \forall\, \, t \in [t_{\mathrm{\small min}}, t_{\mathrm{\small max}}]$
  where $\mathbf{s}$ denotes the ray origin and $\mathbf{r}$ the normalized direction
• The AABB is spanned by two 3D points $\mathbf{a}$ and $\mathbf{b}$
  $\mathbf{a} = (a_x, a_y, a_z)^\top \quad \quad \mathbf{b} = (b_x, b_y, b_z)^\top$
• Inserting $\mathbf{a}$ and $\mathbf{b}$ into the ray equation results in 6 intersection points of the ray with the 6 planes of the AABB:
  $\begin{pmatrix}a_x\\a_y \\a_z\end{pmatrix}= \begin{pmatrix}s_x\\s_y \\s_z\end{pmatrix} + t \,\, \begin{pmatrix}r_x\\r_y \\r_z\end{pmatrix} \quad \quad \quad \begin{pmatrix}b_x\\b_y \\b_z\end{pmatrix}= \begin{pmatrix}s_x\\s_y \\s_z\end{pmatrix} + t \,\, \begin{pmatrix}r_x\\r_y \\r_z\end{pmatrix} $

• If the equations are solved for $t$, we get the distances along the ray for all 6 intersections with the AABB planes:
  $\begin{aligned} t_{a_x} = (a_x - s_x) / r_x & &t_{b_x} = (b_x - s_x) / r_x \\ t_{a_y} = (a_y - s_y) / r_y &&t_{b_y} = (b_y - s_y) / r_y \\ t_{a_z} = (a_z - s_z) / r_z &&t_{b_z} = (b_z - s_z) / r_z\\ \end{aligned}$
• For each dimension $x$, $y$, and $z$ there is an intersection interval:
  • x-interval: $t \in [t_{{\small\mbox{min}}_x}, \,t_{{\small\mbox{max}}_x}] = [\min(t_{a_x}, t_{b_x}), \max(t_{a_x}, t_{b_x})]$
  • y-interval: $t \in [t_{{\small\mbox{min}}_y}, \,t_{{\small\mbox{max}}_y}] = [\min(t_{a_y}, t_{b_y}), \max(t_{a_y}, t_{b_y})]$
  • z-interval: $t \in [t_{{\small\mbox{min}}_z}, \,t_{{\small\mbox{max}}_z}] = [\min(t_{a_z}, t_{b_z}), \max(t_{a_z}, t_{b_z})]$
• There are two potential intersections:
  • the maximum of the minimum values $t_{\mathrm{\small in}} = \max(t_{{\small\mbox{min}}_x}, \,t_{{\small\mbox{min}}_y}, \,t_{{\small\mbox{min}}_z})$
  • the minimum of the maximum values $t_{\mathrm{\small out}} = \min(t_{{\small\mbox{max}}_x}, \,t_{{\small\mbox{max}}_y}, \,t_{{\small\mbox{max}}_z})$

• There can be no, one, or two intersections
  • Two intersections if $t_{\mathrm{\small in}} < t_{\mathrm{\small out}}$:
    Entry at $t_{\mathrm{\small in}}$
    Exit at $t_{\mathrm{\small out}}$
  • One intersection if $t_{\mathrm{\small in}} = t_{\mathrm{\small out}}$
  • No intersection if $t_{\mathrm{\small in}} > t_{\mathrm{\small out}}$
• Typically the ray interval is limited
  $\mathbf{x}(t) = \mathbf{s} + t \,\, \mathbf{r} \quad \forall\, \, t \in [t_{\mathrm{\small min}}, t_{\mathrm{\small max}}]$.
  In this case, it must also be checked whether $t_{\mathrm{\small in}}$ or $t_{\mathrm{\small out}}$ are within the desired interval $[t_{\mathrm{\small min}}, t_{\mathrm{\small max}}]$

• Parameterization of the ray:
  $\mathbf{x}(t) = \mathbf{s} + t \,\, \mathbf{r} \quad \forall\, \, t \in [t_{\mathrm{\small min}}, t_{\mathrm{\small max}}]$
  where $\mathbf{s}$ denotes the ray origin and $\mathbf{r}$ the normalized direction
• The triangle is defined by three points $\mathbf{a}$, $\mathbf{b}$, and $\mathbf{c}$
  

  $\mathbf{a}$

  $\mathbf{b}$

  $\mathbf{c}$

  $\mathbf{s}$

  $\mathbf{r}$
• Plane equation:
  $\mathbf{x}(u, v) = \mathbf{a} + u \,(\mathbf{b} - \mathbf{a}) + v \,(\mathbf{c} - \mathbf{a})$
• Equating gives:
  $\mathbf{a} + u \,(\mathbf{b} - \mathbf{a}) + v \,(\mathbf{c} - \mathbf{a}) = \mathbf{s} + t \,\, \mathbf{r}$
• The equation system in matrix notation
  (3 equations, 3 unknowns):
  $\begin{bmatrix}-\mathbf{r} & \mathbf{b} - \mathbf{a} &\mathbf{c} - \mathbf{a}\end{bmatrix} \begin{pmatrix}t\\u\\v\end{pmatrix} = \mathbf{s} - \mathbf{a}$

• By substitution with $(\mathbf{b} - \mathbf{a}) = \mathbf{i}$, $(\mathbf{c} - \mathbf{a}) = \mathbf{j}$, and $(\mathbf{s} - \mathbf{a}) = \mathbf{k}$:
  $\begin{bmatrix}-\mathbf{r} & \mathbf{b} - \mathbf{a} &\mathbf{c} - \mathbf{a}\end{bmatrix} \begin{pmatrix}t\\u\\v\end{pmatrix} = \mathbf{s} - \mathbf{a} \quad\Leftrightarrow \quad \begin{bmatrix}-\mathbf{r} & \mathbf{i} &\mathbf{j}\end{bmatrix} \begin{pmatrix}t\\u\\v\end{pmatrix} = \mathbf{k} $
• The solution of the linear system of equations is given by (Cramer's rule)
  $\begin{pmatrix}t\\u\\v\end{pmatrix} = \frac{1}{\det(-\mathbf{r}, \mathbf{i}, \mathbf{j})} \begin{bmatrix}\det(\mathbf{k}, \mathbf{i}, \mathbf{j})\\ \det(-\mathbf{r}, \mathbf{k}, \mathbf{j})\\ \det(-\mathbf{r}, \mathbf{i}, \mathbf{k})\\ \end{bmatrix} = \frac{1}{(\mathbf{r} \times \mathbf{j})^\top \mathbf{i}} \begin{bmatrix}(\mathbf{k} \times \mathbf{i})^\top \mathbf{j}\\ (\mathbf{r} \times \mathbf{j})^\top \mathbf{k}\\ (\mathbf{k} \times \mathbf{i})^\top \mathbf{r}\\ \end{bmatrix}$
  with $\det(\mathbf{x}, \mathbf{y}, \mathbf{z}) = (\mathbf{x} \times \mathbf{y})^\top \mathbf{z} = -(\mathbf{x} \times \mathbf{z})^\top \mathbf{y} = -(\mathbf{z} \times \mathbf{y})^\top \mathbf{x}$
• The solution requires only 2 cross products and 4 scalar products

$\mathbf{a}$

$\mathbf{b}$

$\mathbf{c}$

• The plane equation from before
  $\mathbf{x}(u, v) = \mathbf{a} + u \,(\mathbf{b} - \mathbf{a}) + v \,(\mathbf{c} - \mathbf{a})$
  can also be written as
  $\mathbf{x}(u, v) = (1- u - v) \,\mathbf{a} + u \,\mathbf{b} + v \,\mathbf{c}$
• Now, an additional variable is introduced $w = (1 - u - v)$ and the plane equation reads:
  $\mathbf{x}(w, u, v) = w \,\mathbf{a} + u \,\mathbf{b} + v \,\mathbf{c}$
  with the constraint $u + v + w = 1$
• The plane is now over-parameterized (3 instead of 2 parameters) but the constraint ensures that all points $\mathbf{x}(w, u, v)$ are on the plane
• The parameters $w$, $u$, $v$ are called barycentric coordinates

$\mathbf{a}$

$\mathbf{b}$

$\mathbf{c}$

• The barycentric coordinates $u$, $v$, and $w$ are the coefficients for the linear combination of the triangle's vertices:
  $\mathbf{x}(w, u, v) = w \,\mathbf{a} + u \,\mathbf{b} + v \,\mathbf{c}$
  with the constraint $u + v + w = 1$
• If $u$, $v$, and $w \ge 0$, the point $\mathbf{x}(w, u, v)$ lies within the triangle

$\mathbf{x}(w=\frac{1}{3}, u=\frac{1}{3}, v=\frac{1}{3})$

$\mathbf{x}(w=0, u=\frac{1}{2}, v=\frac{1}{2})$

$\mathbf{x}(w=0, u=0, v=1)$

$\mathbf{x}$

$\mathbf{x}$

$\mathbf{x}$

• Barycentric coordinates are often used to interpolate vertex data (such as color, texture coordinate, normal, etc.)
• The built-in intersection shader for triangles also provides barycentric coordinates
  $u$ und $v$ of the intersection point via the
  "hitAttributeEXT vec2 baryCoord" variable
• The $w$ coordinate can be computed with:
  $w = (1 - u - v)$

vec3 barys = vec3(1.0f - baryCoord.x - baryCoord.y, baryCoord.x, baryCoord.y);
vec3 interpColor = blue * barys.x + red * barys.y + green * barys.z;

• Intersection computations for more complex objects can be composed of intersection calculations for simple basic shapes (primitives)
• The combination is typically performed by the following operations:
  • Union $\mathcal{A}\cup \mathcal{B}=\{x\in \mathcal{A} \lor x\in \mathcal{B}\}$
  • Intersection $\mathcal{A} \cap \mathcal{B} =\{x\in \mathcal{A} \land x\in \mathcal{B}\}$
  • Difference $\mathcal{A}\setminus \mathcal{B}=\{x\in \mathcal{A}\land x\notin \mathcal{B}\}$

Union

Intersection

Difference

$\mathbf{s}$

$\mathbf{r}$

$t_{\mathrm{\small in}}$

$t_{\mathrm{\small out}}$

• In order to implement the operations in an intersection shader, the intervals for the ray parameter $t$ are considered, in which $t$ is within the object
• For simple (convex) objects (sphere, cube) this interval can be described by the entry and exit point $[t_{\mathrm{\small in}}, t_{\mathrm{\small out}}]$
  

  $\mathcal{A}$

  $\mathcal{B}$

  $t_{\mathcal{A}, \mathrm{\small in}}$

  $t_{\mathcal{A}, \mathrm{\small out}}$

  $t_{\mathcal{B}, \mathrm{\small in}}$

  $t_{\mathcal{B}, \mathrm{\small out}}$
• Now we can think about how the new intervals are computed after the various CSG operations
• For this purpose, 5 different cases are considered (see next slide)