How to Test the Ray-Triangle Intersection in C++

Syed Hassan Sabeeh Kazmi Feb 02, 2024

C++ C++ Intersection

Normalize the Ray Direction to Test the Ray-Triangle Intersection in C++
Write an Algorithm to Test a Ray-Triangle Intersection in C++

How to Test the Ray-Triangle Intersection in C++

Testing the ray-triangle intersection could require millions of tests and is known to be one of the kernel operations in any ray tracer (requires different function implementation for each geometric primitive type) in C++. This tutorial will teach you how to program ray-triangle intersections in C++.

Triangles, as the rendering primitives, are exceptionally unique and act as one of the rarest denominators between modelers/renderers. Their unique characteristics make them simple, uniform, and easy to tessellate, and for these reasons, C++ use them.

You will learn two primary methods for programming and optimizing ray-triangle intersection. Both are similar, but the second method/approach is more practical than the other as it supports modern compilers and processors.

Normalize the Ray Direction to Test the Ray-Triangle Intersection in C++

It’s possible to normalize the ray direction using a single variable similar to the triangle variables, as the triangle’s vertices are found in an object representing a triangle. It’s great if you are familiar with the Moller-Trumbore ray-triangle algorithm and can use it only if you can handle the complexity.

Always normalize the ray direction and triangle elements before doing cross-product to achieve/get consistent results. To have the collision between ray and triangle, you need to have something like dot((point - planeOrigin), planeNormal) > 0.0 for all the three planes created by triangle sides.

#include <locale.h>
#include <math.h>
#include <stdint.h>
#include <stdio.h>
#include <stdlib.h>
#include <sys/time.h>
#include <time.h>

// vector mathematics declaration for algorithm
struct raytri_intersec {
  float vector_1;
  float vector_2;
  float vector_3;
};

// ray structure for initializing its origin and direction
struct Ray {
  raytri_intersec ray_origin;
  raytri_intersec ray_direction;
};

raytri_intersec intersec_sub(raytri_intersec point_a, raytri_intersec point_b) {
  return {point_a.vector_1 - point_b.vector_1,
          point_a.vector_2 - point_b.vector_2,
          point_a.vector_3 - point_b.vector_3};
}

float intersec_dot(raytri_intersec point_a, raytri_intersec point_b) {
  return point_a.vector_1 * point_b.vector_1 +
         point_a.vector_2 * point_b.vector_2 +
         point_a.vector_3 * point_b.vector_3;
}

float intersec_len(raytri_intersec obj_vec) {
  return sqrt(obj_vec.vector_1 * obj_vec.vector_1 +
              obj_vec.vector_2 * obj_vec.vector_2 +
              obj_vec.vector_3 * obj_vec.vector_3);
}

raytri_intersec intersec_normalize(raytri_intersec obj_vec) {
  float length = intersec_len(obj_vec);
  return {obj_vec.vector_1 / length, obj_vec.vector_2 / length,
          obj_vec.vector_3 / length};
}

raytri_intersec intersec_cross(raytri_intersec point_a,
                               raytri_intersec point_b) {
  return {
      point_a.vector_2 * point_b.vector_3 - point_a.vector_3 * point_b.vector_2,
      point_a.vector_3 * point_b.vector_1 - point_a.vector_1 * point_b.vector_3,
      point_a.vector_1 * point_b.vector_2 -
          point_a.vector_2 * point_b.vector_1};
}

// ray-triangle intersection routine
float raytriangle_intersection(Ray *_ray, raytri_intersec *_vector1,
                               raytri_intersec *_vector2,
                               raytri_intersec *_vector3) {
  raytri_intersec intersec_vec1vec2 = intersec_sub(*_vector2, *_vector1);
  raytri_intersec intersec_vec1vec3 = intersec_sub(*_vector3, *_vector1);

  raytri_intersec intersec_pev =
      intersec_cross(_ray->ray_direction, intersec_vec1vec3);

  float intersec_det = intersec_dot(intersec_vec1vec2, intersec_pev);

  if (intersec_det < 0.000001) return -INFINITY;

  float intersec_invDet = 1.0 / intersec_det;

  raytri_intersec intersec_tvec = intersec_sub(_ray->ray_origin, *_vector1);

  float random_obj =
      intersec_dot(intersec_tvec, intersec_pev) * intersec_invDet;

  if (random_obj < 0 || random_obj > 1) return -INFINITY;

  raytri_intersec intersec_qvec =
      intersec_cross(intersec_tvec, intersec_vec1vec2);

  float random_objec =
      intersec_dot(_ray->ray_direction, intersec_qvec) * intersec_invDet;

  if (random_objec < 0 || random_obj + random_objec > 1) return -INFINITY;

  return intersec_dot(intersec_vec1vec3, intersec_qvec) * intersec_invDet;
}

