迹忆客 专注技术分享

测试光线-三角形相交可能需要数百万次测试，并且被认为是 C++ 中任何光线追踪器的核心操作之一（需要为每种几何基元类型提供不同的函数实现）。 本教程将教您如何用 C++ 编写射线-三角形相交的程序。

三角形作为渲染基元非常独特，是建模者/渲染者之间最稀有的分母之一。 它们独特的特性使它们简单、统一且易于细分，因此 C++ 使用它们。

您将学习编程和优化射线三角形相交的两种主要方法。 两者很相似，但第二种方法比另一种更实用，因为它支持现代编译器和处理器。

在 C++ 中标准化光线方向以测试光线与三角形的相交 可以使用类似于三角形变量的单个变量来标准化光线方向，因为三角形的顶点是在表示三角形的对象中找到的。 如果您熟悉 Moller-Trumbore 射线三角形算法并且只有在能够处理复杂性的情况下才能使用它，那就太好了。

在进行叉积之前，请始终对光线方向和三角形元素进行归一化，以实现/获得一致的结果。 要使射线和三角形之间发生碰撞，对于三角形边创建的所有三个平面，需要有类似 dot((point - planOrigin),planeNormal) > 0.0 的值。

#include <stdio.h>
#include <math.h>
#include <time.h>
#include <locale.h>
#include <stdlib.h>
#include <stdint.h>
#include <sys/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

此外，这种方法是下一种方法的平台。 它使用第一个知识和技术，并进一步改进它以提高性能并实现详细的定制。

在编程中研究射线与三角形相交时，了解射线和三角形的位置非常重要。

射线和三角形可以平行，三角形可以在射线后面，或者表示为交点的 P 可以在三角形内部或外部。 研究案例和几何计算可以引导您编写写入算法。

用 C++ 编写测试射线三角形相交的算法

下面的射线-三角形相交的 C++ 实现是为了获得最佳性能而定制的。 为了确保数值稳定性，我们需要测试代码来消除平行光线，并且必须将行列式与 0 附近的小区间进行比较。

该算法还将反映射线-三角形相交的内-外技术，使其适用于任何凸多边形。 结果向量（作为输出）和凸多边形法线的点积的结果符号包含确定向量边缘点的信息。

// 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 <cstdio>
#include <cstdlib>
#include <fstream>
#include <iostream>
#include <limits>
#include <memory>
#include <random>
#include <utility>
#include <vector>
#include <cstdint>

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

编写基于 Moller Trumbore 技术的算法无需存储平面方程，这可以为三角形网格节省大量内存。 此外，使用这种技术，您将获得最快的射线三角形插入例程，而无需复杂或预先计算的平面方程。