Point triangle distance

介绍

该文章中的c++代码参考了ipc-toolkit的实现。

点到三角形的距离可以用如下定义：

\[ \begin{aligned} \text{distance}(\vec{\mathbf{x_p}}, \vec{\mathbf{x_t}_1}, \vec{\mathbf{x_t}_2}, \vec{\mathbf{x_t}_3}) &= \min_{\beta_1, \beta_2} \left\| \vec{\mathbf{x_p}} - ( \vec{\mathbf{x_t}_1} + \beta_1 (\vec{\mathbf{x_t}_2} - \vec{\mathbf{x_t}_1}) + \beta_2 (\vec{\mathbf{x_t}_3} - \vec{\mathbf{x_t}_1}) ) \right\| \\&s.t. \beta_1 \geq 0, \beta_2 \geq 0, \beta_1 + \beta_2 \leq 1 \end{aligned} \]

这是一个分断连续的函数，实际计算时可以根据点和三角形的位置关系分以下几种情况讨论，首先要将点投影到三角形所在的平面上：

1. 投影后，点在三角形内部，此时距离为点到三角形所在平面的距离。
2. 投影后，点在三角形的某个边朝外的半平面且投影在边上的点在边上，此时距离为点到边的距离。
3. 其他情况，此时距离为点到三角形的三个顶点的最小距离。

因此首先判断点到三角形所在平面的距离的类型。

点到平面的距离类型判断

一个可行的计算方式是把点\(\mathbf{\vec{x_p}}\)分别投影到三角形三条边各自张成的二维空间，得到6个投影后的参数，通过参数的数值判断。因为三角形三条边已知，可以先计算出三角形的法向量，对于任意一条边\(t_1 t_2\)，可以构建一个直角坐标系，三个坐标轴分别是边的方向向量\(\vec{e_0} = \mathbf{\vec{x_{t2}}} - \mathbf{\vec{x_{t1}}}\)，法向量和边的方向向量的叉乘\(e_1 = \vec{e_0} \times \vec{n}\)，法向量\(\vec{n}\)。那么因为\(\vec{e_0}, \vec{e_1}\)张成一个二维空间，可以将点\(\mathbf{\vec{x_p}}\)投影到这个平面上，得到坐标\((\lambda_0^{1,2}, \lambda_1^{1,2})\)。对于边\(t_2 t_3\)和\(t_3 t_1\)，同样可以得到坐标\((\lambda_0^{2,3}, \lambda_1^{2,3})\)和\((\lambda_0^{3,1}, \lambda_1^{3,1})\)。

最小二乘法计算投影坐标

要求解坐标，可以构建一个\(3\times 2\) 的矩阵\(A\)，其中两个列向量分别是是边的方向向量\(\vec{e_0}\)和法向量\(\vec{e_1}\)，计算投影坐标等价于找到该空间内的一个点使得它与\(\mathbf{\vec{x_p}}\)的距离最短。可以通过最小二乘法求解：

\[ \begin{aligned} \min_{\vec{x}} & \left\| A \vec{x} - \mathbf{\vec{x_p}} \right\| ^2 \\ &= \min_{\vec{x}} \vec x^T A^T A \vec x - 2 \mathbf{\vec{x_p}}^T A \vec x + \mathbf{\vec{x_p}}^T \mathbf{\vec{x_p}} \\ s.t. & A = \begin{bmatrix} e_{01} & e_{02} \\ \end{bmatrix}\\ & x = \begin{bmatrix} \lambda_0^{1,2} \\ \lambda_1^{1,2} \end{bmatrix} \end{aligned} \]

令梯度为0，可以得到：

\[ \begin{aligned} \frac{\partial}{\partial \vec{x}} \left( \vec x^T A^T A \vec x - 2 \mathbf{\vec{x_p}}^T A \vec x + \mathbf{\vec{x_p}}^T \mathbf{\vec{x_p}} \right) &= 0 \\ 2 A^T A \vec x - 2 A^T \mathbf{\vec{x_p}} &= 0 \\ A^T A \vec x &= A^T \mathbf{\vec{x_p}} \\ \vec x &= (A^T A)^{-1} A^T \mathbf{\vec{x_p}} \end{aligned} \]

c++代码如下：

template<typename Scalar>
__forceinline__ __device__ glm::tvec2<Scalar> computeProjectedCoordinate(
    const glm::tvec3<Scalar>& p,
    const glm::tvec3<Scalar>& t0,
    const glm::tvec3<Scalar>& t1,
    const glm::tvec3<Scalar>& normal) {
    glm::tmat2x3<Scalar> basis;
    basis[0] = t1 - t0;
    basis[1] = glm::cross(basis[0], normal);
    glm::tmat2x2<Scalar> basisT_basis = glm::tmat2x2<Scalar>(
        glm::dot(basis[0], basis[0]), glm::dot(basis[0], basis[1]),
        glm::dot(basis[1], basis[0]), glm::dot(basis[1], basis[1])
    );
    return glm::inverse(basisT_basis) * glm::tvec2<Scalar>(
        glm::dot(basis[0], p - t0), glm::dot(basis[1], p - t0)
    );
}

