浅析基于二维轮廓线重构表面算法

轮廓线重构算法

　　由一组二维轮廓线重建出物体的三维表面是三维数据场可视化中的一种表面绘制方法。在医学图像可视化以及其他可视化领域中有着广泛的应用。三维表面重建实际上是对物体表面进行三角形划分，从轮廓线的角度出发就是将轮廓线上的顶点按照一定规则进行三角形拼接，从而构成可视的三维物体表面，再利用三维显示技术将其显示出来。本文讨论了一种实现轮廓线重构的简易方法，其关键的步骤体现在相邻轮廓线的编织方法以及顶底层轮廓线的三角化上。并且用c++代码来实现之。



轮廓线数据结构表示

　　首先为了表示输入的轮廓线，需要使用一种数据结构来表示每一层的轮廓线。通常重建的输入会由N个轮廓线组成，每个轮廓线位于不同的高度位置，我们可以把轮廓线认为是以条在X-Y平面上闭合的曲线，而每条轮廓线有着不同的Z值。而二维曲线通常都是用一系列二维点组成的折线近似表示的，这样轮廓线的数据结构可以表示为如下的形式：

struct FloatDouble
{
float X;
float Y;
FloatDouble(float x,float y)
{
X=x;Y=y;
}
FloatDouble()
{
X=0;
Y=0;
}
};

struct Box3Float
{
public:
float Min3[3];
float Max3[3];
Box3Float()
{
Min3[0]=99999;
Min3[1]=99999;
Min3[2]=99999;
Max3[0]=-99999;
Max3[1]=-99999;
Max3[2]=-99999;
}
Box3Float(float minX, float minY, float minZ, float maxX, float maxY, float maxZ)
{
Min3[0] = minX;
Min3[1] = minY;
Min3[2] = minZ;
Max3[0] = maxX;
Max3[1] = maxY;
Max3[2] = maxZ;
}
void UpdateRange(float x, float y, float z)
{
if (x < Min3[0])
Min3[0] = x;
if (y < Min3[1])
Min3[1] = y;
if (z < Min3[2])
Min3[2] = z;
if (x > Max3[0])
Max3[0] = x;
if (y > Max3[1])
Max3[1] = y;
if (z > Max3[2])
Max3[2] = z;
}
float GetXLength()
{
return Max3[0]-Min3[0];
}
float GetYLength()
{
return Max3[1]-Min3[1];
}
float GetZLength()
{
return Max3[2]-Min3[2];
}
};

class ContourLine
{
private:
std::vector<FloatDouble> points;
float z;
public:
ContourLine():z(0) { }
void AddPoint2d(float x,float y)
{
FloatDouble t(x,y);
points.push_back(t);
}
inline std::vector<FloatDouble>& GetPointList()
{
return points;
}
inline int GetLinePointCount() const
{
return points.size();
}
inline void SetZ(float z)
{
this->z=z;
}
inline float GetZ() const
{
return z;
}
Box3Float GetBox()
{
Box3Float box;
for (size_t i = 0; i < points.size(); i++)
{
box.UpdateRange(points[i].X, points[i].Y, 0);
}
return box;
}
};

　　这样std::vector<ContourLine>即可以作为轮廓线重构算法的输入类型，重构算法最终输出代表结果表面的Mesh。Mesh结构在之前的文章中多次出现，就不详细描述了。在介绍轮廓线重构算法中，需要首先解决下面的几个子算法问题，所以为了正确理解算法逻辑，首先需要介绍一下这些子算法。

子算法一：判断折线圈是否为逆时针

　　上面说了一般用来表示平面上一条轮廓线的方法是使用一系列点。按这些点练成折线，最后再将最后一个点与第一个点连起来，便是一条闭合的折线，假设折线是不自交的，则这个线圈一定能够判断它是逆时针还是顺时针。在轮廓线重构算法中，需要保证每一层的轮廓线为是由逆时针的折现圈表示的。所以，这里介绍一种判断折线圈是否为逆时针的方法。

　　首先在计算几何中，利用解析几何的原理我们知道向量的叉积可以用来判断向量拐弯的方向，例如下图中向量A×B若为正数，则意味着B相对A是向逆时针方向旋转的；为0则A,B共线；否则B是相对A顺时针旋转的。而1/2*abs(A×B)则为AB向量构成三角形的面积。

