点、线、面相关的算法（2）

可用射线法来判断一个点是否在多边形的内部： 

射线法就是以这个要判断的点作一射线(为了方便，直接作一水平向右的射线)，数一下线段与 
多边形边的交点数，奇数时就是在多边形内，偶数时就是在多边形外。 

/* 判断线段是否在简单多边形内(注意：如果多边形是凸多边形，下面的算法可以化简) 
     原理： 
必要条件一：线段的两个端点都在多边形内； 
必要条件二：线段和多边形的所有边都不内交； 
用途：1. 判断折线是否在简单多边形内 
        2. 判断简单多边形是否在另一个简单多边形内 
*/ 
bool LinesegInsidePolygon(int vcount,POINT polygon[],LINESEG l) 
{ 
// 判断线端l的端点是否不都在多边形内 
if(!insidepolygon(vcount,polygon,l.s)||!insidepolygon(vcount,polygon,l.e)) 
return false; 
int top=0,i,j; 
POINT PointSet[MAXV],tmp; 
LINESEG s; 

for(i=0;i<vcount;i++) 
{ 
s.s=polygon[i]; 
s.e=polygon[(i+1)%vcount]; 
if(online(s,l.s)) //线段l的起始端点在线段s上 
PointSet[top++]=l.s; 
else if(online(s,l.e)) //线段l的终止端点在线段s上 
PointSet[top++]=l.e; 
else 
{ 
if(online(l,s.s)) //线段s的起始端点在线段l上 
PointSet[top++]=s.s; 
else if(online(l,s.e)) // 线段s的终止端点在线段l上 
PointSet[top++]=s.e; 
else 
{ 
if(intersect(l,s)) // 这个时候如果相交，肯定是内交，返回false 
return false; 
} 
} 
} 

for(i=0;i<top-1;i++) /* 冒泡排序，x坐标小的排在前面；x坐标相同者， 
y坐标小的排在前面 */ 
{ 
for(j=i+1;j<top;j++) 
{ 
if( PointSet[i].x>PointSet[j].x || fabs(PointSet[i].x-PointSet[j].x)<EP && PointSet[i].y>PointSet[j].y ) 
     { 
tmp=PointSet[i]; 
PointSet[i]=PointSet[j]; 
PointSet[j]=tmp; 
} 
} 
} 

for(i=0;i<top-1;i++) 
{ 
tmp.x=(PointSet[i].x+PointSet[i+1].x)/2; //得到两个相邻交点的中点 
tmp.y=(PointSet[i].y+PointSet[i+1].y)/2; 
if(!insidepolygon(vcount,polygon,tmp)) 
return false; 
} 
return true; 
} 

/* 求任意简单多边形polygon的重心 
需要调用下面几个函数： 
void AddPosPart(); 增加右边区域的面积 
void AddNegPart(); 增加左边区域的面积 
void AddRegion(); 增加区域面积 
在使用该程序时，如果把xtr,ytr,wtr,xtl,ytl,wtl设成全局变量就可以使这些函数的形式得到化简,但要注意函数的声明和调用要做相应变化 
*/ 
void AddPosPart(double x, double y, double w, double &xtr, double &ytr, double &wtr) 
{ 
if (abs(wtr + w)<1e-10 ) return; // detect zero regions 
xtr = ( wtr*xtr + w*x ) / ( wtr + w ); 
ytr = ( wtr*ytr + w*y ) / ( wtr + w ); 
wtr = w + wtr; 
return; 
} 
void AddNegPart(double x, ouble y, double w, double &xtl, double &ytl, double &wtl) 
{ 
if ( abs(wtl + w)<1e-10 ) 
return; // detect zero regions 

xtl = ( wtl*xtl + w*x ) / ( wtl + w ); 
ytl = ( wtl*ytl + w*y ) / ( wtl + w ); 
wtl = w + wtl; 
return; 
} 