点到平面的距离类型的三种情况

1. 当\(\lambda_0^{1,2}\)满足\(0 \lt \lambda_0^{1,2} \lt 1\)，\(\lambda_1^{1,2}\)满足\(\lambda_1^{1,2} \geq 0\)时，距离为点到边\(t_1 t_2\)所在直线的距离。其他边以此类推。
2. 如果点对于三条边的坐标都不满足条件1，那么\(\lambda_0^{1,2} \lt 0, \lambda_0^{3,1} \ge 1\)，则距离为点到\(\mathbf{x_{t1}}\)的距离；\(\lambda_0^{2,3} \lt 0, \lambda_0^{1,2} \ge 1\)，则距离为点到\(\mathbf{x_{t2}}\)的距离；\(\lambda_0^{3,1} \lt 0, \lambda_0^{2,3} \ge 1\)，则距离为点到\(\mathbf{x_{t3}}\)的距离。
3. 如果点对于三条边的坐标都不满足条件1和2，那么距离为点到三个平面的最小距离。

c++代码如下：

template<typename Scalar>
__device__ DistanceType point_triangle_distance_type(
    const glm::tvec3<Scalar>& p,
    const glm::tvec3<Scalar>& t0,
    const glm::tvec3<Scalar>& t1,
    const glm::tvec3<Scalar>& t2)
{
    glm::tvec3<Scalar> normal = glm::cross(t1 - t0, t2 - t0);
    glm::tmat3x2<Scalar> param;

    param[0] = computeProjectedCoordinate(p, t0, t1, normal);
    if (param[0][0] > 0.0 && param[0][0] < 1.0 && param[0][1] >= 0.0) {
        return DistanceType::P_E0;
    }

    param[1] = computeProjectedCoordinate(p, t1, t2, normal);
    if (param[1][0] > 0.0 && param[1][0] < 1.0 && param[1][1] >= 0.0) {
        return DistanceType::P_E1;
    }

    param[2] = computeProjectedCoordinate(p, t2, t0, normal);
    if (param[2][0] > 0.0 && param[2][0] < 1.0 && param[2][1] >= 0.0) {
        return DistanceType::P_E2;
    }

    if (param[0][0] <= 0.0 && param[2][0] >= 1.0) {
        return DistanceType::P_T0;
    }
    else if (param[1][0] <= 0.0 && param[0][0] >= 1.0) {
        return DistanceType::P_T1;
    }
    else if (param[2][0] <= 0.0 && param[1][0] >= 1.0) {
        return DistanceType::P_T2;
    }

    return DistanceType::P_T;
}

计算距离

下面给出的c++代码计算的都是距离的平方。

点到点的距离

点到点的距离是两点之间的欧氏距离。

\[ \text{distance}(\vec{\mathbf{x_p}}, \vec{\mathbf{x_q}}) = \left\| \vec{\mathbf{x_p}} - \vec{\mathbf{x_q}} \right\| \]

c++代码如下：

template<typename Scalar>
__device__ Scalar point_point_distance(
    const glm::tvec3<Scalar>& p0,
    const glm::tvec3<Scalar>& p1)
{
    return glm::length2(p0 - p1);
}

点到直线的距离

点到直线的距离可以通过计算点与直线上两点连线得的三角形的面积的2倍除以直线上两点的长度得到。

\[ \text{distance}(\vec{\mathbf{x_p}}, \vec{\mathbf{x_t}_1}, \vec{\mathbf{x_t}_2}) = \frac{\left\| (\vec{\mathbf{x_p}} - \vec{\mathbf{x_t}_1}) \times (\vec{\mathbf{x_p}} - \vec{\mathbf{x_t}_2}) \right\|}{\left\| \vec{\mathbf{x_t}_2} - \vec{\mathbf{x_t}_1} \right\|} \]

c++代码如下：

template<typename Scalar>
__device__ Scalar point_line_distance(
    const glm::tvec3<Scalar>& p,
    const glm::tvec3<Scalar>& e0,
    const glm::tvec3<Scalar>& e1)
{
    return glm::length2(glm::cross(e0 - p, e1 - p)) / glm::length2(e1 - e0);
}

点到平面的距离

点到平面的距离是点到平面的法向量的投影。

\[ \text{distance}(\vec{\mathbf{x_p}}, \vec{\mathbf{x_t}_1}, \vec{\mathbf{x_t}_2}, \vec{\mathbf{x_t}_3}) = \frac{\left| (\vec{\mathbf{x_p}} - \vec{\mathbf{x_t}_1}) \cdot \vec{n} \right|}{\left\| \vec{n} \right\|} \]

c++代码如下：

template<typename Scalar>
__device__ Scalar point_plane_distance(
    const glm::tvec3<Scalar>& p,
    const glm::tvec3<Scalar>& origin,
    const glm::tvec3<Scalar>& normal)
{
    const Scalar point_to_plane = glm::dot(p - origin, normal);
    return point_to_plane * point_to_plane / glm::dot(normal, normal);
}