AB×BC为正	AB×BC为负	逆时针的轮廓线	顺时针的轮廓线

　　因此一种直观的判断由点A1、A2、......An构成的折线圈是否为逆时针的方法是计算三角形(A1,A2,A3)，(A2,A3,A4)，(A3,A4,A5).....(An,A1,A2)的面积之和是否为正值。这里不做原理性证明为什么可以这样，详细的原理可以从如下网址查看：

　　原网址：http://www.mathopenref.com/coordpolygonarea2.html ，打不开的话可见CSDN的备份：http://blog.csdn.net/swfa1/article/details/18146581 。

　　这里只附上相应的代码：

static float Area(std::vector<FloatDouble>& contour)
{
int n = contour.size();
float A = 0.0f;
for (int p = n - 1, q = 0; q < n; p = q++)
{
A += contour[p].X * contour[q].Y - contour[q].X * contour[p].Y;
}
return A * 0.5f;
}

　　上述函数返回true，则折线圈为逆时针。

子算法二： 两层轮廓线编织算法

　　所谓轮廓线编织，就是把相邻两层的轮廓线上的点采用合适的方式连接起来，形成三角网，但这种连接的方式显然是不唯一的。下图显示了两种不同的编织方式形成的三角网，显然，两者存在形态上的差别。

一种好的编织方式	一种看上去不太好的编织方式

　　从直观感觉上看，显然前一个三角网中三角形更为规则（较少的钝角三角形）。所以究竟如何生成好的三角网，需要采用合适的策略。这样的策略在相关领域有很多研究成果，提出了不少的方法，这里介绍一个容易理解的法则，叫做最短对角线原则。下图为了说明这个原则，将轮廓线剪开并展平成直线，轮廓线的点则分布在这条直线上，相隔不等的距离。



　　将轮廓线编织模拟成在两条直线上连接相应的点的连接，采用最短对角线原则连接方式为：若当前Ai已经与Bj相连接，检查由四个点(Ai,Ai+1,Bj+1,Bj)组成的四边形，检测对角线AiBj+1与Ai+1Bj的长度，选择较短的一条作为新三角形的边。

最短对角线连接方式	非最短对角线的连接方式

　　当创建第一条连接线的时候，对于点A0，需要遍历一次B点的数组，寻找到离A0最近的B点，假设为Bk，则需要将B数组调整为以Bk为首元素，与A方向相同(顺逆时针)的新数组。所以这就解释了为何需要先清楚子算法一的内容。这样，在编织活动开始的时候，A0与B0是最近的点对，A0B0是连接的第一条线，之后的循环过程按最短对角线原则进行。在实现过程中需要使用两个指示变量分别标记A数组和B数组中行进到的位置，有点算法基础的同学会发现这与归并排序中的Merge过程有点相似。循环到A、B数组中有一个访问完为止，这时另一个数组余下的点可以与访问完的数组的最后一个点连接成若干三角形。最后还需要注意要将最后的An-1Bm-1与A0B0组成的四边形用最短对角线法再三角化一次。

　　在应用这个原则时，当两层轮廓线投影在平面上的位置偏差较大时可能出现的问题，例如下图的情况，在投影偏差较大的时候，可能会导致所有三角形的顶点都是同一点，如下图所示，这样显然没有达到编织轮廓线的目的，解决的方案是将两层轮廓线的几何中心移在一起，这样就会避免这种现象。



　　综上，编织由点数组A、B表示的轮廓线的过程用文字表述如下：

检查A、B数组的方向，将A，B数组调整为逆时针。
获取A、B的Box范围并保存，然后将A，B平移到使几何中心到原点(0,0)。
为A0寻找B中的最近点，然后以此最近点为B0，将B数组按顺序重新编号。
设置A、B数组的位置指示变量i，j，初始化为0
当i与j均未达到A和B的长度

检查(Ai,Ai+1,Bj+1,Bj)组成的四边形，若AiBj+1短于Ai+1Bj，则创建三角片(Ai,Bj,Bj+1)否则创建三角片(Ai,Bj,Ai+1)

