S01E07:点与平面的矩阵解法

329 阅读4分钟

说明

在 3D 几何中,有时平面是由矩阵来定义的:在局部坐标系中,指定 xy 平面为需要的平面,而局部坐标系的位置和姿态由矩阵定义。对于这样定义的平面,我们可以直接用矩阵解法来求解距离问题,及最近点问题。

而平面如果是由 3 个点,或者点+法线形式定义的,那么也可以强行构造矩阵,来使用矩阵解法,只是在构造过程中会引入精度误差,因此受到精度影响并不是理想的通用解法。

几何

当我们已知平面所在局部坐标系 T 时,要求点在平面的投影点和距离,只要将点的坐标转换到平面所在坐标系中就可以了,转换后的相对关系如图: 此时要求点 B 的投影点,其实就是将点 B 的 z 坐标置为 0 就可以了(因为平面的 z 坐标是 0),而距离则是点 B 的 z 坐标。那么如何将点 B 转换到局部坐标系呢?很简单,将点 B 乘以 T 的逆矩阵就可以了。

如果平面是由点+法线形式定义的,又该如何构造矩阵呢?答案是利用叉乘。比如下图,平面由点 A 和法线向量 AA' 来定义,点 B 位置如图。要构造局部坐标系,将AA'做为 z 轴,只需要求出 x 轴和 y 轴方向即可:

这里我们先构造向量 AB,这样就行成了 AB,AA'两个向量组成的平面,垂直于这个平面的,就是 y 轴;需要注意的是叉乘的顺序问题,AA' 叉乘 AB 才是 y 轴正方向。表示从 AA'向量,旋转到 AB 按右手螺旋法则,得到向量 AC 即为 y 轴正方向。

有了 y 轴后,再利用 y 轴和 z 轴,叉乘得到 x 轴,然后对三个轴做归一化,就可以构造矩阵了。这里同样,根据右手螺旋法则,从 y 轴旋转到 z 轴,得到的正方向是 x 轴正方向。

将归一化后的 x 轴写在矩阵第一列,y 轴写在第二列,z 轴写在第三列,平移(原点 A的坐标)写在第四列,我们就得到这个平面的矩阵。然后与前面保持一致,利用逆矩阵求出点 B 在局部坐标系中的坐标,就可以得到距离和投影点了。

最后,投影点的坐标需要再转换到世界坐标中,所以再乘以平面矩阵即可。而点到平面的距离,点是否在平面上,并不随坐标系改变而改变,无需再做转换。

代码

// 返回值为:投影点坐标,距离,是否在平面上
static func matrixRelationship(point:simd_float3, plane:Plane) -> (simd_float3, Float, Bool) {
    let vector = point - plane.position
    let yAxisVector = cross(plane.normal, vector)
    
    if length_squared(yAxisVector) < 0.0001{
        // 点在平面上的投影点,距离平面原点太近,即 vector 与 plane.normal 几乎共线。
        let distance2 = distance(point, plane.position)
        
        return (plane.position, distance2, distance2 < 0.0001)
    }
    let xAxisVector = cross(yAxisVector, plane.normal)
    if length_squared(xAxisVector) < 0.0001 {
        // 点与平面原点距离过近
        return (plane.position, 0, true)
    }
    let xAxis = simd_float4(normalize(xAxisVector), 0)
    let yAxis = simd_float4(normalize(yAxisVector), 0)
    let zAxis = simd_float4(normalize(plane.normal), 0)
    let origin = simd_float4(plane.position, 1)
    
    let matrix = matrix_float4x4(
        xAxis,
        yAxis,
        zAxis,
        origin
    )
    
    let pointP = simd_float4(point, 1)
    let localP = matrix.inverse * pointP
    // 将 z 坐标置 0,得到投影点,再补 1,构成齐次坐标,方便与矩阵相乘
    let localProject = simd_float4(localP.x, localP.y, 0, 1)
    let projectP = matrix * localProject
    let projectP2 = simd_float3(projectP.x, projectP.y, projectP.z)
    
    let localDistance = localP.z
    
    return (projectP2, localDistance, localDistance < 0.0001)
}

其他

如果平面本身就是由矩阵定义的,那利用矩阵来求解是方便又快速的方法。同样道理,直线,也可以由矩阵定义,也可以利用矩阵来求解,与本文思路相同。矩阵解法通用性很强,而且思路简单。

惟一遗憾的是,如果不是由矩阵直接定义的,那在构造矩阵的过程中,多次的叉乘,多次归一化,以及求解矩阵的逆矩阵,都会引入一定的误差,让精度变差。因此,还要根据具体情况来使用。