Ray

编程入门行业动态更新时间:2024-10-22 11:21:58

Ray

思路概述

射线和三角形求交的思路很简单，求出射线和三角形所在的平面的交点，然后算这个交点的在三角形的重心坐标，重心坐标在[0, 1]之间，说明射线和三角形相交，否则不相交。

假设我么的ray的方程是

这个计算总体来说是两步，第一步求射线和平面的交点，只需要求出t。第二步是求三角中心坐标u和v，我们想找出一个更高效的方法，叫unit triangle intersection method。

简单来说，就是把三角形的三个点变换到(0, 0, 0)，(1, 0, 0)，(0, 1, 0)这三个点构成的空间里，然后射线和xoy平面相交即可，如下图所示：

注意，这个图是右手坐标。

这里用一个Affine transform来做变换，Affine transform的意思是，直线变换后仍然是直线，如果变换后是曲线，就不是Affine transform了，所以我们日常用的平移，缩放，旋转，都是Affine transform。

公式推导

我们设定这个Affine transform的变换过程是T(P)，P是世界坐标中一点。

三角形的三个点是ABC，那么有：

，，，

那么T的逆变换就是：

所以我们可以构造T变换的逆变换为：

最好我们可以用矩阵的形式去表示这个逆变换：

那么变换到Unit triangle空间中的Ray就是

由于射线和xoy平面相交，也就是交点的z坐标为0.

所以t可以很方便的计算：

下面看三角形的重心坐标u v如何计算。

交点在xoy平面上的重心坐标可以表示为：

由于C'就是原点0，A' = (1, 0, 0)，B' = (0, 1, 0)，那么展开向量：

最后得到的uv是：

整个过程计算完毕，是不是很简单呢！

代码实现

第一步：构建矩阵M

Matrix4x4 matrix = new Matrix4x4();
Vector4 col0 = v0 - v2;
col0.w = 0;
matrix.SetColumn(0, col0);
Vector4 col1 = v1 - v2;
col1.w = 0;
matrix.SetColumn(1, col1);
Vector4 col2 = Vector3.Cross(v0 - v2, v1 - v2);
col2.w = 0;
matrix.SetColumn(2, col2);
Vector4 col3 = v2;
col3.w = 1;
matrix.SetColumn(3, col3);
matrix = Matrix4x4.Inverse(matrix);

第二步，把矩阵写入buffer中，这里我不贴代码了，简单来说就是compute shader或者cuda里在求交时要访问的bvh叶子节点。

第三步，就是求交的代码，假设我把矩阵M的每一行传递到buffer里。

我们要求的是O'和D'，那么根据矩阵和向量的乘法，求t的代码如下：

const float4 m0 = Positions[triAddr];     //matrix row 0 
const float4 m1 = Positions[triAddr + 1]; //matrix row 1 
const float4 m2 = Positions[triAddr + 2]; //matrix row 2//Oz is a point, must plus w
float Oz = m2.w + origx * m2.x + origy * m2.y + origz * m2.z;
//Dz is a vector
float invDz = 1.0f / (dirx * m2.x + diry * m2.y + dirz * m2.z);
float t = -Oz * invDz;

下面是求u和v的代码。

//if t is in bounding and less than the ray.tMax
if (t > tmin && t < hitT)
{// Compute and check barycentric u.float Ox = m0.w + origx * m0.x + origy * m0.y + origz * m0.z;float Dx = dirx * m0.x + diry * m0.y + dirz * m0.z;float u = Ox + t * Dx;if (u >= 0.0f){// Compute and check barycentric v.float Oy = m1.w + origx * m1.x + origy * m1.y + origz * m1.z;float Dy = dirx * m1.x + diry * m1.y + dirz * m1.z;float v = Oy + t * Dy;if (v >= 0.0f && u + v <= 1.0f){// Record intersection.// Closest intersection not required => terminate.hitT = t;hitIndex = triAddr;}}
}