若i达到A的长度，则将B中余下的点连A的最后一点构成三角片
若j达到B的长度，则将A中余下的点连B的最后一点构成三角片
为四边形(An-1,Bm-1,A0,B0)构建按最短对角线原则创建三角片。

　　下列图片序列简单的反映了这一过程：

初始四边形
第二次迭代
第三次迭代
	迭代直到有一个点数组访问完毕
迭代完成

　　相应的实现代码，写成一个类形式如下，其中QuadUnit类用来表示一个用来使用最短对角线原则的四边形：

class ContourStitcher
{
private:
struct QuadUnit
{
public:
int UpIndex0;
int UpIndex1;
int DownIndex0;
int DownIndex1;
double DiaU0D1Len;
double DiaU1D0Len;
std::vector<FloatDouble>* lineUp;
std::vector<FloatDouble>* lineDown;
void Init(int upIndex, int downIndex)
{
UpIndex0 = upIndex;
DownIndex0 = downIndex;
UpIndex1 = (upIndex + 1); //%lineUp.Count;
DownIndex1 = (downIndex + 1); //%lineDown.Count;
DiaU0D1Len = GetDLen(UpIndex0, DownIndex1);
DiaU1D0Len = GetDLen(UpIndex1, DownIndex0);
}
void InitLast()
{
UpIndex0 = lineUp->size() - 1;
DownIndex0 = lineDown->size() - 1;
UpIndex1 = 0; //%lineUp.Count;
DownIndex1 = 0; //%lineDown.Count;
DiaU0D1Len = GetDLen(UpIndex0, DownIndex1);
DiaU1D0Len = GetDLen(UpIndex1, DownIndex0);
}
double GetDLen(int index1, int index2)
{
float x0 = (*lineUp)[index1].X;
float y0 = (*lineUp)[index1].Y;
float x1 = (*lineDown)[index2].X;
float y1 = (*lineDown)[index2].Y;
return sqrt((x0 - x1) * (x0 - x1) + (y0 - y1) * (y0 - y1));
}
};
void Transform(std::vector<FloatDouble>& list, float dx, float dy)
{
for (size_t i = 0; i < list.size(); i++)
{
list[i].X += dx;
list[i].Y += dy;
}
}
Point3d GetCenter(Box3Float& box3)
{
return Point3d((box3.Max3[0] + box3.Min3[0]) / 2.0f, (box3.Max3[1] + box3.Min3[1]) / 2.0f, (box3.Max3[2] + box3.Min3[2]) / 2.0f);
}
public:
ContourLine* lineUp;
ContourLine* lineDown;
std::vector<FloatDouble> lineUpProcessed;
std::vector<FloatDouble> lineDownProcessed;
Box3Float boxUp;
Box3Float boxDown;
ContourStitcher(ContourLine* line1, ContourLine* line2)
{
if (line1->GetZ() > line2->GetZ())
{
this->lineUp = line1;
this->lineDown = line2;
}
else
{
this->lineUp = line2;
this->lineDown = line1;
}
lineUpProcessed.reserve(lineUp->GetLinePointCount());
lineDownProcessed.reserve(lineDown->GetLinePointCount());
CopyArray(lineUp->GetPointList(), lineUpProcessed);
CopyArray(lineDown->GetPointList(), lineDownProcessed);
boxUp = lineUp->GetBox();
boxDown = lineDown->GetBox();
Point3d cU = GetCenter(boxUp);
Point3d cD = GetCenter(boxDown);
Transform(lineDownProcessed, -cD.X, -cD.Y);
Transform(lineUpProcessed, -cU.X, -cU.Y);
int indexDown = GetNearIndex();
AdjustDownArray(indexDown);
}
~ContourStitcher()
{
lineUp=NULL;
lineDown=NULL;
}
private:
void CopyArray(std::vector<FloatDouble> &olist, std::vector<FloatDouble> &tarlist)
{
for (size_t i = 0; i < olist.size(); i++)
{
tarlist.push_back(FloatDouble(olist[i].X, olist[i].Y));
}
}
int ContourStitcher::GetNearIndex()
{
int index = -1;
double distense = DBL_MAX;
FloatDouble& p = lineUpProcessed[0];
float x0 = p.X;
float y0 = p.Y;
for (size_t i = 0; i < lineDownProcessed.size(); i++)
{
float x1 = lineDownProcessed[i].X;
float y1 = lineDownProcessed[i].Y;
double dis = sqrt((x0 - x1) * (x0 - x1) + (y0 - y1) * (y0 - y1));
if (dis < distense)
{
distense = dis;
index = i;
}
}
return index;
}
void AdjustDownArray(int indexDown)
{
std::vector<FloatDouble> list;
list.reserve(lineDownProcessed.size());
for (size_t i = 0; i < lineDownProcessed.size(); i++)
{
list.push_back(lineDownProcessed[(indexDown + i) % lineDownProcessed.size()]);
}
lineDownProcessed.swap(list);
}
public:
Mesh* DoStitching()
{
Mesh *m = new Mesh();
int* upMap=new int[lineUpProcessed.size()];
float dx1 = GetCenter(boxUp).X;
float dy1 = GetCenter(boxUp).Y;
float dx2 = GetCenter(boxDown).X;
float dy2 =GetCenter( boxDown).Y;
int* downMap=new int[lineDownProcessed.size()];
for (size_t i = 0; i < lineDownProcessed.size(); i++)
{
Point3d p(lineDownProcessed[i].X + dx2, lineDownProcessed[i].Y + dy2, lineDown->GetZ());
downMap[i] = m->AddVertex(p);
}
for (size_t i = 0; i < lineUpProcessed.size(); i++)
{
Point3d p(lineUpProcessed[i].X + dx1, lineUpProcessed[i].Y + dy1, lineUp->GetZ());
upMap[i] = m->AddVertex(p);
}
size_t upIndex = 0;
size_t downIndex = 0;
QuadUnit quad;
quad.lineDown = &lineDownProcessed;
quad.lineUp = &lineUpProcessed;
while (true)
{
if (upIndex == lineUpProcessed.size() - 1 || downIndex == lineDownProcessed.size() - 1)
{
break;
}
quad.Init(upIndex, downIndex);
if (quad.DiaU0D1Len < quad.DiaU1D0Len)
{
Triangle t(upMap[quad.UpIndex0], downMap[quad.DownIndex0], downMap[quad.DownIndex1]);
m->AddFace(t);
downIndex++;
}
else
{
Triangle t(upMap[quad.UpIndex0], downMap[quad.DownIndex0], upMap[quad.UpIndex1]);
m->AddFace(t);
upIndex++;
}
}
if (upIndex == lineUpProcessed.size() - 1 || downIndex == lineDownProcessed.size() - 1)
{
if (downIndex == lineDownProcessed.size() - 1)
{
int last = lineDownProcessed.size() - 1;
while (upIndex != lineUpProcessed.size() - 1)
{
Triangle t(downMap[last], upMap[upIndex + 1], upMap[upIndex]);
m->AddFace(t);
upIndex++;
}
}
else
{
int last = lineUpProcessed.size() - 1;
while (downIndex != lineDownProcessed.size() - 1)
{
Triangle t(upMap[last], downMap[downIndex], downMap[downIndex + 1]);
m->AddFace(t);
downIndex++;
}
}
}
quad.InitLast();
if (quad.DiaU0D1Len < quad.DiaU1D0Len)
{
Triangle t(upMap[quad.UpIndex0], downMap[quad.DownIndex0], downMap[quad.DownIndex1]);
Triangle t2(upMap[quad.UpIndex0], downMap[quad.DownIndex1], upMap[quad.UpIndex1]);
m->AddFace(t);
m->AddFace(t2);
}
else
{
Triangle t(upMap[quad.UpIndex0], downMap[quad.DownIndex0], upMap[quad.UpIndex1]);
Triangle t2(upMap[quad.UpIndex1], downMap[quad.DownIndex0], downMap[quad.DownIndex1]);
m->AddFace(t);
m->AddFace(t2);
}
delete[] upMap;
delete[] downMap;
return m;
}
};