void AddRegion ( double x1, double y1, double x2, double y2, double &xtr, double &ytr, 
double &wtr, double &xtl, double &ytl, double &wtl ) 
{ 
if ( abs (x1 - x2)< 1e-10 ) 
return; 

if ( x2 > x1 ) 
     { 
AddPosPart ((x2+x1)/2, y1/2, (x2-x1) * y1,xtr,ytr,wtr); /* rectangle 全局 
变量变化处 */ 
AddPosPart ((x1+x2+x2)/3, (y1+y1+y2)/3, (x2-x1)*(y2-y1)/2,xtr,ytr,wtr);    
// triangle 全局变量变化处 
     } 
     else 
     { 
AddNegPart ((x2+x1)/2, y1/2, (x2-x1) * y1,xtl,ytl,wtl);  
// rectangle 全局变量变化处 
AddNegPart ((x1+x2+x2)/3, (y1+y1+y2)/3, (x2-x1)*(y2-y1)/2,xtl,ytl,wtl);  
// triangle    全局变量变化处 
     } 
} 

POINT cg_simple(int vcount,POINT polygon[]) 
{ 
     double xtr,ytr,wtr,xtl,ytl,wtl;        
//注意： 如果把xtr,ytr,wtr,xtl,ytl,wtl改成全局变量后这里要删去 
     POINT p1,p2,tp; 
     xtr = ytr = wtr = 0.0; 
     xtl = ytl = wtl = 0.0; 
     for(int i=0;i<vcount;i++) 
{ 
p1=polygon[i]; 
p2=polygon[(i+1)%vcount]; 
AddRegion(p1.x,p1.y,p2.x,p2.y,xtr,ytr,wtr,xtl,ytl,wtl); //全局变量变化处 
} 
tp.x = (wtr*xtr + wtl*xtl) / (wtr + wtl); 
tp.y = (wtr*ytr + wtl*ytl) / (wtr + wtl); 
return tp; 
} 

// 求凸多边形的重心,要求输入多边形按逆时针排序 
POINT gravitycenter(int vcount,POINT polygon[]) 
{ 
POINT tp; 
double x,y,s,x0,y0,cs,k; 
x=0;y=0;s=0; 
for(int i=1;i<vcount-1;i++) 
{ 
x0=(polygon[0].x+polygon[i].x+polygon[i+1].x)/3; 
y0=(polygon[0].y+polygon[i].y+polygon[i+1].y)/3; //求当前三角形的重心 
cs=multiply(polygon[i],polygon[i+1],polygon[0])/2; 
//三角形面积可以直接利用该公式求解 
if(abs(s)<1e-20) 
{ 
x=x0;y=y0;s+=cs;continue; 
} 
k=cs/s; //求面积比例 
x=(x+k*x0)/(1+k); 
y=(y+k*y0)/(1+k); 
s += cs; 
} 
tp.x=x; 
tp.y=y; 
return tp; 
} 

/* 给定一简单多边形，找出一个肯定在该多边形内的点 
定理1：每个多边形至少有一个凸顶点 
定理2：顶点数>=4的简单多边形至少有一条对角线 
结论： x坐标最大，最小的点肯定是凸顶点 
         y坐标最大，最小的点肯定是凸顶点            
*/ 
POINT a_point_insidepoly(int vcount,POINT polygon[]) 
{ 
     POINT v,a,b,r; 
     int i,index; 
     v=polygon[0]; 
     index=0; 
     for(i=1;i<vcount;i++) //寻找一个凸顶点 
{ 
if(polygon[i].y<v.y) 
{ 
v=polygon[i]; 
index=i; 
} 
} 
a=polygon[(index-1+vcount)%vcount]; //得到v的前一个顶点 
b=polygon[(index+1)%vcount]; //得到v的后一个顶点 
POINT tri[3],q; 
tri[0]=a;tri[1]=v;tri[2]=b; 
double md=INF; 
int in1=index; 
bool bin=false; 
for(i=0;i<vcount;i++) //寻找在三角形avb内且离顶点v最近的顶点q 
{ 
if(i == index)continue; 
if(i == (index-1+vcount)%vcount)continue; 
if(i == (index+1)%vcount)continue; 
if(!InsideConvexPolygon(3,tri,polygon[i]))continue; 
bin=true; 
if(dist(v,polygon[i])<md) 
{ 
q=polygon[i]; 
md=dist(v,q); 
} 
} 
if(!bin) //没有顶点在三角形avb内，返回线段ab中点 
{ 
r.x=(a.x+b.x)/2; 
r.y=(a.y+b.y)/2; 
return r; 
} 
r.x=(v.x+q.x)/2; //返回线段vq的中点 
r.y=(v.y+q.y)/2; 
return r; 
} 

