【總結】計算幾何模板

計算幾何模板

1. 幾何公式

1.1 角：

1.正弦定理：

a/sinA = b/sinB = c/sinC = 2R

2.餘弦定理：

a² = b²+c²-2bc*cosA

b² = a²+c²-2ac*cosB

c² = a²+b²-2ab*cosC

3.倍角公式:

sin(2α)=2sinαcosα

cos(2α)=(cosα)²-1=1-2(sinα)²

tan(2α)=2tanα/[1-(tanα)²]

sin(3α)=3sinα-4(sinα)^3

cos(3α)=4(cosα)^3-3cosα

1.2 三角形：

在△ABC中,設AB=c,AC=b,CB=a, r為內切圓半徑, R為外接圓半徑

1.半周長 P = a+b+c;

2.面積 S =a*b*sin(C)/2;

S = sqrt(P(P-a)(P-b)(P-c))

3.a邊中線長 Ma =(1/2)*sqrt(2b²+2c²-a²)

Ma = (1/2)*sqrt(b²+c²+2bc*cosA)

4.a邊高長 ha =c*sinB=b*sinC

ha = sqrt(b²-(a²+b²-c²)²/(2a)²)

5.A角平分線長 la =2bc*cos(A/2)/(b+c)

la = sqrt{bc[(b+c)²-a²]}/(b+c)

6.內切圓半徑 r = S/P;

r = 4R*sin(A/2)sin(B/2)sin(C/2);

r = P*tan(A/2)tan(B/2)tan(C/2)

7.外接圓半徑 R =a/(2sinA)

R = abc/(4P)

R = abc/[2r(a+b+c)]

1.3 四邊形：

D1,D2為對角線,M對角線中點連線,A為對角線夾角,P為四邊形半周長

1. a^2+b^2+c^2+d^2 = D1^2+D2^2+4M^2

2. 面積 S=D1D2sin(A)/2

(以下對於圓的內接四邊形)

3. ac+bd=D1D2

4. S=sqrt((P-a)(P-b)(P-c)(P-d))

1.4 正n邊形：

R為外接圓半徑,r為內切圓半徑

1. 中心角 A=2PI/n

2. 內角 C=(n-2)PI/n

3. 邊長a=2sqrt(R^2-r^2)=2Rsin(A/2)=2rtan(A/2)

4. 面積S=nar/2=nr^2tan(A/2)=nR^2sin(A)/2=na^2/(4tan(A/2))

1.5 圓：

A為所求的角度

1. 弧長 l=rA

2. 弦長 a=2sqrt(2hr-h^2)=2rsin(A/2)

3. 弓形高h=r-sqrt(r^2-a^2/4)=r(1-cos(A/2))=atan(A/4)/2

4. 扇形面積S1=rl/2=r^2A/2

5. 弓形面積S2=(rl-a(r-h))/2=r^2(A-sin(A))/2

1.6 稜柱：

1. 體積 V=Ah,A為底面積,h為高

2. 側面積 S=lp,l為稜長,p為直截面周長

3. 全面積 T=S+2A

1.7 稜錐：

1. 體積 V=Ah/3,A為底面積,h為高

(以下對正稜錐)

2. 側面積 S=lp/2,l為斜高,p為底面周長

3. 全面積 T=S+A

1.8 稜臺：

1. 體積V=(A1+A2+sqrt(A1A2))h/3,A1.A2為上下底面積,h為高

(以下為正稜臺)

2. 側面積S=(p1+p2)l/2,p1.p2為上下底面周長,l為斜高

3. 全面積 T=S+A1+A2

1.9 圓柱：

1. 側面積 S=2PIrh

2. 全面積 T=2PIr(h+r)

3. 體積 V=PIr^2h

1.10圓錐：

1. 母線l=sqrt(h^2+r^2)

2. 側面積 S=PIrl

3. 全面積 T=PIr(l+r)

4. 體積 V=PIr^2h/3

1.11 圓臺：

1. 母線l=sqrt(h^2+(r1-r2)^2)

2. 側面積 S=PI(r1+r2)l

3. 全面積T=PIr1(l+r1)+PIr2(l+r2)

4. 體積V=PI(r1^2+r2^2+r1r2)h/3

1.12 球：

1. 全面積 T=4PIr^2

2. 體積 V=4PIr^3/3

1.13 球檯：

1. 側面積 S=2PIrh

2. 全面積T=PI(2rh+r1^2+r2^2)

3. 體積V=PIh(3(r1^2+r2^2)+h^2)/6

1.14 球扇形：

1. 全面積T=PIr(2h+r0),h為球冠高,r0為球冠底面半徑

2. 體積 V=2PIr^2h/3

2. 幾何向量表示

2.1：定義

struct Point{ double x,y; }; //點座標

struct V{ Point start,endd; }; //兩點式

struct V2{ double A,B,C; }; //一般式 ax+by+c = 0

2.2：轉化

//兩點式化一般式

V2 V2toV1(V Line) // y=kx+c k=y/x

{

V2temp;

temp.A = Line.endd.y - Line.start.y;

temp.B = Line.start.x - Line.endd.x;

temp.C = Line.endd.x * Line.start.y - Line.start.x * Line.endd.y;

return temp;

}

//兩點式起點化為 0,0

V change(V a)

{

Vb;

b.start.x= b.start.y = 0;

b.endd.x = a.endd.x - a.start.x;

b.endd.y = a.endd.y - a.start.y;

return b;

}

向量運算需要先把兩點式起點化為 0,0

2.3 向量相加：

V add(V a,V b)