　　下图是按此原则编织的两层轮廓线的例图：
子算法三：Mesh合并方法简述

　　在实现了相邻两层轮廓线的三角网格生成操作后，对于已经按Z高度排序的轮廓线C1到Cn，可以使用子算法二对<C1,C2>...<Cn-1,Cn>生成对应的三角网格，其编号为M1,M2,...Mn-1。

　　考虑到每一个Mesh都是由两层顶点组成，所以这些网格合并成一个大的M网格M需要考虑顶点去重的问题，一个简单的方式是为每层顶点建立一个映射数组f，数组存储着每层顶点在M中的索引位置，这样每当需要添加Mi的顶点和三角片进入M的时候，Mi中低一层的顶点由于肯定在Mi-1加入M时已经加入M，则只添加Mi中高一层的顶点即可，在顶点加入时更新其三角片的顶点新索引。这样在所有三角片加入M后，Mesh便生成完成。由于这个子算法属于算法细节，这里就不单独描述伪代码。包含这部分的代码会在总的实现中贴出来。

合并前独立的两个Mesh	合并后，一些顶点共用

子算法四：任意多边形三角化方法

　　在轮廓线层之间的三角网格建立起来后，往往还需要把顶面和底面的轮廓线三角化，使得整个三角网形成封闭的拓扑结构，这部分算法摘抄自John W. Ratcliff在http://www.flipcode.com/archives/Efficient_Polygon_Triangulation.shtml上贴出的代码，其输入为点集表示的折线，输出为三角片集合，三角片顶点索引指向点集，这份代码中的Area函数，即是子算法一中判断折线圈方向的函数，该函数计算一个有符号的面积，返回正则说明折线圈为逆时针方向。

