本篇介绍一种近似平面的算法和实现，要解决的问题是：给定点集，确定尽可能经过这些点集的近似平面。
文章目录：
- 算法介绍
- 源码实现
- 参考
算法介绍
在三维空间上，有一个多边形
，确定经过该多边形的平面，并不能保证多边形上的顶点精确的在同一个平面上。
计算经过多边形的平面，有一种非常简单的方法，就是从中任选三个顶点，使用叉积，计算出法向量，再计算出平面表示中的常量。这种方法存在两个问题：（1）如果所有点不能精确的在同一个平面上，任选的三个点不具代表性，无法代表点集中所有的点；（2）如果选择的三个点，“几乎”（很接近，但是又不是）是同一条直线上，计算叉积的结果会非常的小，则误差就会变得非常的大。
另外一种方法，是从顶点集中任选三个点匹配，那么就会有
个组合，计算出它们的法向量，排除掉长度很小的法向量，归一化所有的法向量，最终取所有向量的平均值。这是一种估计法向量的方法，但是算法的复杂度达到
，效率成了这种方法的瓶颈。
Martin Newell提出了一种计算经过多边形
的平面的估计方法[1]，考虑到所有点对法向量的贡献。设求得的法向量是
，它的计算公式如等式（1）所示：
其中，
表示顶点的个数，
表示第
个点，
，计算出的法向量方向遵守右手法则，如果沿着法向量向下方向观察，顶点是以逆时针方向排序的，如图1所示。
举个例子，给定多边形的三个点
，使用叉积的方法可以求得法向量
现在使用等式（1），同样可以计算出最终的法向量是
。
在求得法向量
后，还需要求得平面经过的一点V，可以通过计算所有点的平均值，采用来求点V，如下所示：
知道了平面的法向量和平面上的点V，就可以很容易计算出平面的表示公式。接下来，介绍该公式的推导过程，解释为什么它能表示成经过经过多边形的“近似”平面。
把一个多边形正交投影到
以及
的平面上，可以得到3个2D多边形，3个2D多边形的面积分别为Ax，Ay，Az。首先需要证明的一个子结论是:法向量
与(Ax，Ay，Az)成比例，即
=k(Ax，Ay，Az)。
1. 假设多边形P是一个三角形（多边形可以表示成多个三角形的和），如图2所示。三角形T所在的平面的法向量是
，是单位长度，三角形面积为Area(T)；三角形
是三角T在另一个平面上的投影，投影平面的法向量是
，也是单位长度，三角形面积为Area(
)。
（a）
，则有
；同理，计算投影后的三角形
，有
，则有
（b）考虑单个向量的投影，如图3所示，向量
投影到另一个平面，得到新的向量
，易知
，同理，可以计算出
；
（c）把
和
代入（a）中的等式，可以得到
由于
，所以可以消除最后一项，则有
（d）把上述等式两边点积上
，可以得到
2. 现在考虑一般多边形的情况，多边形P可以表示成多个三角形的和，即有
3. 多边形在
或者
的平面上的投影，相当于在法向量分别是
，
，
的平面上的投影，则(Ax, Ay, Az) = K(nx, ny, nz)，K是一个常数。法向量可以通过求多边形在多个平面上的投影面积，再将它归一化，就是多边形的法向量。
接下来，考虑计算多边形P在3个平面上的投影面积，以xy平面为例，即
。多边形P上的顶点
，投影到平面上，得到
，由顶点
与在轴对齐直线
构成一个梯形，计算梯形的面积，得到
每条边都可以求得一个面积，一般的等式是
有些梯形的面积可能是正，有些可能是负，所有面积的和，就可以得到图4所示的多边形的面积。即使存在两条边几乎在同一条直线上的情况，但是它们都为最后求得的法向量作出了贡献。
每条边的面积和为
同理，可以计算出多边形P在xz平面上的投影面积
和在yz平面上的投影面积
。
源码实现
基础库源码链接，参见这里，下面是前面所描述的算法的实现。
#include "SrGeometricTools.h"
#include "SrDataType.h"
namespace {
/**
brief 3D plane class.
This is a 3D plane class with public data members.
The line is parameterized as ^n*^X+d=0,in which ^n is the 'normal' data,d is the 'd' data.
The normal isn't normalized.
*/
class SrPlane3D
{
public:
/**
brief Initialize a plane by a set of points.
*/
bool initPlane(const SrPoint3D* point, int numPoint)
{
SrVector3D normal(0, 0, 0);
if (!computeNormal(point, numPoint, normal))
return false;
SrPoint3D avePoint = SrPoint3D(0, 0, 0);
int i;
for (i = 0; i < numPoint; i++)
{
avePoint += point[i];
}
avePoint /= numPoint;
mNormal = normal;
mD = -mNormal.dot(avePoint);
return true;
}
protected:
bool computeNormal(const SrPoint3D* point, int number, SrVector3D& result) const
{
SrVector3D normal(0, 0, 0);
int i;
for (i = 0; i < number; i++)
{
normal.x += (point[i].y - point[(i + 1) % number].y)*(point[i].z + point[(i + 1) % number].z);
normal.y += (point[i].z - point[(i + 1) % number].z)*(point[i].x + point[(i + 1) % number].x);
normal.z += (point[i].x - point[(i + 1) % number].x)*(point[i].y + point[(i + 1) % number].y);
}
if (EQUAL(normal.x, 0) && EQUAL(normal.y, 0) && EQUAL(normal.z, 0))
return false;
result = normal;
return true;
}
public:
SrVector3D mNormal;
SrReal mD;
};
}
参考
[1] Hill, F., and S. Kelley. Computer Graphics Using OpenGL, 3/E, Pearson, 2007.
[2] Philip Schneider, and David H. Eberly. Geometric tools for computer graphics, Morgan Kaufmann, 2002.