Ray-Triangle intersection

思路概述

射线和三角形求交的思路很简单，求出射线和三角形所在的平面的交点，然后算这个交点的在三角形的重心坐标，重心坐标在[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;
	    }
    }
}

参考

http://www.sven-woop.de/papers/2004-GH-SaarCOR.pdf