　　下面简要说明一下该算法的和主要思想：先将点顺序调整为逆时针，然后沿着折线圈，迭代着创建连续三个顶点(Ai,Ai+1,Ai+2)组成的三角片，若该三角片的旋转方向也为顺时针，则将该三角片删去，直到最后剩下一个三角片为止。该算法的代码整合为一个静态类如下：

class PolyTriangulator
{
public:
static bool Process(std::vector<FloatDouble>& contour, std::vector<Triangle>& result)
{
int n = contour.size();
if (n < 3) return false;
int* V = new int
;
if (0.0f < Area(contour))
for (int v = 0; v < n; v++) V[v] = v;
else
for (int v = 0; v < n; v++) V[v] = (n - 1) - v;
int nv = n;
int count = 2 * nv;   /* error detection */
for (int m = 0, v = nv - 1; nv > 2; )
{
if (0 >= (count--))
{
return false;
}
int u = v; if (nv <= u) u = 0;     /* previous */
v = u + 1; if (nv <= v) v = 0;     /* new v    */
int w = v + 1; if (nv <= w) w = 0;     /* next     */
if (Snip(contour, u, v, w, nv, V))
{
int a, b, c, s, t;
a = V[u]; b = V[v]; c = V[w];
Triangle tri(a, b, c);
result.push_back(tri);
m++;
for (s = v, t = v + 1; t < nv; s++, t++) V[s] = V[t]; nv--;
count = 2 * nv;
}
}
delete []V;
return true;
}
static float Area(std::vector<FloatDouble>& contour)
{
int n = contour.size();
float A = 0.0f;
for (int p = n - 1, q = 0; q < n; p = q++)
{
A += contour[p].X * contour[q].Y - contour[q].X * contour[p].Y;
}
return A * 0.5f;
}
private:
static bool InsideTriangle(float Ax, float Ay, float Bx, float By, float Cx, float Cy, float Px, float Py)
{
float ax, ay, bx, by, cx, cy, apx, apy, bpx, bpy, cpx, cpy;
float cCROSSap, bCROSScp, aCROSSbp;

ax = Cx - Bx; ay = Cy - By;
bx = Ax - Cx; by = Ay - Cy;
cx = Bx - Ax; cy = By - Ay;
apx = Px - Ax; apy = Py - Ay;
bpx = Px - Bx; bpy = Py - By;
cpx = Px - Cx; cpy = Py - Cy;

aCROSSbp = ax * bpy - ay * bpx;
cCROSSap = cx * apy - cy * apx;
bCROSScp = bx * cpy - by * cpx;

return ((aCROSSbp >= 0.0f) && (bCROSScp >= 0.0f) && (cCROSSap >= 0.0f));
}
static bool Snip(std::vector<FloatDouble>& contour, int u, int v, int w, int n, int* V)
{
int p;
float Ax, Ay, Bx, By, Cx, Cy, Px, Py;
Ax = contour[V[u]].X;
Ay = contour[V[u]].Y;
Bx = contour[V[v]].X;
By = contour[V[v]].Y;
Cx = contour[V[w]].X;
Cy = contour[V[w]].Y;
if (0.00001f > (((Bx - Ax) * (Cy - Ay)) - ((By - Ay) * (Cx - Ax)))) return false;
for (p = 0; p < n; p++)
{
if ((p == u) || (p == v) || (p == w)) continue;
Px = contour[V[p]].X;
Py = contour[V[p]].Y;
if (InsideTriangle(Ax, Ay, Bx, By, Cx, Cy, Px, Py)) return false;
}
return true;
}
};