{

Vaa = change(a);

Vbb = change(b);

Vc;

c.start.x = c.start.y = 0;

c.endd.x = aa.endd.x+bb.endd.x;

c.endd.y = aa.endd.y+bb.endd.y;

return c;

}

2.4 向量相減：

V dec(V a,V b)

{

Vaa = change(a);

Vbb = change(b);

Vc;

c.start.x = 0;

c.start.y = 0;

c.endd.x = aa.endd.x-bb.endd.x;

c.endd.y = aa.endd.y-bb.endd.y;

return c;

}

2.5 向量點積：

double dmult(V a,V b)

{

double result;

Vaa = change(a);

Vbb = change(b);

result = aa.endd.x*bb.endd.x + aa.endd.y*bb.endd.y;

return result;

}

double dmult(Point a,Point b,Point c) //(ca)•(cb)

{

return (a.x-c.x)*(b.x-c.x)+(a.y-c.y)*(b.y-c.y);

}

可以通過向量點積計算向量夾角

a•b = |a|*|b|*cos(∠)

a•b = x1*x2 + y1*y2;

2.6 向量叉積：

double xmult(V a,V b)

{

double result;

Vaa = change(a);

Vbb = change(b);

result = aa.endd.x*bb.endd.y - aa.endd.y*bb.endd.x;

return result;

}

double xmult(Point a,Point b,Point c) //(ca)×(cb)

{

return (a.x-c.x)*(b.y-c.y)-(b.x-c.x)*(a.y-c.y);

}

叉積可以計算三角形面積(取絕對值)

可以通過向量叉積的符號判斷向量的順逆時針關係

a×b > 0 , a 在 b的順時針方向

a×b < 0 , a 在 b的逆時針方向

a×b = 0 , a 在 b同方向或反方向

a×b = x1*y2 - x2*y1;

3. 點線關係

3.1 點與點問題：

//兩點距離

double dist(Point a,Point b) {

return sqrt((a.x-b.x)*(a.x-b.x)+(a.y-b.y)*(a.y-b.y));

}

//判斷三點共線

bool dot_line(Point a,Point b,Point c)

{

double temp = xmult(b,c,a);

if(fabs(temp) < eps) //滿足向量 AC×AB == 0

return true;

return false;

}

3.2 點與線問題：

//判斷點C是否線上段A，B上(包括端點)

bool dot_seg(Point a,Point b,Point c)

{

//滿足三點共線，並且滿足C在AB構成的矩形內

if(dot_line(a,b,c) && c.x >= min(a.x,b.x) && c.x<= max(a.x,b.x) && c.y >= min(a.y,b.y) && c.y <=max(a.y,b.y))

return true;

return false;

}

//判斷兩點線上的同側，點線上上或異側返回false

bool same_side(Point a,Point b,V line) {

return xmult(line.start,a,line.endd)*xmult(line.start,b,line.endd) >eps;

}

//判斷兩點線上的異側，點線上上或同側返回false

bool opp_side(Point a,Point b,V line) {

return xmult(line.start,a,line.endd)*xmult(line.start,b,line.endd) <-eps;

}

//獲得某點對於線垂線的另個點

Point get_cLine(Point a, V Line)

{

Point t = a;

t.x += Line.start.y-Line.endd.y;

t.y += Line.endd.x-Line.start.x;

return t;

}

//點在直線上的最近點（垂點）

Point ptoline(Point a,V Line)

{

Point t = get_cLine(a, Line); //得到垂線的另個點

VLine2; //垂線

Line2.start = a;

Line2.endd = t;

return intersection(Line2,Line); //兩直線相交的交點

}

//點到直線的距離（點到垂點的距離）

double dist_ptoline(Point a,V Line) {

returnfabs(xmult(a,Line.start,Line.endd)/dist(Line.start,Line.endd)); //叉積求面積/底長

}

//點到線段上的最近點

Point ptoseg(Point a,V Line)

{

Point t = get_cLine(a, Line); //得到垂線的另個點

VLine2;

Line2.start = a;

Line2.endd = t;

if(!seg_Line_intersect(Line,Line2)) //如果直線Line2與線段Line不相交

return dist(a,Line.start)<dist(a,Line.endd)?Line.start:Line.endd;

return ptoline(a,Line); //點在直線上的垂點

}

//點到線段距離

double dist_ptoseg(Point a,V Line) {

Point t = get_cLine(a, Line); //得到垂線的另個點

VLine2;

Line2.start = a;

Line2.endd = t;

if(!seg_Line_intersect(Line,Line2)) //如果直線Line2與線段Line不相交

returndist(a,Line.start)<dist(a,Line.endd)?dist(a,Line.start):dist(a,Line.endd);

return dist_ptoline(a,Line); //點到直線的距離（點到垂點的距離）

}

//向量V以P為頂點逆時針旋轉angle並放大scale倍

Point rotat(Point v,Point p,doubleangle,double scale)

{

Point ret=p;

v.x-=p.x,v.y-=p.y;

p.x=scale*cos(angle);

p.y=scale*sin(angle);

ret.x+=v.x*p.x-v.y*p.y;

ret.y+=v.x*p.y+v.y*p.x;

return ret;

}

//p點關於直線L的對稱點(用一般式,未驗證)

Point symmetricalPointofLine(Point p, V2Line)

{

Point p2;

double d;

d= Line.a * Line.a + Line.b * Line.b;

p2.x = (Line.b * Line.b * p.x - Line.a * Line.a * p.x -

2 * Line.a * Line.b * p.y - 2 * Line.a * Line.c) / d;

p2.y = (Line.a * Line.a * p.y - Line.b * Line.b * p.y -

2 * Line.a * Line.b * p.x - 2 * Line.b * Line.c) / d;

return p2;

}

