Skip to content

Instantly share code, notes, and snippets.

@kennyalive

kennyalive/intersect_triangle_watertight.cpp

Last active November 6, 2021 17:55
Show Gist options

  • Select an option

    Learn more about clone URLs
  • Save kennyalive/c3cc15e2208fdca8a72b9fce4b9874a8 to your computer and use it in GitHub Desktop.

Select an option

Learn more about clone URLs
Save kennyalive/c3cc15e2208fdca8a72b9fce4b9874a8 to your computer and use it in GitHub Desktop.
Download ZIP
watertight triangle intersection with ray origin adjustment
Raw
intersect_triangle_watertight.cpp
constexpr bool perform_expensive_numerical_analysis_to_prevent_intersection_behind_ray_origin = true; // could noticeable slowdown algorithm (like 2x)
//
// Sven Woop, Carsten Benthin, and Ingo Wald, Watertight Ray/Triangle Intersection,
// Journal of Computer Graphics Techniques (JCGT), vol. 2, no. 1, 65-82, 2013
// http://jcgt.org/published/0002/01/05/
//
float intersect_triangle_watertight(const Ray& ray, const Vector3& p0, const Vector3& p1, const Vector3& p2, Vector3* barycentrics)
{
const int kz = ray.direction.abs().max_dimension();
const int kx = (kz == 2 ? 0 : kz + 1);
const int ky = (kz == 0 ? 2 : kz - 1);
Vector3 direction = ray.direction.permutation(kx, ky, kz);
float sx = -direction.x / direction.z;
float sy = -direction.y / direction.z;
float sz = 1.f / direction.z;
const Vector3 p0t = (p0 - ray.origin).permutation(kx, ky, kz);
const Vector3 p1t = (p1 - ray.origin).permutation(kx, ky, kz);
const Vector3 p2t = (p2 - ray.origin).permutation(kx, ky, kz);
float x0 = p0t.x + sx * p0t.z;
float y0 = p0t.y + sy * p0t.z;
float x1 = p1t.x + sx * p1t.z;
float y1 = p1t.y + sy * p1t.z;
float x2 = p2t.x + sx * p2t.z;
float y2 = p2t.y + sy * p2t.z;
float e0 = x1*y2 - y1*x2;
float e1 = x2*y0 - y2*x0;
float e2 = x0*y1 - y0*x1;
if (e0 == 0.f || e1 == 0.f || e2 == 0.f) {
e0 = float(double(x1)*double(y2) - double(y1)*double(x2));
e1 = float(double(x2)*double(y0) - double(y2)*double(x0));
e2 = float(double(x0)*double(y1) - double(y0)*double(x1));
}
if ((e0 < 0 || e1 < 0 || e2 < 0) && (e0 > 0 || e1 > 0 || e2 > 0))
return Infinity;
float det = e0 + e1 + e2;
if (det == 0.f)
return Infinity;
float z0 = sz * p0t.z;
float z1 = sz * p1t.z;
float z2 = sz * p2t.z;
float t_scaled = e0*z0 + e1*z1 + e2*z2;
// This following code is more efficient version of this snippet:
// if (det < 0 && t_scaled > 0 || det > 0 && t_scaled < 0) return Infinity;
// MSVC compiler generates 2 x64 instructions: xor, jl.
// Benchmark consistenly shows 1 cycle gain :)
// TODO: replace memcpy with std::bit_cast if no surprises in generated code
{
uint32_t det_bits;
memcpy(&det_bits, &det, 4);
uint32_t t_scaled_bits;
memcpy(&t_scaled_bits, &t_scaled, 4);
if ((det_bits ^ t_scaled_bits) >> 31)
return Infinity;
}
float inv_det = 1.f / det;
float t = inv_det * t_scaled;
ASSERT(t >= 0);
if constexpr(perform_expensive_numerical_analysis_to_prevent_intersection_behind_ray_origin)
{
float max_xt = Vector3(p0t.x, p1t.x, p2t.x).abs().max_component();
float max_yt = Vector3(p0t.y, p1t.y, p2t.y).abs().max_component();
float max_zt = Vector3(p0t.z, p1t.z, p2t.z).abs().max_component();
float max_x = Vector3(x0, x1, x2).abs().max_component();
float max_y = Vector3(y0, y1, y2).abs().max_component();
float max_z = Vector3(z0, z1, z2).abs().max_component();
float max_e = Vector3(e0, e1, e2).abs().max_component();
float delta_x = gamma(4) * (max_xt + max_zt);
float delta_y = gamma(4) * (max_yt + max_zt);
float delta_z = gamma(3) * max_z;
float delta_xy = gamma(1)*max_x*max_y + delta_x*max_y + delta_y*max_x;
float delta_e = 2*delta_xy + gamma(1)*max_e;
float delta_ez = gamma(1)*max_e*max_z + delta_z*max_e + delta_e*max_z;
float delta_t = (3*delta_ez + gamma(2)*std::abs(t_scaled)) * std::abs(inv_det);
if (delta_t >= t)
return Infinity;
}
barycentrics->x = e0 * inv_det;
barycentrics->y = e1 * inv_det;
barycentrics->z = e2 * inv_det;
return t;
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment