C++ で光線と三角形の交差をテストする
Syed Hassan Sabeeh Kazmi
2023年10月12日
光線と三角形の交差をテストするには、何百万ものテストが必要になる可能性があり、C++ の任意の光線追跡のカーネル操作の 1つとして知られています (ジオメトリ プリミティブ タイプごとに異なる関数の実装が必要です)。 このチュートリアルでは、C++ で光線と三角形の交差をプログラミングする方法を説明します。
レンダリング プリミティブとしての三角形は非常にユニークで、モデラーとレンダーの間で最もまれな分母の 1つとして機能します。 それらの独自の特性により、それらはシンプルで統一され、テッセレーションが容易になります。これらの理由から、C++ はそれらを使用します。
光線と三角形の交差をプログラミングおよび最適化するための 2つの主要な方法を学習します。 どちらも似ていますが、最新のコンパイラとプロセッサをサポートしているため、2 番目の方法/アプローチは他の方法よりも実用的です。
光線の方向を正規化して、C++ で光線と三角形の交差をテストする
三角形の頂点は三角形を表すオブジェクト内にあるため、三角形の変数と同様の単一の変数を使用して光線の方向を正規化することができます。 Moller-Trumbore ray-triangle アルゴリズムに精通していて、複雑さを処理できる場合にのみそれを使用できるのは素晴らしいことです。
一貫した結果を達成/取得するために、外積を実行する前に、常に光線方向と三角形要素を正規化してください。 光線と三角形を衝突させるには、三角形の側面によって作成された 3つの平面すべてに対して dot((point - planeOrigin), planeNormal) > 0.0 のようなものが必要です。
#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
さらに、このアプローチは次の方法のプラットフォームです。 最初の知識と技術を使用し、さらに改善してパフォーマンスを向上させ、詳細なカスタマイズを可能にします。
プログラミングで光線と三角形の交差を調べる場合、光線と三角形の位置を理解することが重要です。
光線と三角形は平行にすることも、三角形を光線の後ろに置くことも、交点として表される 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 <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 ¤t_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 手法に基づくアルゴリズムを作成すると、平面方程式のストレージ要件がなくなり、三角形メッシュのメモリを大幅に節約できます。 さらに、この手法を使用すると、複雑な平面方程式や事前に計算された平面方程式を使用せずに、最速の光線三角形挿入ルーチンを取得できます。
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