/* Test data generation */

raytri_intersec *allocate_triangle(int num_triangles) {
  return (raytri_intersec *)malloc(sizeof(raytri_intersec) * num_triangles * 3);
}

raytri_intersec intersec_random_sphere() {
  double ray_1 = (float)rand() / RAND_MAX;
  double ray_2 = (float)rand() / RAND_MAX;
  double _latitude = acos(2 * ray_1 - 1) - M_PI / 2;
  double _longitude = 2 * M_PI * ray_2;

  return {(float)(cos(_latitude) * cos(_longitude)),
          (float)(cos(_latitude) * sin(_longitude)), (float)sin(_latitude)};
}

raytri_intersec *generate_randomtriangles(int num_triangles) {
  raytri_intersec *_vertices = allocate_triangle(num_triangles);

  for (int i = 0; i < num_triangles; ++i) {
    _vertices[i * 3 + 0] = intersec_random_sphere();
    _vertices[i * 3 + 1] = intersec_random_sphere();
    _vertices[i * 3 + 2] = intersec_random_sphere();
  }

  return _vertices;
}

long intersec_ellapsed_Ms(struct timeval time_sec, struct timeval time_usec) {
  return 1000 * (time_usec.tv_sec - time_sec.tv_sec) +
         (time_usec.tv_usec - time_sec.tv_usec) / 1000;
}

int main() {
  const int number_rays = 1000;
  const int number_triangles = 100 * 1000;

  srand(time(NULL));

  raytri_intersec *_vertices = generate_randomtriangles(number_triangles);
  const int number_vertices = number_triangles * 3;

  int intersec_hit = 0;
  int intersec_miss = 0;

  Ray _ray;

  struct timeval time_start;
  gettimeofday(&time_start, 0);

  for (int ing = 0; ing < number_rays; ++ing) {
    _ray.ray_origin = intersec_random_sphere();
    raytri_intersec point_one = intersec_random_sphere();
    _ray.ray_direction =
        intersec_normalize((intersec_sub(point_one, _ray.ray_origin)));

    for (int temp = 0; temp < number_vertices / 3; ++temp) {
      float _time = raytriangle_intersection(&_ray, &_vertices[temp * 3 + 0],
                                             &_vertices[temp * 3 + 1],
                                             &_vertices[temp * 3 + 2]);
      _time >= 0 ? ++intersec_hit : ++intersec_miss;
    }
  }

  struct timeval time_end;
  gettimeofday(&time_end, 0);

  double time_total =
      (float)intersec_ellapsed_Ms(time_start, time_end) / 1000.0;

  int numberof_tests = number_rays * number_triangles;
  float intersec_hitpercentage =
      ((float)intersec_hit / numberof_tests) * 100.0f;
  float intersec_misspercentage =
      ((float)intersec_miss / numberof_tests) * 100.0f;
  float tests_persecond = (float)numberof_tests / time_total / 1000000.0f;

  setlocale(LC_NUMERIC, "");

  printf("Total intersection tests:  %'10i\n", numberof_tests);
  printf("  Hits:\t\t\t    %'10i (%5.2f%%)\n", intersec_hit,
         intersec_hitpercentage);
  printf("  Misses:\t\t    %'10i (%5.2f%%)\n", intersec_miss,
         intersec_misspercentage);
  printf("\n");
  printf("Total time:\t\t\t%6.2f seconds\n", time_total);
  printf("Millions of tests per second:\t%6.2f\n", tests_persecond);
}

Output:

Total intersection tests:  100,000,000
  Hits:			     4,819,136 ( 4.82%)
  Misses:		    95,180,864 (95.18%)

Total time:			  8.52 seconds
Millions of tests per second:	 11.74

Furthermore, this approach is a platform for the next method. It uses the knowledge and technique from the first one and further improves it to enhance performance and enable in-detail customization.

When studying the ray-triangle intersection in programming, it’s important to understand the ray and triangle positions.

The ray and the triangle can be parallel, the triangle can be behind the ray, or P represented as the point of intersection can either be inside or outside of the triangle. Studying the case and geometric calculations can lead you to program a write algorithm.