3.3 線與線問題：

//兩直線平行

bool parallel(V Line1,V Line2) {

return fabs((Line1.start.x-Line1.endd.x)*(Line2.start.y-Line2.endd.y) -(Line2.start.x-Line2.endd.x)*(Line1.start.y-Line1.endd.y)) < eps;

}

//兩直線垂直

bool perpendicular(V Line1,V Line2) {

return fabs((Line1.start.x-Line1.endd.x)*(Line2.start.x-Line2.endd.x) +(Line1.start.y-Line1.endd.y)*(Line2.start.y-Line2.endd.y)) < eps;

}

//判斷直線Line2與線段Line1相交,包括端點,不包括平行

bool seg_Line_intersect(V Line1,V Line2)

{

return !same_side(Line1.start,Line1.endd,Line2);

}

//判斷兩線段相交，包括端點和部分重合

bool intersect_in(V Line1,V Line2)

{

if(!dot_line(Line1.start,Line1.endd,Line2.start) ||!dot_line(Line1.start,Line1.endd,Line2.endd))

return !same_side(Line1.start,Line1.endd,Line2) &&!same_side(Line2.start,Line2.endd,Line1);

returndot_seg(Line1.start,Line1.endd,Line2.start)||dot_seg(Line1.start,Line1.endd,Line2.endd)||dot_seg(Line2.start,Line2.endd,Line1.start)||dot_seg(Line2.start,Line2.endd,Line1.endd);

}

//判兩線段相交,不包括端點和部分重合

bool intersect_ex(V Line1,V Line2)

{

return opp_side(Line1.start,Line1.endd,Line2) &&opp_side(Line2.start,Line2.endd,Line1);

}

//計算兩直線交點,注意事先判斷直線是否平行!

//線段交點請另外判線段相交(同時還是要判斷是否平行!)

Point intersection(V u,V v)

{

Point ret=u.start;

doublet=((u.start.x-v.start.x)*(v.start.y-v.endd.y)-(u.start.y-v.start.y)*(v.start.x-v.endd.x))

/((u.start.x-u.endd.x)*(v.start.y-v.endd.y)-(u.start.y-u.endd.y)*(v.start.x-v.endd.x));

ret.x+=(u.endd.x-u.start.x)*t;

ret.y+=(u.endd.y-u.start.y)*t;

return ret;

}

4三角形

4.1 面積：

//計算三角形面積,輸入三頂點

double area_triangle(Point p1,Pointp2,Point p3) {

return fabs(xmult(p1,p2,p3))/2;

}

//計算三角形面積,輸入三邊長

double area_triangle(double a,doubleb,double c) {

double s=(a+b+c)/2;

return sqrt(s*(s-a)*(s-b)*(s-c));

}

4.2 三角形與點：

//判斷點p是否在三角形內或上

//點在平面內與三角形三個點構成的三個三角形的面積和與原三角形的面積比較,如果相等，即可判斷在三角形內或上（若需要區分，先判斷點與三條邊）

bool dot_triangle(Point p,Point a,Pointb,Point c)

{

double abx = fabs(xmult(a,b,p));

double acx = fabs(xmult(a,c,p));

double bcx = fabs(xmult(b,c,p));

double abc = fabs(xmult(a,b,c));

if(fabs(abx+acx+bcx - abc) < eps)

return true;

return false;

}

4.3 五個重要點：

//外心三角形三邊垂直平分線的交點

Point circumcenter(Point a,Point b,Point c){

Vu,v; //垂直平分線

u.start.x=(a.x+b.x)/2;

u.start.y=(a.y+b.y)/2;

u.endd.x=u.start.x-a.y+b.y;

u.endd.y=u.start.y+a.x-b.x;

v.start.x=(a.x+c.x)/2;

v.start.y=(a.y+c.y)/2;

v.endd.x=v.start.x-a.y+c.y;

v.endd.y=v.start.y+a.x-c.x;

return intersection(u,v);

}

//內心三個內角的三條角平分線的交點

Point incenter(Point a,Point b,Point c) {

Vu,v; //角平分線

double n,m; //角度

u.start = a;

m=atan2(b.y-a.y,b.x-a.x);

n=atan2(c.y-a.y,c.x-a.x);

u.endd.x=u.start.x+cos((m+n)/2);

u.endd.y=u.start.y+sin((m+n)/2);

v.start=b;

m=atan2(a.y-b.y,a.x-b.x);

n=atan2(c.y-b.y,c.x-b.x);

v.endd.x=v.start.x+cos((m+n)/2);

v.endd.y=v.start.y+sin((m+n)/2);

return intersection(u,v);

}

//垂心三角形的三條高或其延長線相交於一點

Point perpencenter(Point a,Point b,Point c){

Vu,v; //高

u.start=c;

u.endd.x=u.start.x-a.y+b.y;

u.endd.y=u.start.y+a.x-b.x;

v.start=b;

v.endd.x=v.start.x-a.y+c.y;

v.endd.y=v.start.y+a.x-c.x;

return intersection(u,v);

}

//重心是三角形三邊中線的交點

//到三角形三頂點距離的平方和最小的點

//三角形內到三邊距離之積最大的點

Point barycenter(Point a,Point b,Point c)

{

Vu,v; //中線

u.start.x=(a.x+b.x)/2;

u.start.y=(a.y+b.y)/2;

u.endd=c;

v.start.x=(a.x+c.x)/2;

v.start.y=(a.y+c.y)/2;

v.endd=b;

return intersection(u,v);

}

