最近畢設專案中用到了最大包絡體求算演算法，在這裡進行簡單的整理，為了以後更好的理解。

準備知識

• 關於點的定義

//空間上任何一個點資訊
struct Point {
	double x, y, z;
	Point(){} 
	Point(double xx,double yy,double zz):x(xx),y(yy),z(zz){}
	
	//兩向量之差  
    Point operator -(const Point p1)  
    {  
        return Point(x-p1.x,y-p1.y,z-p1.z);  
    } 

	//兩向量之和  
 

    Point operator +(const Point p1)  
    {  
        return Point(x+p1.x,y+p1.y,z+p1.z);  
    }  

    //叉乘  
    Point operator *(const Point p)  
    {  
        return Point(y*p.z-z*p.y,z*p.x-x*p.z,x*p.y-y*p.x);  
    } 
     
	// 數乘
    Point operator *(double d)  
    {  
        return Point(x*
 
d,y*d,z*d);  
    }  

	// 數除
    Point operator / (double d)  
    {  
        return Point(x/d,y/d,z/d);  
    }  

    //點乘  
    double  operator ^(Point p)  
    {  
        return (x*p.x+y*p.y+z*p.z);  
    } 
};

平面凸多邊形求算演算法

給一系列處於同一平面的空間點，然後求出所有隻在最外凸多邊形上的所有點集；其實為實現該目標有多種具體的演算法， 筆者將通過程式碼具體實現的方式將其中一種具體實現。

• 演算法步驟
  
• 程式碼實現

/*求平面內的最大包絡多邊形
引數解釋：平面內所有點資訊，用於儲存多邊形上下兩半的二維陣列，平面的法向量
*/
void dealWith(vector<Point> &allPoints, vector<Point> polygon[2], Point n1) {
	if(allPoints.size() < 2) return;
	Point a, b; //最小和最大兩個極端頂點；
	a.x = allPoints[0].x;
	a.y = allPoints[0].y;
	a.z = allPoints[0].z;
	b.x = allPoints[0].x;
	b.y = allPoints[0].y;
	b.z = allPoints[0].z;
	for(int i=1; i<allPoints.size(); i++) {
		if(a.x - allPoints[i].x > eps) {
			a.x = allPoints[i].x;
			a.y = allPoints[i].y;
			a.z = allPoints[i].z;
		} else if(fabs(a.x - allPoints[i].x) < eps) {
			if(a.y - allPoints[i].y > eps) {
				a.x = allPoints[i].x;
				a.y = allPoints[i].y;
				a.z = allPoints[i].z;
			} else if(fabs(a.y - allPoints[i].y) < eps) {
				if(a.z - allPoints[i].z > eps) {
					a.x = allPoints[i].x;
					a.y = allPoints[i].y;
					a.z = allPoints[i].z;
				} 
			}
		}

		if(allPoints[i].x - b.x > eps) {
			b.x = allPoints[i].x;
			b.y = allPoints[i].y;
			b.z = allPoints[i].z;
		} else if(fabs(b.x - allPoints[i].x) < eps) {
			if(allPoints[i].y - b.y > eps) {
				b.x = allPoints[i].x;
				b.y = allPoints[i].y;
				b.z = allPoints[i].z;
			} else if(fabs(b.y - allPoints[i].y) < eps) {
				if(allPoints[i].z - b.z > eps) {
					b.x = allPoints[i].x;
					b.y = allPoints[i].y;
					b.z = allPoints[i].z;
				} 
			}
		}
	}

	if (fabs(a.x - b.x) + fabs(a.y - b.y) + fabs(a.z - b.z) < eps) {
		polygon[0].push_back(a);
		printf("兩極值點相距過近，返回了直接");
		return;
	}

	polygon[0].push_back(a);
	polygon[0].push_back(b);
	polygon[1].push_back(a);
	polygon[1].push_back(b);
	vector<Point> p1, p2; // p1是直線左邊所有點集合，p2是直線右邊所有點集合
	Point mid ((a.x+b.x)/2, (a.y+b.y)/2, (a.z+b.z)/2); // 線段中點
	Point n2 (b.x-a.x, b.y-a.y, b.z-a.z); //兩個極值點的線段所在的向量
	Point n3 (n1.y*n2.z-n2.y*n1.z, n2.x*n1.z-n1.x*n2.z, n1.x*n2.y-n2.x*n1.y); // 計算所在平面內的線段的法向量
	for (int i = 0; i < allPoints.size(); ++i)
	{
		Point temp (allPoints[i].x-mid.x, allPoints[i].y-mid.y, allPoints[i].z-mid.z); //點集合中任意一個點到直線中點的向量
		double value = n3.x*temp.x + n3.y*temp.y + n3.z*temp.z; //向量和平面內直線法向量的點積
		if(value > eps) p1.push_back(allPoints[i]);
		else if(value < -eps) p2.push_back(allPoints[i]);
	}
	FindPoint2(p1, a, b, mid, polygon[0], n1);
	FindPoint2(p2, a, b, mid, polygon[1], n1);
}