综合算法

　　在了解了这些子算法的过程之后，可以总结出轮廓线算法的总体步骤如下：

为输入轮廓线排序，确保按Z升序排列
调整输入轮廓线，确保为均逆时针
平移所有轮廓线，使其几何中心重合
对所有连续两相邻轮廓线，执行编织算法(子算法二)
合并子算法二生成的Mesh
使用子算法四三角化顶面和地面，并将三角片与上一步的Mesh合并
将轮廓线平移回原来位置。

　　这样一个集成子算法完成轮廓线重建的类代码如下：

class ContourLineSurfaceGenerator
{
private:
static bool CompareContourline(ContourLine* r1,ContourLine* r2)
{
return r1->GetZ() < r2->GetZ();
}
void Transform(std::vector<FloatDouble>& list, float dx, float dy)
{
for (size_t i = 0; i < list.size(); i++)
{
list[i].X += dx;
list[i].Y += dy;
}
}
Point3d GetCenter(Box3Float& box3)
{
return Point3d((box3.Max3[0] + box3.Min3[0]) / 2.0f, (box3.Max3[1] + box3.Min3[1]) / 2.0f, (box3.Max3[2] + box3.Min3[2]) / 2.0f);
}
private:
std::vector<ContourLine*> lines;
public:
ContourLineSurfaceGenerator(std::vector<ContourLine*> lines)
{
this->lines = lines;
}private:
void ReverseNormal(std::vector<Triangle> &list)
{
for (size_t i = 0; i < list.size(); i++)
{
int temp = list[i].P0Index;
list[i].P0Index = list[i].P2Index;
list[i].P2Index = temp;
}
}
bool IsNormalZ2(std::vector<FloatDouble> &vertices, std::vector<Triangle> &faces)
{
if (faces.empty())
{
return true;
}
else
{
Triangle &t = faces[0];
Point3d p0(vertices[t.P0Index].X, vertices[t.P0Index].Y, 0.0f);
Point3d p1(vertices[t.P1Index].X, vertices[t.P1Index].Y, 0.0f);
Point3d p2(vertices[t.P2Index].X, vertices[t.P2Index].Y, 0.0f);

Vector v = Triangle::CaculateNormal(p0,p1,p2);
if (v.Z > 0)
{
return true;
}
else
{
return false;
}
}
}
void FillDownHole(Mesh *m, std::vector<int>* maps, std::vector<FloatDouble> &downContourPrcessed, Box3Float& boxD)
{
Transform(downContourPrcessed, GetCenter(boxD).X, GetCenter(boxD).Y);
std::vector<Triangle> down;
PolyTriangulator::Process(downContourPrcessed, down);
if (IsNormalZ2(downContourPrcessed, down))
{
ReverseNormal(down);
}
for (size_t i = 0; i < down.size(); i++)
{
Triangle t = down[i];
t.P0Index = maps[0][t.P0Index];
t.P1Index = maps[0][t.P1Index];
t.P2Index = maps[0][t.P2Index];
m->AddFace(t);
}
}
void FillUpHole(Mesh *m, std::vector<int>* maps, std::vector<FloatDouble> &upContourPrcessed, Box3Float& boxU)
{
Transform(upContourPrcessed,GetCenter(boxU).X,GetCenter(boxU).Y);
std::vector<Triangle> up;
PolyTriangulator::Process(upContourPrcessed, up);
if (!IsNormalZ2(upContourPrcessed, up))
{
ReverseNormal(up);
}
for (size_t i = 0; i < up.size(); i++)
{
Triangle t = up[i];
t.P0Index = maps[lines.size() - 1][t.P0Index];
t.P1Index = maps[lines.size() - 1][t.P1Index];
t.P2Index = maps[lines.size() - 1][t.P2Index];
m->AddFace(t);
}
}
void ReverseArray(std::vector<FloatDouble> &list)
{
size_t count = list.size();
for (size_t i = 0; i < count / 2; ++i)
{
FloatDouble temp = list[count - i - 1];
list[count - i - 1] = list[i];
list[i] = temp;
}
}
public:
Mesh *GenerateSurface()
{
if (lines.size() <= 1)
{
return NULL;
}
std::sort(lines.begin(), lines.end(), CompareContourline);
for (size_t i = 0; i < lines.size(); i++)
{
std::vector<FloatDouble>& linepoints = lines[i]->GetPointList();
if (0.0f > PolyTriangulator::Area(linepoints))
{
ReverseArray(linepoints);
}
}
Mesh *m = new Mesh();
std::vector<int>* maps=new std::vector<int>[lines.size()];
for (size_t i = 0; i < lines.size(); i++)
{
maps[i].resize(lines[i]->GetLinePointCount(),-1);
}
std::vector<FloatDouble> upContourPrcessed;
Box3Float boxU = lines[lines.size() - 1]->GetBox();
std::vector<FloatDouble> downContourPrcessed;
Box3Float boxD = lines[0]->GetBox();
for (size_t i = 0; i < lines.size() - 1; i++)
{
ContourStitcher cs(lines[i],lines[i + 1]);
if (i == 0)
{
downContourPrcessed.insert(downContourPrcessed.end(),cs.lineDownProcessed.begin(),cs.lineDownProcessed.end());
}
if (i == lines.size() - 2)
{
upContourPrcessed.insert(upContourPrcessed.end(),cs.lineUpProcessed.begin(),cs.lineUpProcessed.end());
}
Mesh *submesh = cs.DoStitching();
size_t Z0Count = lines[i]->GetLinePointCount();
size_t Z2Count = lines[i + 1]->GetLinePointCount();
if (submesh->Vertices.size() != Z0Count + Z2Count)
throw std::exception();
for (size_t j = 0; j < Z0Count; j++)
{
if (maps[i][j] == -1)
{
maps[i][j] = m->AddVertex(submesh->Vertices[j]);
}
}
for (size_t j = 0; j < Z2Count; j++)
{
maps[i + 1][j] = m->AddVertex(submesh->Vertices[j + Z0Count]);
}
for (size_t j = 0; j < submesh->Faces.size(); j++)
{
Triangle t = submesh->Faces[j];
if (t.P0Index < (int)Z0Count)
{
t.P0Index = maps[i][t.P0Index];
}
else
{
t.P0Index = maps[i + 1][t.P0Index - Z0Count];
}

if (t.P1Index < (int)Z0Count)
{
t.P1Index = maps[i][t.P1Index];
}
else
{
t.P1Index = maps[i + 1][t.P1Index - Z0Count];
}

if (t.P2Index <(int)Z0Count)
{
t.P2Index = maps[i][t.P2Index];
}
else
{
t.P2Index = maps[i + 1][t.P2Index - Z0Count];
}
m->AddFace(t);
}
delete submesh;
}
FillUpHole(m,maps,upContourPrcessed, boxU);
FillDownHole(m, maps, downContourPrcessed, boxD);
delete[] maps;
return m;
}

};

实验效果

　　本实验采用的轮廓线一共有四层，每一层轮廓线由如下图四条在取自X-Y平面上格点的点组成。数据的预览图如下：

最下层轮廓	第2层轮廓	第3层轮廓	第4层轮廓

　　生成的结果预览如下：
　　本文的项目工程的下载：http://files.cnblogs.com/chnhideyoshi/BlogAlgCp.rar