//費馬點

//到三角形三頂點距離之和最小的點

Point fermentpoint(Point a,Point b,Point c)

{

Point u,v;

double step=fabs(a.x)+fabs(a.y)+fabs(b.x)+fabs(b.y)+fabs(c.x)+fabs(c.y);

int i,j,k;

u.x=(a.x+b.x+c.x)/3;

u.y=(a.y+b.y+c.y)/3;

while (step>1e-10)

for (k=0; k<10; step/=2,k++)

for (i=-1; i<=1; i++)

for (j=-1; j<=1; j++)

{

v.x=u.x+step*i;

v.y=u.y+step*j;

if(dist(u,a)+dist(u,b)+dist(u,c)>dist(v,a)+dist(v,b)+dist(v,c))

u=v;

}

return u;

}

5多邊形

5.1 簡單多邊形面積：

// p[]為頂點集並是順序儲存的（不是的話，無法計算吧），n為多邊形頂點數

//儲存是從0開始 (未驗證)

double area_poly(Point* p,int n)

{

double ans = 0;

for(int i = 1; i < n-1; i++)

ans += xmult(p[i],p[i+1],p[0]);

return fabs(ans)/2;

}

5.2 判定多邊形：

int _sign(double temp)

{

if(temp < -eps) return 2;

else if(temp > eps) return 1;

else return 0;

}

//判定凸多邊形,頂點按順時針或逆時針給出,允許相鄰邊共線

int is_convex(Point* p,int n)

{

int i,s[3]= {1,1,1};

for (i=0; i<n&&s[1]|s[2]; i++)

s[_sign(xmult(p[(i+1)%n],p[(i+2)%n],p[i]))]=0;

return s[1]|s[2];

}

//判定凸多邊形,頂點按順時針或逆時針給出,不允許相鄰邊共線

int is_convex2(Point* p,int n)

{

int i,s[3]= {1,1,1};

for (i=0; i<n&&s[1]|s[2]; i++)

s[_sign(xmult(p[(i+1)%n],p[(i+2)%n],p[i]))]=0;

return s[0]&&s[1]|s[2];

}

5.3 點與多邊形：

//判點在凸多邊形內或多邊形邊上,頂點按順時針或逆時針給出

int inside_convex(Point q,int n,Point* p)

{

int i,s[3]= {1,1,1};

for (i=0; i<n&&s[1]|s[2]; i++)

s[_sign(xmult(p[(i+1)%n],q,p[i]))]=0;

return s[1]|s[2];

}

//判點在凸多邊形內,頂點按順時針或逆時針給出,在多邊形邊上返回0

int inside_convex_v2(Point q,int n,Point*p)

{

int i,s[3]= {1,1,1};

for (i=0; i<n&&s[0]&&s[1]|s[2]; i++)

s[_sign(xmult(p[(i+1)%n],q,p[i]))]=0;

return s[0]&&s[1]|s[2];

}

//判點在任意多邊形內,頂點按順時針或逆時針給出

//on_edge表示點在多邊形邊上時的返回值,offset為多邊形座標上限 (不是很懂這個)

int inside_polygon(Point q,int n,Point* p,inton_edge=1)

{

Point q2;

int i=0,coun;

while (i<n)

for (coun=i=0,q2.x=rand()/*+offset*/,q2.y=rand()/*+offset*/; i<n;i++)

if(fabs(xmult(q,p[i],p[(i+1)%n]))<eps&&(p[i].x-q.x)*(p[(i+1)%n].x-q.x)<eps&&(p[i].y-q.y)*(p[(i+1)%n].y-q.y)<eps)

return on_edge;

else if (fabs(xmult(q,q2,p[i])) < eps)

break;

else if(xmult(q,p[i],q2)*xmult(q,p[(i+1)%n],q2)<-eps&&xmult(p[i],q,p[(i+1)%n])*xmult(p[i],q2,p[(i+1)%n])<-eps)

coun++;

return coun&1;

}

5.4 多邊形重心：

//多邊形重心

Point barycenter(Point* p,int n) {

Point ret,t;

double t1=0,t2;

int i;

ret.x=ret.y=0;

for (i=1;i<n-1;i++)

if (fabs(t2=xmult(p[0],p[i],p[i+1]))>eps)

{

t=barycenter(p[0],p[i],p[i+1]); //三角形重心

ret.x+=t.x*t2;

ret.y+=t.y*t2;

t1+=t2;

}

if (fabs(t1)>eps)

ret.x/=t1,ret.y/=t1;

return ret;

}

5.5 多邊形與網格：

struct PPoint{ int x,y; };

//多邊形上的網格點個數(PPoint中 x,y 是int行)

int grid_onedge(int n,PPoint* p) {

int i,ret=0;

for (i=0; i<n; i++)

ret+=__gcd(abs(p[i].x-p[(i+1)%n].x),abs(p[i].y-p[(i+1)%n].y));

return ret;

}

//多邊形內的網格點個數

int grid_inside(int n,PPoint* p) {

int i,ret=0;

for (i=0; i<n; i++)

ret+=p[(i+1)%n].y*(p[i].x-p[(i+2)%n].x);

return (abs(ret)-grid_onedge(n,p))/2+1;

}

6圓

//判直線和圓相交,包括相切

int intersect_line_circle(Point c,doubler,Point p1,Point p2)

{

VLine;

Line.start = p1;

Line.endd = p2;

return dist_ptoline(c,Line)<r+eps; //點到直線的最近距離

}