Write an Algorithm to Test a Ray-Triangle Intersection in C++

The following C++ implementation of the ray-triangle intersection is tailored for optimum performance. To ensure numerical stability, we need the test code to eliminate parallel rays and must compare the determinant to a small interval around 0.

This algorithm will also reflect the inside-outside technique of ray-triangle intersection, enabling it for any convex polygon. The resulting sign of the dot product of the resulting vector (as an output) and the convex polygon’s normal holds the information to determine the point of the vector’s edge.

// IMPORTANT [NOTE]
// Include the `geometry.h` header file to your project before running this
// program the `geometry.h` header file link -
// https://drive.google.com/file/d/1qoDbXA_HnhbtS5KeyJxZT8FfwtVAJeO5/view?usp=sharing

// for windows users
#define _USE_MATH_DEFINES

// general header declaration
#include <cmath>
#include <cstdint>
#include <cstdio>
#include <cstdlib>
#include <fstream>
#include <iostream>
#include <limits>
#include <memory>
#include <random>
#include <utility>
#include <vector>

// load the `geometry.h` header file to your project
#include "geometry.h"

constexpr float constant_exprator = 1e-8;

inline float scale_degtorag(const float &current_deg2rad) {
  return current_deg2rad * M_PI / 180;
}

inline float clamp(const float &val_low, const float &val_high,
                   const float &value_current) {
  return std::max(val_low, std::min(val_high, value_current));
}
bool rayTriangleIntersect(const Vec3f &vector_original,
                          const Vec3f &vector_direct, const Vec3f &vector_zero,
                          const Vec3f &vector_one, const Vec3f &vector_two,
                          float &float_val1, float &float_val2,
                          float &float_val3) {
// The `Moller Trumbore` algorithm (to study ray-triangle intersection)
// implementation in C++
#ifdef MOLLER_TRUMBORE
  Vec3f v0v1 = v1 - v0;
  Vec3f v0v2 = v2 - v0;
  Vec3f pvec = dir.crossProduct(v0v2);
  float det = v0v1.dotProduct(pvec);
#ifdef CULLING
  // Important | The triangle will be back-facing in case of a `-ve` (negative)
  // determinant, and if it is close to 0 (zero), it means the ray misses the
  // triangle
  if (det < kEpsilon) return false;
#else
  // Important | The `det` is close to 0 (zero) means the ray and triangle are
  // parallel to each other
  if (fabs(det) < kEpsilon) return false;
#endif
  float invDet = 1 / det;

  Vec3f tvec = orig - v0;
  u = tvec.dotProduct(pvec) * invDet;
  if (u < 0 || u > 1) return false;

  Vec3f qvec = tvec.crossProduct(v0v1);
  v = dir.dotProduct(qvec) * invDet;
  if (v < 0 || u + v > 1) return false;

  t = v0v2.dotProduct(qvec) * invDet;

  return true;
#else
  // in case the compute plane is normal
  Vec3f vector_V1 = vector_one - vector_zero;
  Vec3f vector_V2 = vector_two - vector_zero;
  // it does not require any normalization
  Vec3f crosspro_vector = vector_V1.crossProduct(vector_V2);  // N
  float vector_denom = crosspro_vector.dotProduct(crosspro_vector);

  // Step 1: finding the `P` point of intersection and check if the ray and the
  // plane are parallel
  float dot_raydirection_inter = crosspro_vector.dotProduct(vector_direct);

  // can almost 0 (zero) and return false in case of parallel, which means the
  // ray and triangle do not intersect
  if (fabs(dot_raydirection_inter) < constant_exprator) return false;

  // use the `equation 2` to compute the `direction_vec` parameter
  float direction_vec = -crosspro_vector.dotProduct(vector_zero);

  // use the `equation 3` to compute the `float_val1` parameter
  float_val1 = -(crosspro_vector.dotProduct(vector_original) + direction_vec) /
               dot_raydirection_inter;

  // check if the triangle is behind the ray, and in case it returns `false`, it
  // means `yes`
  if (float_val1 < 0) return false;

  // compute the ray-triangle intersection point `intersection_point` using the
  // `equation 1`
  Vec3f intersection_point = vector_original + float_val1 * vector_direct;

  // Step 2: perform the inside-out test by checking the vector perpendicular to
  // the triangle's plane
  Vec3f insideout_test;

  // zero_edge
  Vec3f zero_edge = vector_one - vector_zero;
  Vec3f vector_point_zero = intersection_point - vector_zero;
  insideout_test = zero_edge.crossProduct(vector_point_zero);
  if (crosspro_vector.dotProduct(insideout_test) < 0)
    return false;  // it returning false means the `insideout_test` point is on
                   // the right side

  // single_edge
  Vec3f single_edge = vector_two - vector_one;
  Vec3f vector_point_single = intersection_point - vector_one;
  insideout_test = single_edge.crossProduct(vector_point_single);
  if ((float_val2 = crosspro_vector.dotProduct(insideout_test)) < 0)
    return false;  // it returning false means the `insideout_test` point is on
                   // the right side

  // double_edge
  Vec3f double_edge = vector_zero - vector_two;
  Vec3f vector_point_double = intersection_point - vector_two;
  insideout_test = double_edge.crossProduct(vector_point_double);
  if ((float_val3 = crosspro_vector.dotProduct(insideout_test)) < 0)
    return false;  // it returning false means the `insideout_test` point is on
                   // the right side

  float_val2 /= vector_denom;
  float_val3 /= vector_denom;

  return true;  // means this ray hits the triangle
#endif
}

