有限元基础知识：如何计算OBB包围盒
之前已经讲过AABB的编码实现，今天给大家讲一下OBB的编码实现，与AABB的编码不同，OBB会更紧贴物体，对于有限元接触仿真来说，会更贴近我们这一团网格，所以可以在前期搜索的时候更好的进行筛选。然而，OBB的计算会相对繁琐很多请看下边
PCA(Principal Component Analysis)算法找到一堆节点的主轴
首先找到一堆节点的的中心点
其中 为网格节点的坐标，具体C++代码如下
using point3D = std::array<double, 3>;
void findCentroid(const std::vector<point3D> & pts, point3D & centroid)
{
int number = pts.size();
centroid[0] = centroid[1] = centroid[2] = 0.0;
for (int i = 0; i< number; i++)
{
centroid[0] += pts[i][0];
centroid[1] += pts[i][1];
centroid[2] += pts[i][2];
}
centroid[0] /= n;
centroid[1] /= n;
centroid[2] /= n;
}
在得到了中心点坐标后，我们计算这一群节点与中心点的协方差矩阵，这个东西，大家顾名思义，可以理解成，这一群节点相较之这个节点有多分散，具体的算法如下
这个代码也很好写，cpp代码如下
void computeCovariance(const std::vector<point3D> & pts)
{
point3D centroid;
findCentroid(pts, centroid);
std::array<point3D, 3> covMat {0};
for (const auto& pt : pts) {
double deltax = pt[0] - centroid[0];
double deltay = pt[1] - centroid[1];
double deltaz = pt[2] - centroid[2];
covMat[0][0] += deltax * deltax;
covMat[1][1] += deltay * deltay;
covMat[2][2] += deltaz * deltaz;
covMat[0][1] += deltax * deltay;
covMat[0][2] += deltax * deltaz;
covMat[1][2] += deltay * deltaz;
}
for (int i = 0; i < 3; ++i) {
for (int j = 0; j < 3; ++j) {
covMat[i][j] /= pts.size();
}
}
covMat[1][0] = covMat[0][1];
covMat[2][0] = covMat[0][2];
covMat[2][1] = covMat[1][2];
}
然后我们通过提取协方差矩阵的特征向量，找到主轴，这个提取特征值的过程我们可以直接调用矩阵求解库中一般都有的特征值求解函数，或者由于这就是个 矩阵的特征向量提取问题，我们可也可以自己写个更简单高效的[1]，这里代码略，大家在很多地方都可以查到，大体上有3种写法，分列如下。
void eigenMat33(const Mat33 & matrix, Mat33& eigenvectors, point3D& eigenvalues) {
// 1. 采用数学库比如 Eigen直接进行计算
// 2. 采用迭代求解特征值、特征向量
// 3. 3x3 矩阵特征值有解析求法，实现一下就行了
}
然后我们就得到了特征值与特征向量，这个特征向量3个数据方差最大的方向，也就是能找到的OBB的两两相对面的方向，现在我们终于可以构造OBB了，非常简单，其实就是将中点到各点的连线往这3个轴上进行投影，而投影嘛就是点乘，大家就记住这个操作就行了，然后简简单单的找个各个方向的最大值最小值，OBB就呼之欲出了，C++代码如下
struct OBB
{
Mat33 axes;
point3D center;
point3D halfExtent; // 半轴长
};
OBB computeOBB(const std::vector<point3D>& pts, const Mat33 & eigenvectors) {
OBB obb;
obb.axes[0] = eigenvectors[0]; // 最大特征值对应方向
obb.axes[1] = eigenvectors[1];
obb.axes[2] = eigenvectors[2];
// 计算各轴投影范围
Vec3 maxProj{-FLT_MAX, -FLT_MAX, -FLT_MAX};
Vec3 minProj{FLT_MAX, FLT_MAX, FLT_MAX};
for (const auto& pt : pts) {
point3D dir = {pt[0], pt[1], pt[2]};
for (int i = 0; i < 3; ++i) {
double proj = dot(dir, obb.axes[i]);
maxProj[i] = std::max(maxProj[i], proj);
minProj[i] = std::min(minProj[i], proj);
}
}
// 计算中心点和半长
obb.halfExtent = {(maxProj[0]-minProj[0])/ 2.0,
(maxProj[1]-minProj[1])/ 2.0, (maxProj[2]-minProj[2])/ 2.0};
for (int i = 0; i< 3; i++)
{
obb.center += 0.5 * (maxProj[0] + maxProj[0])*obb.axes[i];
}
return obb;
}
以上的代码就通过3个方向的投影计算了各个方向的半轴长度，与OBB的中心，注意哦，这里OBB的中心并不是刚刚计算出来的那个 , 一个简单的案例就是下边这个三角形：
就这样我们就得到了一片Mesh的OBB，然而这种方法其实在有限元中用的不多，大家可以看到这个计算OBB的过程已经很繁琐了，然而其在游戏、CAD中用的还是不少的，但本着大家学就有用的思想，技多不压身，后面我可以单独讲一讲为啥他们在这些行业会有用。
另外就是大家可能会在别的地方看到不同的OBB计算方法，但是基本思想是一致的，但是OBB的计算方法如果展开说我可以直接开一系列，最常用的一个变形就是先计算一堆节点的凸包（convex hull）再用其来计算OBB，但既然我还没讲过convex hull先就此略过，以后再补。
然后还有更麻烦的呢，OBB的接触检测又需要用到一个算法，分离轴定理（SAT），其实也可以直接用GJK算法，这些也挺繁琐，有待下回分解。
参考文献：
[1] Siddique, Abu Bakar, and Tariq A. Khraishi. “Eigenvalues and Eigenvectors for 3x3 Symmetric Matrices: An Analytical Approach.” Journal of Advances in Mathematics and Computer Science 35.7 (2020): 106-118