//判線段和圓相交,包括端點和相切

int intersect_seg_circle(Point c,doubler,Point p1,Point p2)

{

double t1=dist(c,p1)-r,t2=dist(c,p2)-r;

Point t=c;

if (t1<eps||t2<eps)

return t1>-eps||t2>-eps;

t.x+=p1.y-p2.y;

t.y+=p2.x-p1.x;

VLine;

Line.start = p1;

Line.endd = p2;

returnxmult(p1,c,t)*xmult(p2,c,t)<eps&&dist_ptoline(c,Line)-r<eps;

}

//判圓和圓相交,包括相切

int intersect_circle_circle(Point c1,doubler1,Point c2,double r2) {

return dist(c1,c2)<r1+r2+eps&&dist(c1,c2)>fabs(r1-r2)-eps;

}

//計算圓上到點p最近點,如p與圓心重合,返回p本身

Point dot_to_circle(Point c,double r,Pointp) {

Point u,v;

if (dist(p,c)<eps)

return p;

u.x=c.x+r*fabs(c.x-p.x)/dist(c,p);

u.y=c.y+r*fabs(c.y-p.y)/dist(c,p)*((c.x-p.x)*(c.y-p.y)<0?-1:1);

v.x=c.x-r*fabs(c.x-p.x)/dist(c,p);

v.y=c.y-r*fabs(c.y-p.y)/dist(c,p)*((c.x-p.x)*(c.y-p.y)<0?-1:1);

return dist(u,p)<dist(v,p)?u:v;

}

//計算直線與圓的交點,保證直線與圓有交點

//計算線段與圓的交點可用這個函式後判點是否線上段上

void intersection_line_circle(Pointc,double r,V Line,Point& p1,Point& p2) {

Point p=c;

double t;

p.x+=Line.start.y-Line.endd.y;

p.y+=Line.endd.x-Line.start.x;

VLine2;

Line2.start = p;

Line2.endd = c;

p=intersection(Line2,Line);

t=sqrt(r*r-dist(p,c)*dist(p,c))/dist(Line.start,Line.endd);

p1.x=p.x+(Line.endd.x-Line.start.x)*t;

p1.y=p.y+(Line.endd.y-Line.start.y)*t;

p2.x=p.x-(Line.endd.x-Line.start.x)*t;

p2.y=p.y-(Line.endd.y-Line.start.y)*t;

}

//計算圓與圓的交點,保證圓與圓有交點,圓心不重合

void intersection_circle_circle(Pointc1,double r1,Point c2,double r2,Point& p1,Point& p2) {

Point u,v;

double t;

t=(1+(r1*r1-r2*r2)/dist(c1,c2)/dist(c1,c2))/2;

u.x=c1.x+(c2.x-c1.x)*t;

u.y=c1.y+(c2.y-c1.y)*t;

v.x=u.x+c1.y-c2.y;

v.y=u.y-c1.x+c2.x;

VLine;

Line.start = u;

Line.endd = v;

intersection_line_circle(c1,r1,Line,p1,p2);

}

//將向量p逆時針旋轉angle角度

Point Rotate(Point p,double angle) {

Point res;

res.x=p.x*cos(angle)-p.y*sin(angle);

res.y=p.x*sin(angle)+p.y*cos(angle);

return res;

}

//求圓外一點對圓(o,r)的兩個切點result1和result2

void TangentPoint_PC(Point poi,Pointo,double r,Point &result1,Point &result2) {

double line=sqrt((poi.x-o.x)*(poi.x-o.x)+(poi.y-o.y)*(poi.y-o.y));

double angle=acos(r/line);

Point unitvector,lin;

lin.x=poi.x-o.x;

lin.y=poi.y-o.y;

unitvector.x=lin.x/sqrt(lin.x*lin.x+lin.y*lin.y)*r;

unitvector.y=lin.y/sqrt(lin.x*lin.x+lin.y*lin.y)*r;

result1=Rotate(unitvector,-angle);

result2=Rotate(unitvector,angle);

result1.x+=o.x;

result1.y+=o.y;

result2.x+=o.x;

result2.y+=o.y;

return ;

}

7三維幾何（未）

//三維幾何函式庫

#define zero(x)(((x)>0?(x):-(x))<eps)

struct point3 {

double x,y,z;

};

struct line3 {

point3 a,b;

};

struct plane3 {

point3 a,b,c;

};

//計算dot product U .V

double dmult(point3 u,point3 v) {

return u.x*v.x+u.y*v.y+u.z*v.z;

}

//計算cross product Ux V

point3 xmult(point3 u,point3 v) {

point3 ret;

ret.x=u.y*v.z-v.y*u.z;

ret.y=u.z*v.x-u.x*v.z;

ret.z=u.x*v.y-u.y*v.x;

return ret;

}

//向量差 U - V

point3 subt(point3 u,point3 v) {

point3 ret;

ret.x=u.x-v.x;

ret.y=u.y-v.y;

ret.z=u.z-v.z;

return ret;

}

//取平面法向量

point3 pvec(plane3 s) {

return xmult(subt(s.a,s.b),subt(s.b,s.c));

}

point3 pvec(point3 s1,point3 s2,point3 s3) {

return xmult(subt(s1,s2),subt(s2,s3));

}

//兩點距離,單引數取向量大小

double dist(point3 p1,point3 p2) {

return sqrt((p1.x-p2.x)*(p1.x-p2.x)+(p1.y-p2.y)*(p1.y-p2.y)+(p1.z-p2.z)*(p1.z-p2.z));

}

//向量大小