/* 求从多边形外一点p出发到一个简单多边形的切线,如果存在返回切点,其中rp点是右切点,lp是左切点 
注意：p点一定要在多边形外 
输入顶点序列是逆时针排列 
原 理:如果点在多边形内肯定无切线;凸多边形有唯一的两个切点,凹多边形就可能有多于两个的切点; 
如果polygon是凸多边形，切点只有两个只要找到就可以,可以化简此算法 
如果是凹多边形还有一种算法可以求解:先求凹多边形的凸包,然后求凸包的切线 
*/ 
void pointtangentpoly(int vcount,POINT polygon[],POINT p,POINT &rp,POINT &lp) 
{ 
LINESEG ep,en; 
bool blp,bln; 
rp=polygon[0]; 
lp=polygon[0]; 
for(int i=1;i<vcount;i++) 
{ 
ep.s=polygon[(i+vcount-1)%vcount]; 
ep.e=polygon[i]; 
en.s=polygon[i]; 
en.e=polygon[(i+1)%vcount]; 
blp=multiply(ep.e,p,ep.s)>=0;                  // p is to the left of pre edge 
bln=multiply(en.e,p,en.s)>=0;                  // p is to the left of next edge 
          if(!blp&&bln) 
{ 
if(multiply(polygon[i],rp,p)>0)             // polygon[i] is above rp 
rp=polygon[i]; 
} 
if(blp&&!bln) 
{ 
if(multiply(lp,polygon[i],p)>0)             // polygon[i] is below lp 
lp=polygon[i]; 
} 
} 
return ; 
} 

// 如果多边形polygon的核存在，返回true，返回核上的一点p.顶点按逆时针方向输入  
bool core_exist(int vcount,POINT polygon[],POINT &p) 
{ 
int i,j,k; 
LINESEG l; 
LINE lineset[MAXV]; 
for(i=0;i<vcount;i++) 
{ 
lineset[i]=makeline(polygon[i],polygon[(i+1)%vcount]); 
} 
for(i=0;i<vcount;i++) 
{ 
for(j=0;j<vcount;j++) 
{ 
if(i == j)continue; 
if(lineintersect(lineset[i],lineset[j],p)) 
{ 
for(k=0;k<vcount;k++) 
{ 
l.s=polygon[k]; 
l.e=polygon[(k+1)%vcount]; 
if(multiply(p,l.e,l.s)>0)      
//多边形顶点按逆时针方向排列，核肯定在每条边的左侧或边上 
break; 
} 
if(k == vcount)               //找到了一个核上的点 
break; 
} 
} 
if(j<vcount) 
break; 
} 
if(i<vcount) 
return true; 
else 
return false; 
} 



/*************************\ 
* * 
* 圆的基本运算 * 
* * 
\*************************/ 


/* 返回值： 点p在圆内(包括边界)时，返回true 
用途： 因为圆为凸集，所以判断点集，折线，多边形是否在圆内时，只需要逐一判断点是否在圆内即可。 
*/ 
bool point_in_circle(POINT o,double r,POINT p) 
{ 
double d2=(p.x-o.x)*(p.x-o.x)+(p.y-o.y)*(p.y-o.y); 
double r2=r*r; 
return d2<r2||abs(d2-r2)<EP; 
} 


/* 用 途：求不共线的三点确定一个圆 
输 入：三个点p1,p2,p3 
返回值：如果三点共线，返回false；反之，返回true。圆心由q返回，半径由r返回 
*/ 
bool cocircle(POINT p1,POINT p2,POINT p3,POINT &q,double &r) 
{ 
double x12=p2.x-p1.x; 
double y12=p2.y-p1.y; 
double x13=p3.x-p1.x; 
double y13=p3.y-p1.y; 
double z2=x12*(p1.x+p2.x)+y12*(p1.y+p2.y); 
double z3=x13*(p1.x+p3.x)+y13*(p1.y+p3.y); 
double d=2.0*(x12*(p3.y-p2.y)-y12*(p3.x-p2.x)); 
if(abs(d)<EP) //共线，圆不存在 
return false; 
q.x=(y13*z2-y12*z3)/d; 
q.y=(x12*z3-x13*z2)/d; 
r=dist(p1,q); 
return true; 
} 
int line_circle(LINE l,POINT o,double r,POINT &p1,POINT &p2) 
{ 



return true; 
} 


