四叉树剖分平面区域
四叉树是剖分任意平面区域的重要方法，其核心思想是通过递归算法，离散区域为四边形，四边形边长满足单元尺寸要求。其算法稳定，思路清晰。因为单元形状被四边形限定，生成的网格质量较好。
1. 四叉树
四叉树的结构如下图，每个父节点只有四个子节点。采用链表实现，如果某个父节点没有子节点，四个指针为空；子节点同时存储其父信息，如果父信息为空指针，那么其是根节点。四叉树的数据结构大致如下：
为了满足划分网格，四叉树结构需满足细分四边形、查找相邻四边形的功能。部分代码如下，
struct Vertex //节点坐标{Vertex() :x(0.), y(0.) {}Vertex(double x_,double y_):x(x_),y(y_){}Vertex(const Vertex& other) :x(other.x), y(other.y) {}Vertex getMid(const Vertex& other) {return Vertex((x+other.x)/2,(y+other.y)/2);}double distance(const Vertex& other) {return std::pow(pow(x-other.x,2)+pow(y-other.y,2),0.5);}double x;double y;};//节点struct QUAD {QUAD() {}QUAD(const Vertex& v0, const Vertex& v1, const Vertex& v2 , const Vertex& v3){nodes[0] = v0;nodes[1] = v1;nodes[2] = v2;nodes[3] = v3;sons[0] = nullptr;sons[1] = nullptr;sons[2] = nullptr;sons[3] = nullptr;}void operator=(const QUAD& quad);Vertex nodes[4]; //存储四个节点QUAD* parent; //存储父节点，用于查找相邻单元QUAD* sons[4]; //存储四个子节点，分别是SE，NE，NW，SW；int level;//记录四边形的level。初始level为0，细化依次level++//QUAD* SE; //South East//QUAD* NE; //North East//QUAD* NW; //North West//QUAD* SW; //South West};class QuadsTree{public:QuadsTree();void addRoot(QUAD* quad) {Root = quad;Root->parent = nullptr;}bool splitQUAD(QUAD* current); //添加现有节点的四个子四边形void refine(double len); //如果四边形对角线长度大于len，细分void refine();//查找四边形的邻居。void getNorth(QUAD* current,QUAD* neighbor1, QUAD* neighbor2);//void getSouth(QUAD* current);//void getWest(QUAD* current);//void getEast(QUAD* current);void print(QUAD* current, std::vector<TopoDS_Shape>* shapes); //Root-> TL,TR,BL,BRvoid getQuads(std::vector<TopoDS_Shape>* shapes);public:QUAD* Root;void PreOrder(QUAD* current);void add(QUAD* current, int pos, QUAD* quad); //1-4----> BR,TR,TL,BLvoid refinement(QUAD* current, double len);void refinement(QUAD* current);double spaceFunction(QUAD* current);};
细分四边形的代码较为简单，当QUAD的对角线长度大于spaceFunction函数返回的尺寸时，将四边形二分。refinement函数采用递归算法，即可从Root开始，细化所有不满足的网格。
void QuadsTree::refinement(QUAD* current) {if (current->sons[0] == nullptr) {//计算对角线长度，细化。要满足“level restriction”，此处还应该考虑临近四边形的levelVertex v0 = current->nodes[0];Vertex v3 = current->nodes[3];double lenV0V3 = v0.distance(v3);double len = spaceFunction(current);if (lenV0V3 > len) {splitQUAD(current);}else {return;}}refinement(current->sons[0]);refinement(current->sons[1]);refinement(current->sons[2]);refinement(current->sons[3]);}
如果四边形指针为空，直接返回。 如果当前四边形的父是根节点。如果Current是NE，NW，北边没有相邻四边形，返回。 如果当前四边形是SW，SE，其北边为同级的NW，NE或者其子节点。此时只需要访问到Current的父一级就OK。
if ((current == parent->sons[0] ) || (current == parent->sons[3])) { //如果Current是 SW，SE；其邻居为同一级的NW，NE或其子节点QUAD* temp=nullptr;if ((current == parent->sons[0])) temp = parent->sons[1];if ((current == parent->sons[3])) temp = parent->sons[2];if (!temp->sons[0]) { //没有子节点neighbor1 = temp;neighbor2 = nullptr;return;}else {neighbor1 = temp->sons[0];neighbor2 = temp->sons[3];return;}}
2. 网格划分
使用四叉树进行网格划分的一般步骤如下，
生成平面几何图形的《Bounding Box》
将根据合适的Space function将BoundingBox细分为合适的尺寸。
移除边界外和靠近边界的四边形。
边界恢复。
连接四边形节点，形成连续的区域划分。
以下对这几方面内容将详细介绍，
生成二维图形的Bounding Box较为简单，参考《Bounding Box》，生成AABB形式的。
将区域细分为合适的尺寸，当四边形对角线长度大于中心点的SpaceFunction值时，即细分。并且，当我level值大于相邻四边形的level值时，必须先将相邻四边形细分，再细分此四边形。
当四边形的某一部分在区域外部时，删除。判断四个角点是否都在四边形内部，否则，删除。判断点在区域的内部、外部，参考《带孔平面的网格划分策略》所述的交点个数办法。删除后如下图，
2.4 边界恢复
如上图，删除四边形之后，原始图形边界与四边形边界有未填充的区域。边界恢复的目的是补全这部分区域，采用波前法（Advanced Front Technology），参考《AFT（波前法） 三角网格剖分》。此步的难点在于，查找四边形的边界，并且使其有向，形成闭合的“AFT前沿”。由于此时的原始边界为外边界（逆时针），四边形边界均为内边界（顺时针）。关于查找四边形边界，较为复杂，目前我并没有很清晰的思路。有以下两种思考：
使用四边形邻居关系查找：以某一段边界出发（某个没有邻居的边），沿着其邻居查找。如果某个四边形两条边都没有邻居，都是边界；如果某个四边形的相邻四边形中都是不是边界，在其邻居的邻居中查找，直到找到与前一段边界共节点的边界。如果边界节点与起始点重合，一段完整的边界结束。继续寻找下一段边界，直到所有四边形都找过。对闭合的四边形进行绕向整理，使其顺时针排序。
先找到所有边界，再找环，再绕向。所有没有邻居的边即为边界，以一条边为起始，依次找其共节点下一条边，直到找到闭合的环；以同样的方法，找到下一条环；直到所有边界都找完为止。
此两种方法的复杂度都不小，我并没有编码实现，以后有机会，会做些深入学习。找到四边形边界环如下左图，通过AFT划分之后的模型如下右图。
此六种类别的识别，通过其每条边相邻四边形的情况定义。
第一种：每条边没有邻居或只有一个邻居
某条边有2个邻居，并且其他边没有邻居或只有一个邻居
两条边有两个邻居，并且这两条边相邻。其余两条边没有邻居或只有一个邻居
两条边有两个邻居，并且这两条边相对。其余两条边没有邻居或只有一个邻居。
三条边有2个邻居，另一条边没有边界或只有一条边界。
四条边各有两个边界
参考：
基于自适应笛卡尔网格的虚拟单元方法研究