double vlen(point3 p) {

return sqrt(p.x*p.x+p.y*p.y+p.z*p.z);

}

//判三點共線

int dots_inline(point3 p1,point3 p2,point3p3) {

return vlen(xmult(subt(p1,p2),subt(p2,p3)))<eps;

}

//判四點共面

int dots_onplane(point3 a,point3 b,point3c,point3 d) {

return zero(dmult(pvec(a,b,c),subt(d,a)));

}

//判點是否線上段上,包括端點和共線

int dot_online_in(point3 p,line3 l) {

returnzero(vlen(xmult(subt(p,l.a),subt(p,l.b))))&&(l.a.x-p.x)*(l.b.x-p.x)<eps&&

(l.a.y-p.y)*(l.b.y-p.y)<eps&&(l.a.z-p.z)*(l.b.z-p.z)<eps;

}

int dot_online_in(point3 p,point3 l1,point3l2) {

returnzero(vlen(xmult(subt(p,l1),subt(p,l2))))&&(l1.x-p.x)*(l2.x-p.x)<eps&&

(l1.y-p.y)*(l2.y-p.y)<eps&&(l1.z-p.z)*(l2.z-p.z)<eps;

}

//判點是否線上段上,不包括端點

int dot_online_ex(point3 p,line3 l) {

returndot_online_in(p,l)&&(!zero(p.x-l.a.x)||!zero(p.y-l.a.y)||!zero(p.z-l.a.z))&&

(!zero(p.x-l.b.x)||!zero(p.y-l.b.y)||!zero(p.z-l.b.z));

}

int dot_online_ex(point3 p,point3 l1,point3l2) {

returndot_online_in(p,l1,l2)&&(!zero(p.x-l1.x)||!zero(p.y-l1.y)||!zero(p.z-l1.z))&&

(!zero(p.x-l2.x)||!zero(p.y-l2.y)||!zero(p.z-l2.z));

}

//判點是否在空間三角形上,包括邊界,三點共線無意義

int dot_inplane_in(point3 p,plane3 s) {

return zero(vlen(xmult(subt(s.a,s.b),subt(s.a,s.c)))-vlen(xmult(subt(p,s.a),subt(p,s.b)))-

vlen(xmult(subt(p,s.b),subt(p,s.c)))-vlen(xmult(subt(p,s.c),subt(p,s.a))));

}

int dot_inplane_in(point3 p,point3s1,point3 s2,point3 s3) {

return zero(vlen(xmult(subt(s1,s2),subt(s1,s3)))-vlen(xmult(subt(p,s1),subt(p,s2)))-

vlen(xmult(subt(p,s2),subt(p,s3)))-vlen(xmult(subt(p,s3),subt(p,s1))));

}

//判點是否在空間三角形上,不包括邊界,三點共線無意義

int dot_inplane_ex(point3 p,plane3 s) {

returndot_inplane_in(p,s)&&vlen(xmult(subt(p,s.a),subt(p,s.b)))>eps&&

vlen(xmult(subt(p,s.b),subt(p,s.c)))>eps&&vlen(xmult(subt(p,s.c),subt(p,s.a)))>eps;

}

int dot_inplane_ex(point3 p,point3s1,point3 s2,point3 s3) {

return dot_inplane_in(p,s1,s2,s3)&&vlen(xmult(subt(p,s1),subt(p,s2)))>eps&&

vlen(xmult(subt(p,s2),subt(p,s3)))>eps&&vlen(xmult(subt(p,s3),subt(p,s1)))>eps;

}

//判兩點線上段同側,點線上段上返回0,不共面無意義

int same_side(point3 p1,point3 p2,line3 l) {

return dmult(xmult(subt(l.a,l.b),subt(p1,l.b)),xmult(subt(l.a,l.b),subt(p2,l.b)))>eps;

}

int same_side(point3 p1,point3 p2,point3l1,point3 l2) {

returndmult(xmult(subt(l1,l2),subt(p1,l2)),xmult(subt(l1,l2),subt(p2,l2)))>eps;

}

//判兩點線上段異側,點線上段上返回0,不共面無意義

int opposite_side(point3 p1,point3 p2,line3l) {

returndmult(xmult(subt(l.a,l.b),subt(p1,l.b)),xmult(subt(l.a,l.b),subt(p2,l.b)))<-eps;

}

int opposite_side(point3 p1,point3p2,point3 l1,point3 l2) {

returndmult(xmult(subt(l1,l2),subt(p1,l2)),xmult(subt(l1,l2),subt(p2,l2)))<-eps;

}

//判兩點在平面同側,點在平面上返回0

int same_side(point3 p1,point3 p2,plane3 s){

return dmult(pvec(s),subt(p1,s.a))*dmult(pvec(s),subt(p2,s.a))>eps;

}

int same_side(point3 p1,point3 p2,point3s1,point3 s2,point3 s3) {

returndmult(pvec(s1,s2,s3),subt(p1,s1))*dmult(pvec(s1,s2,s3),subt(p2,s1))>eps;

}

//判兩點在平面異側,點在平面上返回0

int opposite_side(point3 p1,point3p2,plane3 s) {

return dmult(pvec(s),subt(p1,s.a))*dmult(pvec(s),subt(p2,s.a))<-eps;

}

int opposite_side(point3 p1,point3 p2,point3s1,point3 s2,point3 s3) {

returndmult(pvec(s1,s2,s3),subt(p1,s1))*dmult(pvec(s1,s2,s3),subt(p2,s1))<-eps;

}

//判兩直線平行

int parallel(line3 u,line3 v) {

return vlen(xmult(subt(u.a,u.b),subt(v.a,v.b)))<eps;

}

