Probar la intersección Ray-Triangle en C++

#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);
}

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

// 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;
}

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