//平面求包主題演算法
void FindPoint2(vector<Point> &p, Point a, Point b, Point mid, vector<Point> &polygon, Point &n) {
	if (p.size() == 0)
		return;
	Point pmax;
	pmax.x = p[0].x;
	pmax.y = p[0].y;
	pmax.z = p[0].z;
	double k, d;
	k = (b.y - a.y) / (b.x - a.x);
	d = a.y - k * a.x;
	double maxDis = DistanceOfPointToLine(&a, &b, &pmax), maxMid = distanceOfTwoPoints(pmax, mid);
	double newdist;
	for (int i = 1; i < p.size(); ++i)
	{
		newdist = DistanceOfPointToLine(&a, &b, &p[i]);
		if (newdist - maxDis > eps)
		{
			pmax.x = p[i].x;
			pmax.y = p[i].y;
			pmax.z = p[i].z;
			maxDis = newdist;
		}
		else if (fabs(newdist - maxDis) < eps)
		{	//選擇距離線段ab中點最近的那個
			double dis1 = distanceOfTwoPoints(p[i], mid);
			if (dis1 < maxMid)
			{
				pmax.x = p[i].x;
				pmax.y = p[i].y;
				pmax.z = p[i].z;
				maxMid = dis1;
			}
		}
	}
	polygon.push_back(pmax);

	Point mid1 ((pmax.x+a.x)/2, (pmax.y+a.y)/2, (pmax.z+a.z)/2);
	Point mid2 ((pmax.x+b.x)/2, (pmax.y+b.y)/2, (pmax.z+b.z)/2);
	Point v1 (mid1.x-mid.x, mid1.y-mid.y, mid1.z-mid.z);
	Point v2 (mid2.x-mid.x, mid2.y-mid.y, mid2.z-mid.z);
	Point l1 (pmax.x-a.x, pmax.y-a.y, pmax.z-a.z); //兩個極值點的線段所在的向量
	Point n1 (n.y*l1.z-l1.y*n.z, l1.x*n.z-n.x*l1.z, n.x*l1.y-l1.x*n.y); // 計算所在平面內的線段的法向量
	Point l2 (pmax.x-b.x, pmax.y-b.y, pmax.z-b.z); //兩個極值點的線段所在的向量
	Point n2 (n.y*l2.z-l2.y*n.z, l2.x*n.z-n.x*l2.z, n.x*l2.y-l2.x*n.y); // 計算所在平面內的線段的法向量
	if(v1.x*n1.x+v1.y*n1.y+v1.z*n1.z < -eps) {
		n1.x *= -1;
		n1.y *= -1;
		n1.z *= -1;
	}
	double len = sqrt(n1.x*n1.x+n1.y*n1.y+n1.z*n1.z);
	n1.x /= len;
	n1.y /= len;
	n1.z /= len;
	if(v2.x*n2.x+v2.y*n2.y+v2.z*n2.z < -eps) {
		n2.x *= -1;
		n2.y *= -1;
		n2.z *= -1;
	}
	len = sqrt(n2.x*n2.x+n2.y*n2.y+n2.z*n2.z);
	n2.x /= len;
	n2.y /= len;
	n2.z /= len;

	/* 找出各自符合滿足 Pmax,Pa 和 Pmax,Pb 的點 */
	vector<Point> p1, p2;
	for (int i = 0; i < p.size(); ++i)
	{
		Point temp1 (p[i].x-mid1.x, p[i].y-mid1.y, p[i].z-mid1.z);
		double value = temp1.x*n1.x+temp1.y*n1.y+temp1.z*n1.z;
		if(value > eps) p1.push_back(p[i]);
		else {
			Point temp2 (p[i].x-mid2.x, p[i].y-mid2.y, p[i].z-mid2.z);
			value = temp2.x*n2.x+temp2.y*n2.y+temp2.z*n2.z;
			if(value > eps) p2.push_back(p[i]);
		}

	}
	/* 遞迴尋找Pmax */
	FindPoint2(p1, pmax, a, mid1, polygon, n);
	FindPoint2(p2, pmax, b, mid2, polygon, n);
}

空間求凸包絡體演算法

空間凸包演算法，是給定一系列三維空間點，然後求出最小凸包絡體，其中凸包絡體的頂點都來自給定的點，並且任意點都在凸包中。

• 演算法步驟
  
• 實現程式碼

struct CH3D  
{  
    struct face  
    {  
        //表示凸包一個面上的三個點的編號  
        int a,b,c;  
        //表示該面是否屬於最終凸包上的面  
        bool ok;  
    };  

    //初始頂點數  
    int n;  

    //初始頂點  
    Point P[MAXN];  

    //凸包表面的三角形數  
    int num;  

    //凸包表面的三角形  
    face F[8*MAXN];  

    //凸包表面的三角形  
	//g[i][j]儲存的是第i個點連線到第j個點的有向向量所在的在F陣列中的三角面的序號
    int g[MAXN][MAXN]; 

	//共麵點集合，一維是集合數，二維是共面的點數
	vector<set<int>> count;
	vector<Point> polygons[MAXN][2];

    //向量長度  
    double vlen(Point a)  
    {  
        return sqrt(a.x*a.x+a.y*a.y+a.z*a.z);  
    }  

    //叉乘  
    Point cross(const Point &a,const Point &b,const Point &c)  
    {  
        return Point((b.y-a.y)*(c.z-a.z)-(b.z-a.z)*(c.y-a.y),  
                     (b.z-a.z)*(c.x-a.x)-(b.x-a.x)*(c.z-a.z),  
                     (b.x-a.x)*(c.y-