int parallel(point3 u1,point3 u2,point3v1,point3 v2) {

return vlen(xmult(subt(u1,u2),subt(v1,v2)))<eps;

}

//判兩平面平行

int parallel(plane3 u,plane3 v) {

return vlen(xmult(pvec(u),pvec(v)))<eps;

}

int parallel(point3 u1,point3 u2,point3 u3,point3v1,point3 v2,point3 v3) {

return vlen(xmult(pvec(u1,u2,u3),pvec(v1,v2,v3)))<eps;

}

//判直線與平面平行

int parallel(line3 l,plane3 s) {

return zero(dmult(subt(l.a,l.b),pvec(s)));

}

int parallel(point3 l1,point3 l2,point3s1,point3 s2,point3 s3) {

return zero(dmult(subt(l1,l2),pvec(s1,s2,s3)));

}

//判兩直線垂直

int perpendicular(line3 u,line3 v) {

return zero(dmult(subt(u.a,u.b),subt(v.a,v.b)));

}

int perpendicular(point3 u1,point3u2,point3 v1,point3 v2) {

return zero(dmult(subt(u1,u2),subt(v1,v2)));

}

//判兩平面垂直

int perpendicular(plane3 u,plane3 v)

{

return zero(dmult(pvec(u),pvec(v)));

}

int perpendicular(point3 u1,point3u2,point3 u3,point3 v1,point3 v2,point3 v3)

{

return zero(dmult(pvec(u1,u2,u3),pvec(v1,v2,v3)));

}

//判直線與平面平行

int perpendicular(line3 l,plane3 s)

{

return vlen(xmult(subt(l.a,l.b),pvec(s)))<eps;

}

int perpendicular(point3 l1,point3 l2,point3s1,point3 s2,point3 s3)

{

return vlen(xmult(subt(l1,l2),pvec(s1,s2,s3)))<eps;

}

//判兩線段相交,包括端點和部分重合

int intersect_in(line3 u,line3 v)

{

if (!dots_onplane(u.a,u.b,v.a,v.b))

return 0;

if (!dots_inline(u.a,u.b,v.a)||!dots_inline(u.a,u.b,v.b))

return !same_side(u.a,u.b,v)&&!same_side(v.a,v.b,u);

return dot_online_in(u.a,v)||dot_online_in(u.b,v)||dot_online_in(v.a,u)||dot_online_in(v.b,u);

}

int intersect_in(point3 u1,point3 u2,point3v1,point3 v2)

{

if (!dots_onplane(u1,u2,v1,v2))

return 0;

if (!dots_inline(u1,u2,v1)||!dots_inline(u1,u2,v2))

return!same_side(u1,u2,v1,v2)&&!same_side(v1,v2,u1,u2);

returndot_online_in(u1,v1,v2)||dot_online_in(u2,v1,v2)||dot_online_in(v1,u1,u2)||dot_online_in(v2,u1,u2);

}

//判兩線段相交,不包括端點和部分重合

int intersect_ex(line3 u,line3 v) {

returndots_onplane(u.a,u.b,v.a,v.b)&&opposite_side(u.a,u.b,v)&&opposite_side(v.a,v.b,u);

}

int intersect_ex(point3 u1,point3 u2,point3v1,point3 v2) {

returndots_onplane(u1,u2,v1,v2)&&opposite_side(u1,u2,v1,v2)&&opposite_side(v1,v2,u1,u2);

}

//判線段與空間三角形相交,包括交於邊界和(部分)包含

int intersect_in(line3 l,plane3 s)

{

return!same_side(l.a,l.b,s)&&!same_side(s.a,s.b,l.a,l.b,s.c)&&

!same_side(s.b,s.c,l.a,l.b,s.a)&&!same_side(s.c,s.a,l.a,l.b,s.b);

}

int intersect_in(point3 l1,point3 l2,point3s1,point3 s2,point3 s3)

{

return!same_side(l1,l2,s1,s2,s3)&&!same_side(s1,s2,l1,l2,s3)&&

!same_side(s2,s3,l1,l2,s1)&&!same_side(s3,s1,l1,l2,s2);

}

//判線段與空間三角形相交,不包括交於邊界和(部分)包含

int intersect_ex(line3 l,plane3 s)

{

return opposite_side(l.a,l.b,s)&&opposite_side(s.a,s.b,l.a,l.b,s.c)&&

opposite_side(s.b,s.c,l.a,l.b,s.a)&&opposite_side(s.c,s.a,l.a,l.b,s.b);

}

int intersect_ex(point3 l1,point3 l2,point3s1,point3 s2,point3 s3)

{

returnopposite_side(l1,l2,s1,s2,s3)&&opposite_side(s1,s2,l1,l2,s3)&&

opposite_side(s2,s3,l1,l2,s1)&&opposite_side(s3,s1,l1,l2,s2);

}

//計算兩直線交點,注意事先判斷直線是否共面和平行!

//線段交點請另外判線段相交(同時還是要判斷是否平行!)

point3 intersection(line3 u,line3 v)

{

point3 ret=u.a;

double t=((u.a.x-v.a.x)*(v.a.y-v.b.y)-(u.a.y-v.a.y)*(v.a.x-v.b.x))

/((u.a.x-u.b.x)*(v.a.y-v.b.y)-(u.a.y-u.b.y)*(v.a.x-v.b.x));

ret.x+=(u.b.x-u.a.x)*t;

ret.y+=(u.b.y-u.a.y)*t;

ret.z+=(u.b.z-u.a.z)*t;

return ret;

}

