2020-05-27计算机图形学 0次访问

Ray-tracing (3)—Intersection

ps: Thumbnail 《月色真美》。

在Ray-tracing (1)—Overview文章中，我们已经知道我们需要遍历场景中所有的三角面找出和光线相交的那一个，那么本文将讨论的是如何进行三角形的求交从而找到交点。

三角形求交的一般方法

我们可以把三角形的求交分成两个步骤:

求光线和三角形所在平面的交点;
测试交点是否在三角形平面内。

光线和平面求交

上图的平面方程:
$$(p - p\prime) \cdot N = 0.$$$p$代表平面上任意一点，$p\prime$是屏幕上已知的一点，$N$是平面法线。

光线方程(在Ray-tracing (2)—Generating Camera Rays文章中有介绍):
$$r(t) = o + td, 0\leq t < \infty$$o代表的是光线的原点,d代表的是光线的方向,t代表的是沿着光线方向的距离。
把光线方程代入平面方程有:
$$(p - p\prime) \cdot N = (o + td - p\prime) \cdot N = 0$$所以可得:$$t = \frac{(p\prime - o) \cdot N}{d \cdot N}, 0\leq t < \infty$$解出来t之后就可以根据光线方程得到交点的坐标了。

交点是否在三角形内

上面只是求得了光线和三角形所在平面的交点，但是交点在不在三角形中还得讨论一下。

如上图所示，分别让三角形三边的向量和三点连向Q的向量做叉乘，如果得到叉乘得到的向量都朝着相同方向那么在三角形内，否则在三角形外。上图的话显然Q点就在三角形外。

　　经过上两步的判断，现在就可以求得光线和三角形的交点了。

Möller Trumbore Algorithm

上面的一般方法有点繁琐，这里介绍一种更快更好的算法Möller Trumbore Algorithm(简称MT算法)。

推导过程

首先对于三角形内任意一点P，根据重心坐标(barycentric coordinates)我们可以得到P和三角形三顶点的关系:$$P = wA + uB + vC$$ 由于w,u,v满足$$w + u + v = 1$$所以有: $$P = (1 - u - v)A + uB + vC$$展开可得:$$P=A - uA - vA + uB + vC = A + u(B - A) + v(C - A)$$注意这里(B - A)和(C - A)分别表示的是△ABC的边AB和AC所对应的向量。假设光线和三角形的交点是P，那么P可以被写作:$$P=O+tD$$把该等式带入得:$$O+tD = A + u(B - A) + v(C - A)$$ $$O-A = -tD+u(B-A)+v(C-A)$$由线性代数的知识我们可以知道上面的方程可以用矩阵的形式表示:

\begin{array}{l} \begin{bmatrix} -D & (B-A) & (C-A) \end{bmatrix} \begin{bmatrix} t\\u\\v \end{bmatrix} =O-A \end{array}

我们可以解上面的方程:

[\begin{array}{r} t \\ u \\ v \end{array}] = \frac{1}{[| \begin{array}{r} - D & E_{1} & E_{2} \end{array} |]} [\begin{matrix} | \begin{matrix} T & E_{1} & E_{2} \end{matrix} | \\ | \begin{matrix} - D & T & E_{2} \end{matrix} | \\ | \begin{matrix} - D & E_{1} & T \end{matrix} | \end{matrix}],

根据克拉默法则(Cramer's Rule),每个未知数等于用常数项替换对应系数行列式某一列的行列式比上系数行列式。注意这里我们记

T = O - A

,

E_{1} = B - A

,

E_{2} = C - A

。根据标量三重积(scalar triple product)，我们知道(下面ABC是三个3维向量):$$|A B C| = -(A \times C) \cdot B = -(C \times B) \cdot A.$$所以最终结果为:

[\begin{array}{r} t \\ u \\ v \end{array}] = \frac{1}{(D \times E_{2}) \cdot E_{1}} [\begin{matrix} (T \times E_{1}) \cdot E_{2} \\ (D \times E_{2}) \cdot T \\ (T \times E_{1}) \cdot D \end{matrix}] == \frac{1}{P \cdot E_{1}} [\begin{matrix} Q \cdot E_{2} \\ P \cdot T \\ Q \cdot D \end{matrix}],

我们只用把对应的值带入上述矩阵，就可以解出对应未知量的值了。
那么如何判断光线是和三角形相交的呢?
1)如果开启了背面剔除

对于上面结论中的(D×E1)·E2这一项对应到上图，D对应光线方向，E1对应edge1,E2对应edge2,光线方向和E1的叉乘为pvec，于是可以知道，当光线方向和三角形的法线(E1 × E2)想对(反)的时候(D×E1)·E2 > 0，因此如果不满足这个条件，如下图所示:

当光线方向和发现方向相同的时候，(D×E1)·E2 < 0,此时光线就没有打到三角形!因为此时没有照到三角形的正面而照到的是背面。

或者当(D×E1)·E2 = 0时，即光线擦着三角面过去时我们也认为没有打到三角形！
2)如果不满足重心坐标的条件
即当我们解出变量u和v的时候，根据重心坐标可知如果点在三角形内，那么u+v应该是大于0小于等于1的，如果不满足这个条件，那么点在三角形外。

伪代码

下面实现了MT算法的光线与三角形求交。

Triangle.hpp rayTriangleIntersect

bool rayTriangleIntersect( 
    const Vec3f &orig, const Vec3f &dir, 
    const Vec3f &v0, const Vec3f &v1, const Vec3f &v2, 
    float &t, float &u, float &v) 
{ 
#ifdef MOLLER_TRUMBORE 
    Vec3f v0v1 = v1 - v0; 
    Vec3f v0v2 = v2 - v0; 
    Vec3f pvec = dir.crossProduct(v0v2); 
    float det = v0v1.dotProduct(pvec); 
#ifdef CULLING 
    // if the determinant is negative the triangle is backfacing
    // if the determinant is close to 0, the ray misses the triangle
    if (det < kEpsilon) return false; 
#else 
    // ray and triangle are parallel if det is close to 0
    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 
    ... 
#endif 
}