/**************************\ 
* * 
* 矩形的基本运算 * 
* * 
\**************************/ 


/* 
说明：因为矩形的特殊性，常用算法可以化简： 
1.判断矩形是否包含点 
只要判断该点的横坐标和纵坐标是否夹在矩形的左右边和上下边之间。 
2.判断线段、折线、多边形是否在矩形中 
因为矩形是个凸集，所以只要判断所有端点是否都在矩形中就可以了。 
3.判断圆是否在矩形中 
圆在矩形中的充要条件是：圆心在矩形中且圆的半径小于等于圆心到矩形四边的距离的最小值。 
*/ 
// 已知矩形的三个顶点(a,b,c)，计算第四个顶点d的坐标. 注意：已知的三个顶点可以是无序的 
POINT rect4th(POINT a,POINT b,POINT c) 
{ 
POINT d; 
if(abs(dotmultiply(a,b,c))<EP) // 说明c点是直角拐角处 
{ 
d.x=a.x+b.x-c.x; 
d.y=a.y+b.y-c.y; 
} 
if(abs(dotmultiply(a,c,b))<EP) // 说明b点是直角拐角处 
{ 
d.x=a.x+c.x-b.x; 
d.y=a.y+c.y-b.x; 
} 
if(abs(dotmultiply(c,b,a))<EP) // 说明a点是直角拐角处 
{ 
d.x=c.x+b.x-a.x; 
d.y=c.y+b.y-a.y; 
} 
return d; 
} 



/*************************\ 
* * 
* 常用算法的描述 * 
* * 
\*************************/ 

/* 尚未实现的算法： 
1. 求包含点集的最小圆 
2. 求多边形的交 
3. 简单多边形的三角剖分 
4. 寻找包含点集的最小矩形 
5. 折线的化简 
6. 判断矩形是否在矩形中 
7. 判断矩形能否放在矩形中 
8. 矩形并的面积与周长 
9. 矩形并的轮廓 
10.矩形并的闭包 
11.矩形的交 
12.点集中的最近点对 
13.多边形的并 
14.圆的交与并 
15.直线与圆的关系 
16.线段与圆的关系 
17.求多边形的核监视摄象机 
18.求点集中不相交点对 railwai 
*/ 

/* 寻找包含点集的最小矩形 
原理：该矩形至少一条边与点集的凸壳的某条边共线 
First take the convex hull of the points. Let the resulting convex 
polygon be P. It has been known for some time that the minimum 
area rectangle enclosing P must have one rectangle side flush with 
(i.e., collinear with and overlapping) one edge of P. This geometric 
fact was used by Godfried Toussaint to develop the "rotating calipers" 
algorithm in a hard-to-find 1983 paper, "Solving Geometric Problems 
with the Rotating Calipers" (Proc. IEEE MELECON). The algorithm 
rotates a surrounding rectangle from one flush edge to the next, 
keeping track of the minimum area for each edge. It achieves O(n) 
time (after hull computation). See the "Rotating Calipers Homepage" 
http://www.cs.mcgill.ca/~orm/rotcal.frame.html for a description 
and applet. 
*/ 

/* 折线的化简 伪码如下： 
Input: tol = the approximation tolerance 
L = {V0,V1,...,Vn-1} is any n-vertex polyline 

Set start = 0; 
Set k = 0; 
Set W0 = V0; 
for each vertex Vi (i=1,n-1) 
{ 
if Vi is within tol from Vstart 
then ignore it, and continue with the next vertex 

Vi is further than tol away from Vstart 
so add it as a new vertex of the reduced polyline 
Increment k++; 
Set Wk = Vi; 
Set start = i; as the new initial vertex 
} 

Output: W = {W0,W1,...,Wk-1} = the k-vertex simplified polyline 
*/ 





/********************\ 
* * 
* 补充 * 
* * 
\********************/ 

两圆关系： 