point3 intersection(point3 u1,point3u2,point3 v1,point3 v2)

{

point3 ret=u1;

double t=((u1.x-v1.x)*(v1.y-v2.y)-(u1.y-v1.y)*(v1.x-v2.x))

/((u1.x-u2.x)*(v1.y-v2.y)-(u1.y-u2.y)*(v1.x-v2.x));

ret.x+=(u2.x-u1.x)*t;

ret.y+=(u2.y-u1.y)*t;

ret.z+=(u2.z-u1.z)*t;

return ret;

}

//計算直線與平面交點,注意事先判斷是否平行,並保證三點不共線!

//線段和空間三角形交點請另外判斷

point3 intersection(line3 l,plane3 s)

{

point3 ret=pvec(s);

double t=(ret.x*(s.a.x-l.a.x)+ret.y*(s.a.y-l.a.y)+ret.z*(s.a.z-l.a.z))/

(ret.x*(l.b.x-l.a.x)+ret.y*(l.b.y-l.a.y)+ret.z*(l.b.z-l.a.z));

ret.x=l.a.x+(l.b.x-l.a.x)*t;

ret.y=l.a.y+(l.b.y-l.a.y)*t;

ret.z=l.a.z+(l.b.z-l.a.z)*t;

return ret;

}

point3 intersection(point3 l1,point3l2,point3 s1,point3 s2,point3 s3)

{

point3 ret=pvec(s1,s2,s3);

double t=(ret.x*(s1.x-l1.x)+ret.y*(s1.y-l1.y)+ret.z*(s1.z-l1.z))/

(ret.x*(l2.x-l1.x)+ret.y*(l2.y-l1.y)+ret.z*(l2.z-l1.z));

ret.x=l1.x+(l2.x-l1.x)*t;

ret.y=l1.y+(l2.y-l1.y)*t;

ret.z=l1.z+(l2.z-l1.z)*t;

return ret;

}

//計算兩平面交線,注意事先判斷是否平行,並保證三點不共線!

line3 intersection(plane3 u,plane3 v)

{

line3 ret;

ret.a=parallel(v.a,v.b,u.a,u.b,u.c)?intersection(v.b,v.c,u.a,u.b,u.c):intersection(v.a,v.b,u.a,u.b,u.c);

ret.b=parallel(v.c,v.a,u.a,u.b,u.c)?intersection(v.b,v.c,u.a,u.b,u.c):intersection(v.c,v.a,u.a,u.b,u.c);

return ret;

}

line3 intersection(point3 u1,point3u2,point3 u3,point3 v1,point3 v2,point3 v3)

{

line3 ret;

ret.a=parallel(v1,v2,u1,u2,u3)?intersection(v2,v3,u1,u2,u3):intersection(v1,v2,u1,u2,u3);

ret.b=parallel(v3,v1,u1,u2,u3)?intersection(v2,v3,u1,u2,u3):intersection(v3,v1,u1,u2,u3);

return ret;

}

//點到直線距離

double ptoline(point3 p,line3 l) {

return vlen(xmult(subt(p,l.a),subt(l.b,l.a)))/dist(l.a,l.b);

}

double ptoline(point3 p,point3 l1,point3l2) {

return vlen(xmult(subt(p,l1),subt(l2,l1)))/dist(l1,l2);

}

//點到平面距離

double ptoplane(point3 p,plane3 s)

{

return fabs(dmult(pvec(s),subt(p,s.a)))/vlen(pvec(s));

}

double ptoplane(point3 p,point3 s1,point3s2,point3 s3)

{

return fabs(dmult(pvec(s1,s2,s3),subt(p,s1)))/vlen(pvec(s1,s2,s3));

}

//直線到直線距離

double linetoline(line3 u,line3 v)

{

point3 n=xmult(subt(u.a,u.b),subt(v.a,v.b));

return fabs(dmult(subt(u.a,v.a),n))/vlen(n);

}

double linetoline(point3 u1,point3u2,point3 v1,point3 v2)

{

point3 n=xmult(subt(u1,u2),subt(v1,v2));

return fabs(dmult(subt(u1,v1),n))/vlen(n);

}

//兩直線夾角cos值

double angle_cos(line3 u,line3 v) {

returndmult(subt(u.a,u.b),subt(v.a,v.b))/vlen(subt(u.a,u.b))/vlen(subt(v.a,v.b));

}

double angle_cos(point3 u1,point3 u2,point3v1,point3 v2) {

returndmult(subt(u1,u2),subt(v1,v2))/vlen(subt(u1,u2))/vlen(subt(v1,v2));

}

//兩平面夾角cos值

double angle_cos(plane3 u,plane3 v) {

return dmult(pvec(u),pvec(v))/vlen(pvec(u))/vlen(pvec(v));

}

double angle_cos(point3 u1,point3 u2,point3u3,point3 v1,point3 v2,point3 v3) {

returndmult(pvec(u1,u2,u3),pvec(v1,v2,v3))/vlen(pvec(u1,u2,u3))/vlen(pvec(v1,v2,v3));

}

//直線平面夾角sin值

double angle_sin(line3 l,plane3 s) {

return dmult(subt(l.a,l.b),pvec(s))/vlen(subt(l.a,l.b))/vlen(pvec(s));

}

double angle_sin(point3 l1,point3 l2,point3s1,point3 s2,point3 s3) {

returndmult(subt(l1,l2),pvec(s1,s2,s3))/vlen(subt(l1,l2))/vlen(pvec(s1,s2,s3));

}

【總結】計算幾何模板

相關文章