int main(int argc, char **argv) {
  Vec3f first_vec(-1, -1, -5);
  Vec3f second_vec(1, -1, -5);
  Vec3f third_vec(0, 1, -5);

  const uint32_t visible_width = 640;
  const uint32_t visible_height = 480;
  Vec3f inter_columns[3] = {{0.6, 0.4, 0.1}, {0.1, 0.5, 0.3}, {0.1, 0.3, 0.7}};
  Vec3f *intersecframe_buffer = new Vec3f[visible_width * visible_height];
  Vec3f *intersecPoint_pixel = intersecframe_buffer;
  float field_of_view = 51.52;
  float intersec_scale = tan(scale_degtorag(field_of_view * 0.5));
  float visible_image_ratio = visible_width / (float)visible_height;
  Vec3f ray_origin(0);
  for (uint32_t temp = 0; temp < visible_height; ++temp) {
    for (uint32_t ing = 0; ing < visible_width; ++ing) {
      // compute primary ray
      float x_axis = (2 * (ing + 0.5) / (float)visible_width - 1) *
                     visible_image_ratio * intersec_scale;
      float y_axis =
          (1 - 2 * (temp + 0.5) / (float)visible_height) * intersec_scale;
      Vec3f intersec_direction(x_axis, y_axis, -1);
      intersec_direction.normalize();
      float a, b, c;
      if (rayTriangleIntersect(ray_origin, intersec_direction, first_vec,
                               second_vec, third_vec, a, b, c)) {
        *intersecPoint_pixel = b * inter_columns[0] + c * inter_columns[1] +
                               (1 - b - c) * inter_columns[2];

        // comment the following line of code if you do not want to visualize
        // the row bary-centric coordinates
        *intersecPoint_pixel = Vec3f(b, c, 1 - b - c);
      }
      intersecPoint_pixel++;
    }
  }

  // Save the result to a PPM image (keep these flags if you compile under
  // Windows)
  std::ofstream obj_ofstream("./out.ppm", std::ios::out | std::ios::binary);
  obj_ofstream << "P6\n" << visible_width << " " << visible_height << "\n255\n";

  std::cout << "The ray-triangle intersection result: " << intersecPoint_pixel
            << "\n"
            << "Width | " << visible_width << "\nHeight | " << visible_height;

  for (uint32_t imp = 0; imp < visible_height * visible_width; ++imp) {
    char red = (char)(255 * clamp(0, 1, intersecframe_buffer[imp].x));
    char green = (char)(255 * clamp(0, 1, intersecframe_buffer[imp].y));
    char blue = (char)(255 * clamp(0, 1, intersecframe_buffer[imp].z));
    obj_ofstream << red << green << blue;
  }

  obj_ofstream.close();

  delete[] intersecframe_buffer;

  return 0;
}

Output:

The ray-triangle intersection result: 0x1d1cdb5d040
Width | 640
Height | 480

Writing an algorithm based on the Moller Trumbore technique enables no storage requirement for the plane equation, which can save significant memory for triangle meshes. Furthermore, using this technique, you will get the fastest ray-triangle insertion routine without having complex or pre-computed plane equations.

Author: Syed Hassan Sabeeh Kazmi

Hassan is a Software Engineer with a well-developed set of programming skills. He uses his knowledge and writing capabilities to produce interesting-to-read technical articles.

GitHub