/* 两圆： 
相离： return 1； 
外切： return 2； 
相交： return 3； 
内切： return 4； 
内含： return 5； 
*/ 
int CircleRelation(POINT p1, double r1, POINT p2, double r2) 
{ 
double d = sqrt( (p1.x-p2.x)*(p1.x-p2.x)+(p1.y-p2.y)*(p1.y-p2.y) ); 


if( fabs(d-r1-r2) < EP ) // 必须保证前两个if先被判定！ 
return 2; 
if( fabs(d-fabs(r1-r2)) < EP ) 
return 4; 
if( d > r1+r2 ) 
        return 1; 
     if( d < fabs(r1-r2) ) 
return 5; 
if( fabs(r1-r2) < d && d < r1+r2 ) 
return 3; 
return 0; // indicate an error! 
} 


判断圆是否在矩形内： 

// 判定圆是否在矩形内，是就返回true（设矩形水平，且其四个顶点由左上开始按顺时针排列） 
// 调用ptoldist函数，在第4页 
bool CircleRecRelation(POINT pc, double r, POINT pr1, POINT pr2, POINT pr3, POINT pr4) 
{ 
if( pr1.x < pc.x && pc.x < pr2.x && pr3.y < pc.y && pc.y < pr2.y ) 
{ 
LINESEG line1(pr1, pr2); 
LINESEG line2(pr2, pr3); 
LINESEG line3(pr3, pr4); 
LINESEG line4(pr4, pr1); 
if( r<ptoldist(pc,line1) && r<ptoldist(pc,line2) && 
r<ptoldist(pc,line3) && r<ptoldist(pc,line4) ) 
return true; 
} 
return false; 
} 


点到平面的距离： 

// 点到平面的距离,平面用一般式表示ax+by+cz+d=0 
double P2planeDist(double x, double y, double z, double a, double b, double c, double d) 
{ 
return fabs(a*x+b*y+c*z+d) / sqrt(a*a+b*b+c*c); 
} 


点是否在直线同侧： 

// 两个点是否在直线同侧，是则返回true 
bool SameSide(POINT p1, POINT p2, LINE line) 
{ 
return (line.a * p1.x + line.b * p1.y + line.c) * 
(line.a * p2.x + line.b * p2.y + line.c) > 0; 
} 


镜面反射线： 

// 已知入射线、镜面，求反射线。 
// a1,b1,c1为镜面直线方程(a1 x + b1 y + c1 = 0 ,下同)系数;  
a2,b2,c2为入射光直线方程系数;  
a,b,c为反射光直线方程系数. 
// 光是有方向的，使用时注意：入射光向量:<-b2,a2>；反射光向量:<b,-a>. 
// 不要忘记结果中可能会有"negative zeros" 

void reflect(double a1,double b1,double c1,double a2,double b2,double c2,double &a,double &b,double &c) 
{ 
     double n,m; 
     double tpb,tpa; 
     tpb=b1*b2+a1*a2; 
     tpa=a2*b1-a1*b2; 
     m=(tpb*b1+tpa*a1)/(b1*b1+a1*a1); 
     n=(tpa*b1-tpb*a1)/(b1*b1+a1*a1); 
     if(fabs(a1*b2-a2*b1)<1e-20) 
{ 
a=a2;b=b2;c=c2; 
return; 
} 
double xx,yy; //(xx,yy)是入射线与镜面的交点。 
xx=(b1*c2-b2*c1)/(a1*b2-a2*b1); 
yy=(a2*c1-a1*c2)/(a1*b2-a2*b1); 
a=n; 
b=-m; 
c=m*yy-xx*n; 
} 


矩形包含： 

bool r2inr1(double A,double B,double C,double D) // 矩形2（C，D）是否在1（A，B）内 
{ 
double X,Y,L,K,DMax; 
if (A < B) 
{ 
double tmp = A; 
A = B; 
B = tmp; 
} 
if (C < D) 
{ 
double tmp = C; 
C = D; 
D = tmp; 
} 

if (A > C && B > D)                   // trivial case  
        return true; 
      else 
        if (D >= B) 
          return false; 
        else 
        { 
          X = sqrt(A * A + B * B);           // outer rectangle's diagonal  
          Y = sqrt(C * C + D * D);           // inner rectangle's diagonal  
          if (Y < B) // check for marginal conditions 
return true; // the inner rectangle can freely rotate inside 
else 
if (Y > X) 
              return false; 
            else 
            { 
              L = (B - sqrt(Y * Y - A * A)) / 2; 
              K = (A - sqrt(Y * Y - B * B)) / 2; 
              DMax = sqrt(L * L + K * K); 
              if (D >= DMax) 
                return false; 
              else 
                return true; 
            } 
        } 
}

