这篇博客整理于2016年，主要介绍了计算几何的基础概念和常用算法模板，包括点线段的关系判断、平面几何变换、多边形处理等方面的知识。
【整理于2016年】
计算几何模板
const double eps = 1e-8;
const double pi = acos(-1.0);
using namespace std;
//采用eps的精度判断大/小于零
int epssgn(double x)
{
if (fabs(x)<eps) return 0;
else return x<0?-1:1;
}
//二维点类
struct Point
{
double x,y;
Point(double a=0,double b=0){x=a;y=b;}
};
typedef Point Vector;//矢量可坐标化表示，故用同一个结构体
//二维直线类，一般方程ax+by+c=0
struct Line
{
double a,b,c,angle;
Point p1,p2;
Line(Point s,Point e)
{
a=s.y-e.y;
b=e.x-s.x;
c=s.x*e.y-e.x*s.y;
angle=atan2(e.y-s.y,e.x-s.x);
p1=s;p2=e;
}
Line(){}
};
//二维线段类
struct Segment
{
Point s,e;
Segment(Point a,Point b){s=a;e=b;}
Segment(double x1,double y1,double x2,double y2)
{
s=Point(x1,y1);
e=Point(x2,y2);
}
Segment(){}
};
//向量的加减及数乘
Vector operator + (Point a,Point b)
{
return Vector(a.x+b.x,a.y+b.y);
}
Vector operator - (Point a,Point b)
{
return Vector(a.x-b.x,a.y-b.y);
}
Vector operator * (Point a,double k)
{
return Vector(a.x*k,a.y*k);
}
Vector operator / (Point a,double k)
{
return Vector(a.x/k,a.y/k);
}
//求向量的模（长度）
double len(Vector a)
{
return sqrt(a.x*a.x+a.y*a.y);
}
//得到sp-op和ep-op的叉积
//>0时ep在opsp的逆时针方向，<0时顺时针，=0时共线
double Cross(Point sp, Point ep, Point op)
{
return (sp.x-op.x)*(ep.y-op.y)-(ep.x-op.x)*(sp.y-op.y);
}
//两向量求叉积，求三角形面积需要除以2
double Cross(Vector a,Vector b)
{
return a.x*b.y-b.x*a.y;
}
//两向量求点积
double Dot(Vector a,Vector b)
{
return a.x*b.x+a.y*b.y;
}
//求最大最小值
double max(double a,double b)
{
if (a<b) return b; else return a;
}
double min(double a,double b)
{
if (a>b) return b; else return a;
}
//三点共线
bool Inline(Point a,Point b,Point c)
{
return !sgn(Cross(b,c,a));
}
//点a在线段bc上
bool PointOnSeg(const Point &a, const Point &b, const Point &c)
{
return Inline(b,c,a) && Dot(a-c,a-b)<eps;
}
//求两点之间的直线距离
double dis(Point a,Point b)
{
return sqrt((a.x-b.x)*(a.x-b.x)+(a.y-b.y)*(a.y-b.y));
}
//判断两直线是否平行
int Parallel(Line l1,Line l2)
{
return (fabs(l1.a*l2.b-l2.a*l1.b)<eps);
}
//判断两直线是否相等
int LineEqual(Line l1,Line l2)
{
if (!Parallel(l1,l2)) return 0;
else return (fabs(l1.a*l2.c-l2.a*l1.c)<eps && fabs(l1.b*l2.c-l2.b*l1.c)<eps);
}
//求点C到线段AB的距离
double PointToSegDist(Point A,Point B,Point C)
{
if (dis(A,B)<eps) return dis(B,C);
if (epssgn(Dot(B-A,C-A))<0) return dis(A,C);
if (epssgn(Dot(A-B,C-B))<0) return dis(B,C);
return fabs(Cross(B-A,C-A))/dis(B,A);
}
//求线段两端AB到另一线段CD的距离
double TwoSegMinDist(Point A,Point B,Point C,Point D)
{
return min(min(PointToSegDist(A,B,C),PointToSegDist(A,B,D)),
min(PointToSegDist(C,D,A),PointToSegDist(C,D,B)));
}
Point SymPoint(Point p,Line l) //求二维平面上点p关于直线p1p2的对称点
{
Point result;
double a=l.p2.x-l.p1.x;
double b=l.p2.y-l.p1.y;
double t=((p.x-l.p1.x)*a+(p.y-l.p1.y)*b)/(a*a+b*b);
result.x=2*l.p1.x+2*a*t-p.x;
result.y=2*l.p1.y+2*b*t-p.y;
return result;
}
//直接求两线段ab和cd的交点（确定相交的情况下）
Point LineCross(const Point &a, const Point &b, const Point &c, const Point &d)
{
double u = Cross(b,c,a), v = Cross(a,d,b);
return Point((c.x * v + d.x * u)/(u+v), (c.y * v + d.y * u)/(u+v));
}
//判断线段s1e1与s2e2是否相交（含端点）
//不含端点的话将下面的<=改成<
int IsSegmentIntersect(Point s1, Point e1, Point s2, Point e2)
{
if( min(s1.x,e1.x)<= max(s2.x,e2.x) &&
min(s1.y,e1.y)<= max(s2.y,e2.y) &&
min(s2.x,e2.x)<= max(s1.x,e1.x) &&
min(s2.y,e2.y)<= max(s1.y,e1.y) &&
Cross(s2,e2,s1)*Cross(s2,e2,e1)<=0 &&
Cross(s1,e1,s2)*Cross(s1,e1,e2)<=0)
return 1;
return 0;
}
//知道直线上两点p1p2，判断直线与线段se是否相交，含端点
int IsLineIntersectSegment(Point p1,Point p2,Point s,Point e)
{
if (Cross(p1,p2,s)*Cross(p1,p2,e)>eps) return 0;
else return 1;
}
int IsLineIntersectSegment(Line l1,Point s,Point e)
{
if (Cross(l1.p1,l1.p2,s)*Cross(l1.p1,l1.p2,e)>eps) return 0;
else return 1;
}
//求两条直线l1和l2的交点
Point GetIntersect(Line l1, Line l2)
{
Point res;
res.x=(l1.b*l2.c-l2.b*l1.c)/(l1.a*l2.b-l2.a*l1.b);
res.y=(l1.c*l2.a-l2.c*l1.a)/(l1.a*l2.b-l2.a*l1.b);
return res;
}
//求质量均匀分布的多边形的重心
Point BaryCenter(Point *p,int n)
{
Point res(0,0);
double s=0.0,t;
int i;
for (i=1;i<n-1;i++)
{
t=Cross(p[i]-p[0],p[i+1]-p[0])/2; //分割成三角形，算面积
s+=t;
res.x+=(p[0].x+p[i].x+p[i+1].x)*t; //将面积作为重量放在三角形的重心上
res.y+=(p[0].y+p[i].y+p[i+1].y)*t; //质量*坐标，三角形重心所需的除以3放在后面
}
res.x/=(3*s);
res.y/=(3*s);
return res;
}
//求多边形面积，要求点集按逆时针顺序
double ConvexPolygonArea(Point *p,int n)
{
int i;
double area=0;
for (i=1;i<n-1;i++) area+=Cross(p[i]-p[0],p[i+1]-p[0]);
return area/2;
}
//圆及相关
struct circle
{
Point o;
double r;
};
double dis(circle a, circle b)
{
return sqrt((a.o.x - b.o.x) * (a.o.x - b.o.x) + (a.o.y - b.o.y) * (a.o.y - b.o.y));
}
//两圆相交面积：CircleIntersectArea
double CIA(circle a, circle b)
{
double d = dis(a,b);
if(d >= a.r + b.r) return 0;
if(d <= fabs(a.r - b.r))
{
double r = a.r < b.r ? a.r : b.r;
return pi*r*r;
}
double ang1 = acos((a.r * a.r + d * d - b.r * b.r)/2. / a.r / d );
double ang2 = acos((b.r * b.r + d * d - a.r * a.r)/2. / b.r / d );
double ret = ang1 * a.r * a.r + ang2 * b.r * b.r - d * a.r * sin(ang1);
return ret;
}
//两圆环相交面积，不需其他模板只需一个公式：
//考虑圆环A和B，相交面积一致性地表达为
// S = A大B大 - A小B大 - A大B小 + A小B小
//所以用四遍CIA函数即可
//一般不用考虑距离判断相交，这个公式有一般性
//求以原点为圆心，过a、b，半径为r的优弧扇形面积（即圆心角小于180°）
//如果是以其他点为圆心，事先做一下坐标平移然后再调用
double SectorArea(Point a,Point b,double r)
{
double theta=atan2(a.y,a.x)-atan2(b.y,b.x);
while (theta<=0) theta+=2*pi;
while (theta>=2*pi) theta-=2*pi;
theta=min(theta,2*pi-theta);
return r*r*theta/2;
}
//另一个版本：求以r点为圆心，R为半径，两端顶点为a,b的优弧扇形面积（不用特意平移处理到原点形式了）
double sectorarea(const Point &r, const Point &a, const Point &b,double R)
{
double A2 = Dot(r-a,r-a), B2 = Dot(r-b,r-b), C2 = Dot(a-b,a-b);
return R * R * acos((A2 + B2 - C2)*0.5 / sqrt(A2) / sqrt(B2) ) * 0.5;
}
// 线段与圆相交（判定）
//形参分别为：线段端点a，线段端点b，圆心r，圆半径R，相交点p1,p2
//返回值rtos,目测是r到线段的最小距离
double LineCrossCircle(const Point &a, const Point &b, const Point &r, double R, Point &p1, Point &p2)
{
Point fp = LineCross(r, Point(r.x + a.y - b.y, r.y + b.x -a.x), a, b);
double rtol = dis(r,fp);
double rtos = PointOnSeg(fp,a,b)? rtol : min(dis(a,r),dis(b,r));
double atob = dis(a,b);
double fptoe = sqrt(R * R - rtol * rtol)/atob;
if(rtos > R -eps) return rtos;
p1 = fp + (a-b)*fptoe;
p2 = fp + (b-a)*fptoe;
return rtos;
}
//求以原点为圆心半径为r的园，与原点、a、b组成的三角形的重叠面积
// Circle & Triangle Intersect Area.
// circle : center(0,0),radius=r;
// triangle : a,b,and(0,0);
// maybe you need transform the coordinates to use this function.
double CTIA(Point a,Point b,double r)
{
Point p[2];
int num=0;
int ina=(epssgn(len(a)-r)<0);
int inb=(epssgn(len(b)-r)<0);
if (ina)
{
if (inb) return fabs(Cross(a,b))/2.0; //情形1：三角形完全在圆内：直接求三角形面积
else //情形2：一个点在圆内，一个点在圆外：分割成一个完全在圆内的三角形和一个扇形
{
CircleCrossLine(a,b,Point(0,0),r,p,num);
return SectorArea(b,p[0],r)+fabs(Cross(a,p[0]))/2.0;
}
}
else
{
CircleCrossLine(a,b,Point(0,0),r,p,num);
//情形2：一个点在圆内，一个点在圆外：分割成一个完全在圆内的三角形和一个扇形
if (inb) return SectorArea(p[0],a,r)+fabs(Cross(p[0],b))/2.0;
else
{
//情形4：两个点都在圆外，但是两点的连线与圆有两个焦点：分割成一个完全在圆内的三角形和两个扇形
if (num==2) return SectorArea(a,p[0],r)+SectorArea(p[1],b,r)+fabs(Cross(p[0],p[1]))/2.0;
//情形3：两个点都在圆外，并且相交处为一个扇形：直接求扇形面积
else return SectorArea(a,b,r);
}
}
}
//另一个版本，求以r点为圆心，R为半径的圆 同 三角形arb的相交面积
//意义同上，不再需要预先平移处理了
//这里调用的是免平移版的sectorarea函数
double CTIA(const Point &r,const Point &a, const Point &b, double R)
{
double adis = dis(r,a),bdis = dis(r,b);
if(adis < R + eps && bdis < R +eps)
return Cross(a,b,r)*0.5;
Point ta,tb;
if(Inline(r,a,b)) return 0.0;
double rtos = LineCrossCircle(a,b,r,R,ta,tb);
if(rtos > R - eps)
return sectorarea(r,a,b,R);
if(adis < R + eps)
return Cross(a,tb,r)*0.5 + sectorarea(r,tb,b,R);
if(bdis < R + eps)
return Cross(ta,b,r)*0.5 + sectorarea(r,a,ta,R);
return Cross(ta,tb,r)*0.5 + sectorarea(r,tb,b,R) + sectorarea(r,a, ta, R);
}
//圆与多边形的面积交
//思路是将多边形三角剖分，转化为若干个圆与三角形面积交的问题，再根据叉积性质求和即可，渐进复杂度O(n)，n为多边形规模
// Circle & Polygon Intersect Area.
double CPIA(int n, Point r, double R)
{
int i;
double ret = 0, is_clock_t;
for( i = 0; i < n; i++)
{
is_clock_t = sgn(Cross(p[i],p[(i+1)%n],r));
if( is_clock_t < 0)
{
ret -= CTIA(r, p[(i+1)%n], p[i], R);
}
else ret += CTIA(r, p[i], p[(i+1)%n], R);
}
return fabs(ret);
}
//求两圆公切线（分类讨论了六种可能情形，具体分析见LRJ训练指南）
//形参：圆A，圆B，a,b分别为切点对，即求出一条切线，存下对应的两圆切点
//返回值：公切线数目，-1代表两圆重合，切线无数
//若要求公切线方程 ①若两圆切点不同，直接查询a,b数组，得到两点式 ②若共切点（圆内外切），直接修改本函数，在求出这些切点的时候记下斜率，做点斜式
int getTangents(circle A, circle B, Point *a, Point *b,int &k)
{
int cnt = 0;
if(A.r < B.r)
{
flag = 1;
swap(A,B); swap(a,b);
}
double d = norm(A.o - B.o);
double rdiff = A.r - B.r;
double rsum = A.r + B.r;
if(sgn(d - rdiff)<0) return 0;//两圆内含
double base = angle(B.o - A.o);
if(sgn(d) == 0 && sgn(rdiff) ==0) return -1;//重合
if(sgn(d - rdiff) == 0)//内切，有一条外公切线
{
a[cnt] = A.point(base);
b[cnt] = B.point(base);
cnt ++;
return 1;
}
//普遍意义下的外公切线
//两个if用于判断两圆圆心是否在切线同侧（即外公切线），原版上没有这个判断，但是被某道反演题(hdu4773)卡了
double ang = acos(rdiff/d);
Point m,n;
m = A.point(base + ang);
n = B.point(base + ang);
if(sgn(Cross(m,c[0].o,n)) == sgn(Cross(m,p,n)) && sgn(Cross(m,c[1].o,n)) == sgn(Cross(m,p,n)))
{
a[cnt] = m;
b[cnt] = n;
//add(a[cnt],b[cnt],k);
cnt ++;
}
m = A.point(base - ang);
m = B.point(base - ang);
if(sgn(Cross(m,c[0].o,n)) == sgn(Cross(m,p,n)) && sgn(Cross(m,c[1].o,n)) == sgn(Cross(m,p,n)))
{
a[cnt] = m;
b[cnt] = n;
//add(a[cnt],b[cnt],k);
cnt ++;
}
if(sgn(d - rsum) == 0)//外切，一条内公切线
{
a[cnt] = A.point(base);
b[cnt] = B.point(base + pi);
cnt++;
}
else if(sgn(d-rsum)>0)//外离，两条内公切线
{
double ang_in = acos(rsum/d);
a[cnt] = A.point(base + ang_in);
b[cnt] = B.point(base + ang_in + pi);
cnt ++;
a[cnt] = A.point(base - ang_in);
b[cnt] = B.point(base - ang_in + pi);
cnt++;
}
return cnt;
}
//圆的反演
/*反演的定义：
已知一圆C，圆心为O，半径为r，如果P与P’在过圆心O的直线上，且OP*OP’=r^2，则称P与P’关于O互为反演。O为反演中心，r为反演半径
反演的性质：
除反演中心外，平面上的每一个点都只有唯一的反演点，且这种关系是对称的，位于反演圆上的点，保持在原处，位于反演圆外部的点，变为圆内部的点，位于反演圆内部的点，变为圆外部的点。 举个最简单的例子，区间(0,1]以1为反演半径，那么反演后的区间就是[1,+∞]，这就是一维反演，而圆的反演是二维反演。
任意一条不过反演中心的直线，它的反形是经过反演中心的圆，反之亦然，特别地，过反演中心相交的圆，变为不过反演中心的相交直线
*/
//hdu4773原题代码：求圆，过点p，又与两圆相切。重点看Inverse函数和add函数，分别是把圆反演和把直线反演
//纯粹是数学题，没有复杂度
int sgn(double x)
{
if(fabs(x)<eps)return 0;
if(x > 0)return 1;
else return -1;
}
struct Point
{
double x,y;
Point(double a = 0, double b = 0){x = a; y = b;}
void Input()
{
scanf("%lf%lf",&x,&y);
}
Point trans()
{
return Point(-y, x);
}
void Output()
{
printf("%.8lf %.8lf\n",x,y);
}
};
typedef Point Vector;
double norm(Vector a)
{
return sqrt(a.x*a.x + a.y*a.y);
}
double angle(Vector v)
{
return atan2(v.y, v.x);
}
Vector operator + (Point a, Point b)
{
return Vector(a.x + b.x, a.y + b.y);
}
Vector operator - (Point a, Point b)
{
return Vector(a.x - b.x, a.y - b.y);
}
Point operator * (Point a, double b)
{
return Point(a.x*b, a.y*b);
}
double Cross(Point sp, Point ep, Point op)
{
return (sp.x - op.x)*(ep.y - op.y) - (ep.x - op.x)*(sp.y - op.y);
}
struct circle
{
Point o;
double r;
circle(Point c = Point(0,0), double r = 0):o(c),r(r){}
Point point(double a)
{
return Point(o.x+r*cos(a), o.y+r*sin(a));
}
void Input()
{
o.Input();
scanf("%lf",&r);
}
void Output()
{
printf("%.8lf %.8lf %.8lf\n",o.x,o.y,r);
}
};
int t;
circle c[20];
Point p;
bool flag = 0;
void add(Point a, Point b, int &k)
{
double t = Cross(a,p,b);
if(t < 0) t = -t;
double d = norm(a - b);
t /= d;
if(t > eps)//ab不经过反演中心p
{
double w = 0.5/t;
Point dir = (b-a).trans();
Point a1 = p + dir * (w/d);
Point b1 = p - dir * (w/d);
if(fabs(Cross(a,b,a1)) < fabs(Cross(a,b,b1)))
c[k++] = circle(a1,w);
else c[k++] = circle(b1,w);
}
}
int getTangents(circle A, circle B, Point *a, Point *b,int &k)
{
int cnt = 0;
if(A.r < B.r)
{
flag = 1;
swap(A,B); swap(a,b);
}
double d = norm(A.o - B.o);
double rdiff = A.r - B.r;
double rsum = A.r + B.r;
if(sgn(d - rdiff)<0) return 0;
double base = angle(B.o - A.o);
if(sgn(d) == 0 && sgn(rdiff) ==0) return -1;
if(sgn(d - rdiff) == 0)
{
a[cnt] = A.point(base);
b[cnt] = B.point(base);
cnt ++;
return 1;
}
double ang = acos(rdiff/d);
a[cnt] = A.point(base + ang);
b[cnt] = B.point(base + ang);
if(sgn(Cross(a[cnt],c[0].o,b[cnt])) == sgn(Cross(a[cnt],p,b[cnt])) && sgn(Cross(a[cnt],c[1].o,b[cnt])) == sgn(Cross(a[cnt],p,b[cnt])))
add(a[cnt],b[cnt],k);
cnt ++;
a[cnt] = A.point(base - ang);
b[cnt] = B.point(base - ang);
if(sgn(Cross(a[cnt],c[0].o,b[cnt])) == sgn(Cross(a[cnt],p,b[cnt])) && sgn(Cross(a[cnt],c[1].o,b[cnt])) == sgn(Cross(a[cnt],p,b[cnt])))
add(a[cnt],b[cnt],k);
cnt ++;
if(sgn(d - rsum) == 0)
{
a[cnt] = A.point(base);
b[cnt] = B.point(base + pi);
cnt++;
}
else if(sgn(d-rsum)>0)
{
double ang_in = acos(rsum/d);
a[cnt] = A.point(base + ang_in);
b[cnt] = B.point(base + ang_in + pi);
cnt ++;
a[cnt] = A.point(base - ang_in);
b[cnt] = B.point(base - ang_in + pi);
cnt++;
}
return cnt;
}
circle Inverse(circle C)
{
circle T;
double t = norm(C.o - p);
double x = 1.0 / (t - C.r);
double y = 1.0 / (t + C.r);
T.r = (x - y)/2.0;
double s = (x + y)/2.0;
T.o = p + (C.o - p)*(s/t);
return T;
}
int solve()
{
c[0] = Inverse(c[0]);
c[1] = Inverse(c[1]);
int k = 2;
Point tmpa[20];
Point tmpb[20];
int tt = getTangents(c[0],c[1],tmpa,tmpb,k);
return k-2;
}
int main()
{
scanf("%d",&t);
while(t--)
{
c[0].Input();c[1].Input();
p.Input();
int num = solve();
printf("%d\n",num);
for(int i = 0; i < num ; i++)
c[i+2].Output();
}
}
//三角形相关
//求质量均匀分布的三角形的重心
Point TriangleMassCenter(Point a,Point b,Point c)
{
return (a+b+c)/3.0;
}
//求三角形的外心
Point TriangleCircumCenter(Point a,Point b,Point c)
{
Point res;
double a1=atan2(b.y-a.y,b.x-a.x)+pi/2;
double a2=atan2(c.y-b.y,c.x-b.x)+pi/2;
double ax=(a.x+b.x)/2;
double ay=(a.y+b.y)/2;
double bx=(c.x+b.x)/2;
double by=(c.y+b.y)/2;
double r1=(sin(a2)*(ax-bx)+cos(a2)*(by-ay))/(sin(a1)*cos(a2)-sin(a2)*cos(a1));
return Point(ax+r1*cos(a1),ay+r1*sin(a1));
}
//求三角形的垂心
Point TriangleOrthoCenter(Point a,Point b,Point c)
{
return TriangleMassCenter(a,b,c)*3.0-TriangleCircumCenter(a,b,c)*2.0;
}
//求三角形的内心
Point TriangleInnerCenter(Point a,Point b,Point c)
{
Point res;
double la=len(b-c);
double lb=len(a-c);
double lc=len(a-b);
res.x=(la*a.x+lb*b.x+lc*c.x)/(la+lb+lc);
res.y=(la*a.y+lb*b.y+lc*c.y)/(la+lb+lc);
}
/*------Graham求凸包----*/
int cmp(Point &a,Point &b)
{
if (a.y==b.y) return (a.x<b.x);
return (a.y<b.y);
}
//对所有的点进行一次排序
void QSort(Point p[],int l,int r)
{
int i=l,j=r;
Point mid=p[(l+r)/2],swap;
while (i<=j)
{
while (cmp(p[i],mid)) i++;
while (cmp(mid,p[j])) j--;
if (i<=j)
{
swap=p[i];
p[i]=p[j];
p[j]=swap;
i++;j--;
}
}
if (i<r) QSort(p,i,r);
if (l<j) QSort(p,l,j);
}
//n为原图的点数，p[]为原图的点（0~n-1），top为凸包点的数量（0~top-1），res[]为凸包点的下标，。
int Graham(Point p[],int n,int res[])
{
int i,len,top;
top=1;
QSort(p,0,n-1);
for (i=0;i<3;i++) res[i]=i;
for (i=2;i<n;i++)
{
while (top && epssgn(Cross(p[i],p[res[top]],p[res[top-1]]))>=0) top--;
res[++top]=i;
}
len=top;
res[++top]=n-2;
for (i=n-3;i>=0;i--)
{
while (top!=len && epssgn(Cross(p[i], p[res[top]], p[res[top-1]]))>=0) top--;
res[++top]=i;
}
return top;
}
//下面是凸包相关功能函数
//判断点x是否在凸包p[res[1~chnum]]中，返回1则为在内部或边界上
int PointInConvexHull(Point p[],int res[],int chnum,Point x)
{
//找一个凸包内的点
Point g=(p[res[0]]+p[res[chnum/3]]+p[res[2*chnum/3]])/3.0;
int l=0,r=chnum,mid;
//二分凸包
while (l+1<r)
{
mid=(l+r)/2;
if (epssgn(Cross(p[res[l]]-g,p[res[mid]]-g))>0)
{
if (epssgn(Cross(p[res[l]]-g,x-g))>=0
&& epssgn(Cross(p[res[mid]]-g,x-g))<0) r=mid;
else l=mid;
}
else
{
if (epssgn(Cross(p[res[l]]-g,x-g))<0
&& epssgn(Cross(p[res[mid]]-g,x-g))>=0) l=mid;
else r=mid;
}
}
r%=chnum;
if (epssgn(Cross(p[res[r]]-x,p[res[l]]-x))==-1) return 1; else return 0;
}
//旋转卡壳求凸包p[res[1~chnum]]的直径，对踵点数量appnum，存储在app[][2]中
double Diameter(Point p[],int res[],int chnum,int app[][2],int &appnum) //AntiPodal Point
{
int i,j;
double ret=0,nowlen;
res[chnum]=res[0];
j=1;
appnum=0;
for (i=0;i<chnum;i++)
{
while (Cross(p[res[i]]-p[res[i+1]],p[res[j+1]]-p[res[i+1]])
<Cross(p[res[i]]-p[res[i+1]],p[res[j]]-p[res[i+1]]))
{
j++;
j%=chnum;
}
app[appnum][0]=res[i];
app[appnum][1]=res[j];
appnum++;
nowlen=dis(p[res[i]],p[res[j]]);
if (nowlen>ret) ret=nowlen;
nowlen=dis(p[res[i+1]],p[res[j+1]]);
if (nowlen>ret) ret=nowlen;
}
return ret;
}
//旋转卡壳求凸包p[res[1~chnum]]内最大面积三角形
double ConvexHullMaxTriangleArea(Point p[],int res[],int chnum)
{
int i,j,k;
double area=0,tmp;
k=2;j=1;
res[chnum]=res[0];
for (i=0;i<chnum;i++)
{
//卡住i,j，若向前旋转k能令面积更大的话就一直转一下去
while (fabs(Cross(p[res[j]]-p[res[i]],p[res[(k+1)%chnum]]-p[res[i]]))
>fabs(Cross(p[res[j]]-p[res[i]],p[res[k]] -p[res[i]])))
k=(k+1)%chnum;
tmp=fabs(Cross(p[res[j]]-p[res[i]],p[res[k]]-p[res[i]]));
if (tmp>area) area=tmp;
//卡住i,k，若向前旋转j能令面积更大的话就一直转一下去
while (fabs(Cross(p[res[(j+1)%chnum]]-p[res[i]],p[res[k]]-p[res[i]]))
>fabs(Cross(p[res[j]]- p[res[i]],p[res[k]]-p[res[i]])))
j=(j+1)%chnum;
tmp=fabs(Cross(p[res[j]]-p[res[i]],p[res[k]]-p[res[i]]));
if (tmp>area) area=tmp;
//j,k转到最大位置，向前转i(i++)
}
return area/2;
}
//求两凸包p、q间的最小距离，注意调用的时候要交换参数位置调用两次
//卡住p（p步进），旋转q，求能得到的最小距离
double TwoConvexHullMinDist(Point P[],Point Q[],int n,int m)
{
int YMinP=0,YMaxQ=0,i;
double tmp,ans=999999999;
for (i=0;i<n;i++) if (P[i].y<P[YMinP].y) YMinP=i;
for (i=0;i<m;i++) if (Q[i].y>Q[YMaxQ].y) YMaxQ=i;
P[n]=P[0];
Q[m]=Q[0];
for (i=0;i<n;i++)
{
//注意，tmp保存的是>比较的结果
while (tmp=Cross(Q[YMaxQ+1]-P[YMinP+1],P[YMinP]-P[YMinP+1])
>Cross(Q[YMaxQ] -P[YMinP+1],P[YMinP]-P[YMinP+1]))
YMaxQ=(YMaxQ+1)%m;
if (tmp<0) ans=min(ans,PointToSegDist(P[YMinP],P[YMinP+1],Q[YMaxQ]));
else ans=min(ans,TwoSegMinDist (P[YMinP],P[YMinP+1],Q[YMaxQ],Q[YMaxQ+1]));
YMinP=(YMinP+1)%n;
}
return ans;
}
/*------Graham求凸包----*/
/*-----半平面交O(nlogn)模板-----*/
int cmp(const Line & l1,const Line & l2)
{
int d=epssgn(l1.angle-l2.angle);
if (!d) return (epssgn(Cross(l2.p1-l1.p1,l2.p2-l1.p1))>0); //极角相同时，将更靠半平面里面的放在前面
return d<0;
}
void QSort(Line L[],int l,int r)
{
int i=l,j=r;
Line swap,mid=L[(l+r)/2];
while (i<=j)
{
while (cmp(L[i],mid)) i++;
while (cmp(mid,L[j])) j--;
if (i<=j)
{
swap=L[i];
L[i]=L[j];
L[j]=swap;
i++;j--;
}
}
if (i<r) QSort(L,i,r);
if (l<j) QSort(L,l,j);
}
//判断l1与l2的交点是否在半平面hpl外
int IntersectionOutOfHalfPlane(Line &hpl,Line &l1,Line &l2)
{
Point p=GetIntersect(l1,l2);
return (epssgn(Cross(hpl.p1-p,hpl.p2-p))<0);
}
//半平面hpl向内推进dis长度
Line HalfPlaneMoveIn(Line &hpl,double &dis)
{
double dx=hpl.p1.x-hpl.p2.x;
double dy=hpl.p1.y-hpl.p2.y;
double ll=len(hpl.p1-hpl.p2);
Point pa=Point(dis*dy/ll+hpl.p1.x,hpl.p1.y-dis*dx/ll);
Point pb=Point(dis*dy/ll+hpl.p2.x,hpl.p2.y-dis*dx/ll);
return Line(pa,pb);
}
//求n个半平面l的半平面交，得到的交点储存在p中，交点数目返回到pn
//可以将一个多边形每一条边看成半平面，求出来的交就是多边形的核，要求pn>=3
void HalfPlaneIntersect(Line l[],int n,Point p[],int &pn)
{
int i,j;
int dq[MaxN],top,bot;
//排序是在满足所有半平面A*x+B*y+C>0或（<,<=,>=）,
//也就是所有半平面的符号均相同的情况下对极角进行排序。
QSort(l,0,n-1);
//极角相同时，只保留最靠里面的那条
for (i=j=0;i<n;i++) if (epssgn(l[i].angle-l[j].angle)>0) l[++j]=l[i];
n=j+1;
dq[0]=0; //双端队列
dq[1]=1;
top=1; //顶部和底部
bot=0;
for (i=2;i<n;i++)
{
//当栈顶的两条直线交点在当前半平面外部时，弹栈
while (top>bot && IntersectionOutOfHalfPlane(l[i],l[dq[top]],l[dq[top-1]])) top--;
//由于求的是一个凸多边形，所以当半平面转过接近一圈时，某个半平面满足上一个while的条件后，
//它又会影响到底部的两条直线，当底部的两条直线的交点，在当前的半平面外部时，底部弹栈
while (top>bot && IntersectionOutOfHalfPlane(l[i],l[dq[bot]],l[dq[bot+1]])) bot++;
dq[++top]=i; //当前半平面入栈
}
//当最顶部的两条直线的交点不在最底部的半平面内时，顶部的那个半平面是多余的，顶部弹栈
while (top>bot && IntersectionOutOfHalfPlane(l[dq[bot]],l[dq[top]],l[dq[top-1]])) top--;
//当最底部的两条直线的交点不在最顶部的半平面内时，底部的那个半平面是多余的，底部弹栈
while (top>bot && IntersectionOutOfHalfPlane(l[dq[top]],l[dq[bot]],l[dq[bot+1]])) bot++;
dq[++top]=dq[bot]; //将最底部的半平面放到最顶部来，方便下面求顶点
for (pn=0,i=bot;i<top;i++,pn++) p[pn]=GetIntersect(l[dq[i+1]],l[dq[i]]);
}
/*-----半平面交O(nlogn)模板-----*/
/*-----半平面交O(n^2)模板-----*/
//p是现在切出来的半平面的点，pn是点的数量，需要按顺时针或者逆时针排序
//新的半平面直线为ax+by+c=0
void HalfPlaneCut(Point p[],int &pn,double a,double b,double c)
{
int cnt=0,i;
Point tp[MaxN]; //temp_p
//现在交出来的点在新的半平面内，保留
for (i=1;i<=pn;i++) if (epssgn(a*p[i].x+b*p[i].y+c)>=0) tp[++cnt]=p[i];
else //否则如果其前后的点在半平面内，重新切割
{
if (epssgn(a*p[i-1].x+b*p[i-1].y+c)>0)
tp[++cnt]=GetIntersect(Line(p[i-1],p[i]),Line(a,b,c));
if (epssgn(a*p[i+1].x+b*p[i+1].y+c)>0)
tp[++cnt]=GetIntersect(Line(p[i],p[i+1]),Line(a,b,c));
}
pn=cnt;
for (i=1;i<=cnt;i++) p[i]=tp[i];
p[0]=p[cnt];
p[cnt+1]=p[1];
}
/*-----半平面交O(n^2)模板-----*/
/*
* *
* 三维部分 *
* *
*/
/***********基础*************/
int dcmp(double x)
{
if(fabs(x) < eps) return 0;
if(x < 0)return -1;
return 1;
}
struct Point3
{
double x, y, z;
Point3(double x=0, double y=0, double z=0):x(x),y(y),z(z) { }
};
typedef Point3 Vector3;
Vector3 operator + (const Vector3& A, const Vector3& B)
{
return Vector3(A.x+B.x, A.y+B.y, A.z+B.z);
}
Vector3 operator - (const Point3& A, const Point3& B)
{
return Vector3(A.x-B.x, A.y-B.y, A.z-B.z);
}
Vector3 operator * (const Vector3& A, double p)
{
return Vector3(A.x*p, A.y*p, A.z*p);
}
Vector3 operator / (const Vector3& A, double p)
{
return Vector3(A.x/p, A.y/p, A.z/p);
}
double Dot(const Vector3& A, const Vector3& B)
{
return A.x*B.x + A.y*B.y + A.z*B.z;
}
double norm(const Vector3& A)
{
return sqrt(Dot(A, A));
}
double Angle(const Vector3& A, const Vector3& B)
{
return acos(Dot(A, B) / norm(A) / norm(B));
}
Vector3 Cross(const Vector3& A, const Vector3& B)
{
return Vector3(A.y*B.z - A.z*B.y, A.z*B.x - A.x*B.z, A.x*B.y - A.y*B.x);
}
//空间三角形面积的2倍，只需求面积的话除以2.0就可以了
double Area2(const Point3& A, const Point3& B, const Point3& C)
{
return norm(Cross(B-A, C-A));
}
//空间四点形成的四面体的体积的6倍，只需求体积的话除以6.0就可以了
double Volume6(const Point3& A, const Point3& B, const Point3& C, const Point3& D)
{
return Dot(D-A, Cross(B-A, C-A));
}
// 四面体的重心
Point3 Centroid(const Point3& A, const Point3& B, const Point3& C, const Point3& D)
{
return (A + B + C + D)/4.0;
}
/************三维点线面*************/
// 点p到平面p0-n的距离
// n是平面法向量，且必须为单位向量
// n不是单位向量的话事先做单位化即可
double DistanceToPlane(const Point3& p, const Point3& p0, const Vector3& n)
{
return fabs(Dot(p-p0, n)); // 如果不取绝对值，得到的是有向距离
}
// 点p在平面p0-n上的投影
// n是平面法向量，且必须为单位向量
// n不是单位向量的话事先做单位化即可
Point3 GetPlaneProjection(const Point3& p, const Point3& p0, const Vector3& n)
{
return p-n*Dot(p-p0, n);
}
//直线p1-p2 与平面p0-n的交点
//内含直线平行于平面和重合的情况，按需改写即可
Point3 LinePlaneIntersection(Point3 p1, Point3 p2, Point3 p0, Vector3 n)
{
vector3 v = p2-p1;
double t = (Dot(n, p0-p1) / Dot(n, p2-p1));//分母为0，直线与平面平行或在平面上
return p1 + v*t; //如果是线段 判断t是否在0~1之间
}
// 点P到直线AB的距离
double DistanceToLine(const Point3& P, const Point3& A, const Point3& B)
{
Vector3 v1 = B - A, v2 = P - A;
return norm(Cross(v1, v2)) / norm(v1);
}
//点p到线段ab的距离
double DistanceToSeg(Point3 p, Point3 a, Point3 b)
{
if(a == b) return norm(p-a);
Vector3 v1 = b-a, v2 = p-a, v3 = p-b;
if(dcmp(Dot(v1, v2)) < 0) return norm(v2);
else if(dcmp(Dot(v1, v3)) > 0) return norm(v3);
else return norm(Cross(v1, v2)) / norm(v1);
}
//求异面直线 p1+s*u与p2+t*v的公垂线对应的s 如果平行|重合，返回false
//不懂。。
bool LineDistance3D(Point3 p1, Vector3 u, Point3 p2, Vector3 v, double& s)
{
double b = Dot(u, u) * Dot(v, v) - Dot(u, v) * Dot(u, v);
if(dcmp(b) == 0) return false;
double a = Dot(u, v) * Dot(v, p1-p2) - Dot(v, v) * Dot(u, p1-p2);
s = a/b;
return true;
}
// p1和p2是否在线段a-b的同侧
bool SameSide(const Point3& p1, const Point3& p2, const Point3& a, const Point3& b)
{
return dcmp(Dot(Cross(b-a, p1-a), Cross(b-a, p2-a))) >= 0;
}
// 点P在三角形P0, P1, P2中
bool PointInTri(const Point3& P, const Point3& P0, const Point3& P1, const Point3& P2)
{
return SameSide(P, P0, P1, P2) && SameSide(P, P1, P0, P2) && SameSide(P, P2, P0, P1);
}
// 三角形P0P1P2是否和线段AB相交
bool TriSegIntersection(const Point3& P0, const Point3& P1, const Point3& P2, const Point3& A, const Point3& B, Point3& P)
{
Vector3 n = Cross(P1-P0, P2-P0);
if(dcmp(Dot(n, B-A)) == 0) return false; // 线段A-B和平面P0P1P2平行或共面
else // 平面A和直线P1-P2有惟一交点
{
double t = Dot(n, P0-A) / Dot(n, B-A);
if(dcmp(t) < 0 || dcmp(t-1) > 0) return false; // 不在线段AB上
P = A + (B-A)*t; // 交点
return PointInTri(P, P0, P1, P2);
}
}
//空间两三角形是否相交
bool TriTriIntersection(Point3* T1, Point3* T2)
{
Point3 P;
for(int i = 0; i < 3; i++)
{
if(TriSegIntersection(T1[0], T1[1], T1[2], T2[i], T2[(i+1)%3], P)) return true;
if(TriSegIntersection(T2[0], T2[1], T2[2], T1[i], T1[(i+1)%3], P)) return true;
}
return false;
}
//空间两直线上最近点对 返回最近距离 点对保存在ans1 ans2中
double SegSegDistance(Point3 a1, Point3 b1, Point3 a2, Point b2, Point& ans1, Point& ans2)
{
Vector v1 = (a1-b1), v2 = (a2-b2);
Vector N = Cross(v1, v2);
Vector ab = (a1-a2);
double ans = Dot(N, ab) / norm(N);
Point p1 = a1, p2 = a2;
Vector d1 = b1-a1, d2 = b2-a2;
double t1, t2;
t1 = Dot((Cross(p2-p1, d2)), Cross(d1, d2));
t2 = Dot((Cross(p2-p1, d1)), Cross(d1, d2));
double dd = norm((Cross(d1, d2)));
t1 /= dd*dd;
t2 /= dd*dd;
ans1 = (a1 + (b1-a1)*t1);
ans2 = (a2 + (b2-a2)*t2);
return fabs(ans);
}
// 判断P是否在三角形A, B, C中，并且到三条边的距离都至少为mindist。保证P, A, B, C共面
bool InsideWithMinDistance(const Point3& P, const Point3& A, const Point3& B, const Point3& C, double mindist)
{
if(!PointInTri(P, A, B, C)) return false;
if(DistanceToLine(P, A, B) < mindist) return false;
if(DistanceToLine(P, B, C) < mindist) return false;
if(DistanceToLine(P, C, A) < mindist) return false;
return true;
}
// 判断P是否在凸四边形ABCD（顺时针或逆时针）中，并且到四条边的距离都至少为mindist。保证P, A, B, C, D共面
bool InsideWithMinDistance(const Point3& P, const Point3& A, const Point3& B, const Point3& C, const Point3& D, double mindist)
{
if(!PointInTri(P, A, B, C)) return false;
if(!PointInTri(P, C, D, A)) return false;
if(DistanceToLine(P, A, B) < mindist) return false;
if(DistanceToLine(P, B, C) < mindist) return false;
if(DistanceToLine(P, C, D) < mindist) return false;
if(DistanceToLine(P, D, A) < mindist) return false;
return true;
}
/*************凸包相关问题*******************/
//加干扰
double rand01()
{
return rand() / (double)RAND_MAX;
}
double randeps()
{
return (rand01() - 0.5) * eps;
}
Point3 add_noise(const Point3& p)
{
return Point3(p.x + randeps(), p.y + randeps(), p.z + randeps());
}
struct Face
{
int v[3];
Face(int a, int b, int c)
{
v[0] = a;
v[1] = b;
v[2] = c;
}
Vector3 Normal(const vector<Point3>& P) const
{
return Cross(P[v[1]]-P[v[0]], P[v[2]]-P[v[0]]);
}
// f是否能看见P[i]
int CanSee(const vector<Point3>& P, int i) const
{
return Dot(P[i]-P[v[0]], Normal(P)) > 0;
}
};
// 增量法求三维凸包
// 注意：没有考虑各种特殊情况（如四点共面）。实践中，请在调用前对输入点进行微小扰动
vector<Face> CH3D(const vector<Point3>& P)
{
int n = P.size();
vector<vector<int> > vis(n);
for(int i = 0; i < n; i++) vis[i].resize(n);
vector<Face> cur;
cur.push_back(Face(0, 1, 2)); // 由于已经进行扰动，前三个点不共线
cur.push_back(Face(2, 1, 0));
for(int i = 3; i < n; i++)
{
vector<Face> next;
// 计算每条边的“左面”的可见性
for(int j = 0; j < cur.size(); j++)
{
Face& f = cur[j];
int res = f.CanSee(P, i);
if(!res) next.push_back(f);
for(int k = 0; k < 3; k++) vis[f.v[k]][f.v[(k+1)%3]] = res;
}
for(int j = 0; j < cur.size(); j++)
for(int k = 0; k < 3; k++)
{
int a = cur[j].v[k], b = cur[j].v[(k+1)%3];
if(vis[a][b] != vis[b][a] && vis[a][b]) // (a,b)是分界线，左边对P[i]可见
next.push_back(Face(a, b, i));
}
cur = next;
}
return cur;
}
struct ConvexPolyhedron
{
int n;
vector<Point3> P, P2;
vector<Face> faces;
bool read()
{
if(scanf("%d", &n) != 1) return false;
P.resize(n);
P2.resize(n);
for(int i = 0; i < n; i++)
{
P[i] = read_point3();
P2[i] = add_noise(P[i]);
}
faces = CH3D(P2);
return true;
}
//三维凸包重心
Point3 centroid()
{
Point3 C = P[0];
double totv = 0;
Point3 tot(0,0,0);
for(int i = 0; i < faces.size(); i++)
{
Point3 p1 = P[faces[i].v[0]], p2 = P[faces[i].v[1]], p3 = P[faces[i].v[2]];
double v = -Volume6(p1, p2, p3, C);
totv += v;
tot = tot + Centroid(p1, p2, p3, C)*v;
}
return tot / totv;
}
//凸包重心到表面最近距离
double mindist(Point3 C)
{
double ans = 1e30;
for(int i = 0; i < faces.size(); i++)
{
Point3 p1 = P[faces[i].v[0]], p2 = P[faces[i].v[1]], p3 = P[faces[i].v[2]];
ans = min(ans, fabs(-Volume6(p1, p2, p3, C) / Area2(p1, p2, p3)));
}
return ans;
}
};
//然后我发现以上三维凸包模板没有写别的函数，比如求凸包表面多边形数目等，我对比了这个模板和网上更全一点但是没有距离函数的模板，发现统一起来的难度蛮大。在此先贴上把后者补充进来，以备所需。
//下面这套模板同样是增量法求三维凸包，原理同上，不过实现方式有所区别，提供的功能侧重于凸包外表的几何性质：凸包表面积，凸包体积，凸包表面多边形个数等
#define PR (1e-8)
//三维点结构
struct TPoint
{
double x,y,z;
TPoint() {}
TPoint(double _x,double _y,double _z):x(_x),y(_y),z(_z) {}
TPoint operator - (const TPoint p)
{
return TPoint(x-p.x,y-p.y,z-p.z);
}
TPoint operator * (const TPoint p)
{
return TPoint(y*p.z-z*p.y,z*p.x-x*p.z,x*p.y-y*p.x);
}//叉积
double operator ^ (const TPoint p)
{
return x*p.x+y*p.y+z*p.z;
}//点积
};
struct fac
{
int a,b,c;//凸包一个面上的三个点的编号
bool ok;//该面是否是最终凸包中的面
};
struct T3dull
{
int n;//初始点数
TPoint ply[N];//初始点
int trianglecnt;//凸包上三角形数
fac tri[N];//凸包三角形
int vis[N][N];//点i到点j是属于哪个面
//两点长度
double dist(TPoint a)
{
return sqrt(a.x*a.x+a.y*a.y+a.z*a.z);
}
//三角形面积*2
double area(TPoint a,TPoint b,TPoint c)
{
return dist((b-a)*(c-a));
}
//四面体有向体积*6
double volume(TPoint a,TPoint b,TPoint c,TPoint d)
{
return (b-a)*(c-a)^(d-a);
}
//点p到面f距离，返回值为正的几何意义是p在f同侧
//这个函数的返回值还没有完全搞懂，仅限模板引用
double ptoplane(TPoint &p,fac &f)
{
TPoint m=ply[f.b]-ply[f.a],n=ply[f.c]-ply[f.a],t=p-ply[f.a];
return (m*n)^t;
}
void deal(int p,int a,int b)
{
int f=vis[a][b];
fac add;
if(tri[f].ok)
{
if((ptoplane(ply[p],tri[f]))>PR) dfs(p,f);
else
{
add.a=b,add.b=a,add.c=p,add.ok=1;
vis[p][b]=vis[a][p]=vis[b][a]=trianglecnt;
tri[trianglecnt++]=add;
}
}
}
//维护凸包，如果点p在凸包外更新凸包
void dfs(int p,int cnt)
{
tri[cnt].ok=0;
deal(p,tri[cnt].b,tri[cnt].a);
deal(p,tri[cnt].c,tri[cnt].b);
deal(p,tri[cnt].a,tri[cnt].c);
}
//判断两个面是否为同一面
bool same(int s,int e)
{
TPoint a=ply[tri[s].a],b=ply[tri[s].b],c=ply[tri[s].c];
return fabs(volume(a,b,c,ply[tri[e].a]))<PR&&fabs(volume(a,b,c,ply[tri[e].b]))<PR&&fabs(volume(a,b,c,ply[tri[e].c]))<PR;
}
//构建凸包，考虑了特殊情况
void construct()
{
int i,j;
trianglecnt=0;
if(n<4) return ;
bool tmp=1;
for(int i=1; i<n; i++) //前两点不共点
{
if((dist(ply[0]-ply[i]))>PR)
{
swap(ply[1],ply[i]);
tmp=0;
break;
}
}
if(tmp) return ;
tmp=1;
for(i=2; i<n; i++) //前三点不共线
{
if((dist((ply[0]-ply[1])*(ply[1]-ply[i])))>PR)
{
swap(ply[2],ply[i]);
tmp=0;
break;
}
}
if(tmp) return ;
tmp=1;
for(i=3; i<n; i++) //前四点不共面
{
if(fabs((ply[0]-ply[1])*(ply[1]-ply[2])^(ply[0]-ply[i]))>PR)
{
swap(ply[3],ply[i]);
tmp=0;
break;
}
}
if(tmp) return ;
fac add;
for(int i=0; i<4; i++) //构建初始四面体
{
add.a=(i+1)%4,add.b=(i+2)%4,add.c=(i+3)%4,add.ok=1;
if((ptoplane(ply[i],add))>0)
swap(add.b,add.c);
vis[add.a][add.b]=vis[add.b][add.c]=vis[add.c][add.a]=trianglecnt;
tri[trianglecnt++]=add;
}
for(int i=4; i<n; i++) //构建更新凸包
{
for(j=0; j<trianglecnt; j++)
{
if(tri[j].ok&&(ptoplane(ply[i],tri[j]))>PR)
{
dfs(i,j);
break;
}
}
}
int cnt=trianglecnt;
trianglecnt=0;
for(i=0; i<cnt; i++)
{
if(tri[i].ok)
tri[trianglecnt++]=tri[i];
}
}
//凸包表面积
double area()
{
double ret=0;
for(int i=0; i<trianglecnt; i++)
ret+=area(ply[tri[i].a],ply[tri[i].b],ply[tri[i].c]);
return ret/2.0;
}
//凸包体积
double volume()
{
TPoint p(0,0,0);
double ret=0;
for(int i=0; i<trianglecnt; i++)
ret+=volume(p,ply[tri[i].a],ply[tri[i].b],ply[tri[i].c]);
return fabs(ret/6);
}
//凸包表面三角形数
int facetri()
{
return trianglecnt;
}
//表面多边形数
int facepolygon()
{
int ans=0,i,j,k;
for(i=0; i<trianglecnt; i++)
{
for(j=0,k=1; j<i; j++)
{
if(same(i,j))
{
k=0;
break;
}
}
ans+=k;
}
return ans;
}
} hull;
扫一扫
1355
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