两圆交点： 

// 两圆已经相交（相切） 
void    c2point(POINT p1,double r1,POINT p2,double r2,POINT &rp1,POINT &rp2) 
{ 
double a,b,r; 
a=p2.x-p1.x; 
b=p2.y-p1.y; 
r=(a*a+b*b+r1*r1-r2*r2)/2; 
if(a==0&&b!=0) 
{ 
rp1.y=rp2.y=r/b; 
rp1.x=sqrt(r1*r1-rp1.y*rp1.y); 
rp2.x=-rp1.x; 
} 
else if(a!=0&&b==0) 
{ 
rp1.x=rp2.x=r/a; 
rp1.y=sqrt(r1*r1-rp1.x*rp2.x); 
rp2.y=-rp1.y; 
} 
else if(a!=0&&b!=0) 
{ 
double delta; 
delta=b*b*r*r-(a*a+b*b)*(r*r-r1*r1*a*a); 
rp1.y=(b*r+sqrt(delta))/(a*a+b*b); 
rp2.y=(b*r-sqrt(delta))/(a*a+b*b); 
rp1.x=(r-b*rp1.y)/a; 
rp2.x=(r-b*rp2.y)/a; 
} 

rp1.x+=p1.x; 
rp1.y+=p1.y; 
rp2.x+=p1.x; 
rp2.y+=p1.y; 
} 


两圆公共面积： 

// 必须保证相交 
double c2area(POINT p1,double r1,POINT p2,double r2) 
{ 
POINT rp1,rp2; 
c2point(p1,r1,p2,r2,rp1,rp2); 

if(r1>r2) //保证r2>r1 
{ 
swap(p1,p2); 
swap(r1,r2); 
} 
double a,b,rr; 
a=p1.x-p2.x; 
b=p1.y-p2.y; 
rr=sqrt(a*a+b*b); 

double dx1,dy1,dx2,dy2; 
double sita1,sita2; 
dx1=rp1.x-p1.x; 
dy1=rp1.y-p1.y; 
dx2=rp2.x-p1.x; 
dy2=rp2.y-p1.y; 
sita1=acos((dx1*dx2+dy1*dy2)/r1/r1); 

dx1=rp1.x-p2.x; 
dy1=rp1.y-p2.y; 
dx2=rp2.x-p2.x; 
dy2=rp2.y-p2.y; 
sita2=acos((dx1*dx2+dy1*dy2)/r2/r2); 
double s=0; 
if(rr<r2) //相交弧为优弧 
s=r1*r1*(PI-sita1/2+sin(sita1)/2)+r2*r2*(sita2-sin(sita2))/2; 
else //相交弧为劣弧 
s=(r1*r1*(sita1-sin(sita1))+r2*r2*(sita2-sin(sita2)))/2; 

return s; 
} 

圆和直线关系： 

//0----相离 1----相切 2----相交 
int clpoint(POINT p,double r,double a,double b,double c,POINT &rp1,POINT &rp2) 
{ 
int res=0; 

c=c+a*p.x+b*p.y; 
double tmp; 
if(a==0&&b!=0) 
{ 
tmp=-c/b; 
if(r*r<tmp*tmp) 
res=0; 
else if(r*r==tmp*tmp) 
{ 
res=1; 
rp1.y=tmp; 
rp1.x=0; 
} 
else 
{ 
res=2; 
rp1.y=rp2.y=tmp; 
rp1.x=sqrt(r*r-tmp*tmp); 
rp2.x=-rp1.x; 
} 
} 
else if(a!=0&&b==0) 
{ 
tmp=-c/a; 
if(r*r<tmp*tmp) 
res=0; 
else if(r*r==tmp*tmp) 
{ 
res=1; 
rp1.x=tmp; 
rp1.y=0; 
} 
else 
{ 
res=2; 
rp1.x=rp2.x=tmp; 
rp1.y=sqrt(r*r-tmp*tmp); 
rp2.y=-rp1.y; 
} 
} 
else if(a!=0&&b!=0) 
{ 
double delta; 
delta=b*b*c*c-(a*a+b*b)*(c*c-a*a*r*r); 
if(delta<0) 
res=0; 
else if(delta==0) 
{ 
res=1; 
rp1.y=-b*c/(a*a+b*b); 
rp1.x=(-c-b*rp1.y)/a; 
} 
else 
{ 
res=2; 
rp1.y=(-b*c+sqrt(delta))/(a*a+b*b); 
rp2.y=(-b*c-sqrt(delta))/(a*a+b*b); 
rp1.x=(-c-b*rp1.y)/a; 
rp2.x=(-c-b*rp2.y)/a; 
} 
} 
rp1.x+=p.x; 
rp1.y+=p.y; 
rp2.x+=p.x; 
rp2.y+=p.y; 
return res; 
} 


内切圆： 

void incircle(POINT p1,POINT p2,POINT p3,POINT &rp,double &r) 
{ 
double dx31,dy31,dx21,dy21,d31,d21,a1,b1,c1; 
dx31=p3.x-p1.x; 
dy31=p3.y-p1.y; 
dx21=p2.x-p1.x; 
dy21=p2.y-p1.y; 

d31=sqrt(dx31*dx31+dy31*dy31); 
d21=sqrt(dx21*dx21+dy21*dy21); 
a1=dx31*d21-dx21*d31; 
b1=dy31*d21-dy21*d31; 
c1=a1*p1.x+b1*p1.y; 

double dx32,dy32,dx12,dy12,d32,d12,a2,b2,c2; 
dx32=p3.x-p2.x; 
dy32=p3.y-p2.y; 
dx12=-dx21; 
dy12=-dy21; 

d32=sqrt(dx32*dx32+dy32*dy32); 
d12=d21; 
a2=dx12*d32-dx32*d12; 
b2=dy12*d32-dy32*d12; 
c2=a2*p2.x+b2*p2.y; 

rp.x=(c1*b2-c2*b1)/(a1*b2-a2*b1); 
rp.y=(c2*a1-c1*a2)/(a1*b2-a2*b1); 
r=fabs(dy21*rp.x-dx21*rp.y+dx21*p1.y-dy21*p1.x)/d21; 
} 


求切点： 

// p---圆心坐标， r---圆半径， sp---圆外一点， rp1,rp2---切点坐标 
void cutpoint(POINT p,double r,POINT sp,POINT &rp1,POINT &rp2) 
{ 
POINT p2; 
p2.x=(p.x+sp.x)/2; 
p2.y=(p.y+sp.y)/2; 

double dx2,dy2,r2; 
dx2=p2.x-p.x; 
dy2=p2.y-p.y; 
r2=sqrt(dx2*dx2+dy2*dy2); 
c2point(p,r,p2,r2,rp1,rp2); 
} 


线段的左右旋： 

/* l2在l1的左/右方向（l1为基准线） 
返回 0： 重合； 
返回 1： 右旋； 
返回 –1： 左旋； 
*/ 
int rotat(LINESEG l1,LINESEG l2) 
{ 
double dx1,dx2,dy1,dy2; 
dx1=l1.s.x-l1.e.x; 
dy1=l1.s.y-l1.e.y; 
dx2=l2.s.x-l2.e.x; 
dy2=l2.s.y-l2.e.y; 

double d; 
d=dx1*dy2-dx2*dy1; 
if(d==0) 
return 0; 
else if(d>0) 
        return -1; 
     else 
        return 1; 
} 


公式： 

球坐标公式： 
直角坐标为 P(x, y, z) 时，对应的球坐标是(rsinφcosθ, rsinφsinθ, rcosφ),其中φ是向量OP与Z轴的夹角，范围[0，π]；是OP在XOY面上的投影到X轴的旋角，范围[0，2π]  

直线的一般方程转化成向量方程： 
ax+by+c=0 
x-x0       y-y0 
     ------ = ------- // (x0,y0)为直线上一点，m,n为向量 
m          n 
转换关系： 
a=n；b=-m；c=m·y0-n·x0； 
m=-b; n=a; 

三点平面方程： 
三点为P1，P2，P3 
设向量    M1=P2-P1; M2=P3-P1; 
平面法向量：    M=M1 x M2 （） 
平面方程：      M.i(x-P1.x)+M.j(y-P1.y)+M.k(z-P1.z)=0