计算几何模板

正文

 

㈠ 点的基本运算 
1. 平面上两点之间距离 1 
2. 判断两点是否重合 1 
3. 矢量叉乘 1 
4. 矢量点乘 2 
5. 判断点是否在线段上 2 
6. 求一点饶某点旋转后的坐标 2 
7. 求矢量夹角 2

㈡ 线段及直线的基本运算 
1. 点与线段的关系 3 
2. 求点到线段所在直线垂线的垂足 4 
3. 点到线段的最近点 4 
4. 点到线段所在直线的距离 4 
5. 点到折线集的最近距离 4 
6. 判断圆是否在多边形内 5 
7. 求矢量夹角余弦 5 
8. 求线段之间的夹角 5 
9. 判断线段是否相交 6 
10.判断线段是否相交但不交在端点处 6 
11.求线段所在直线的方程 6 
12.求直线的斜率 7 
13.求直线的倾斜角 7 
14.求点关于某直线的对称点 7 
15.判断两条直线是否相交及求直线交点 7 
16.判断线段是否相交,如果相交返回交点 7

㈢ 多边形常用算法模块 
1. 判断多边形是否简单多边形 8 
2. 检查多边形顶点的凸凹性 9 
3. 判断多边形是否凸多边形 9 
4. 求多边形面积 9 
5. 判断多边形顶点的排列方向,方法一 10 
6. 判断多边形顶点的排列方向,方法二 10 
7. 射线法判断点是否在多边形内 10 
8. 判断点是否在凸多边形内 11 
9. 寻找点集的graham算法 12 
10.寻找点集凸包的卷包裹法 13 
11.判断线段是否在多边形内 14 
12.求简单多边形的重心 15 
13.求凸多边形的重心 17 
14.求肯定在给定多边形内的一个点 17 
15.求从多边形外一点出发到该多边形的切线 18 
16.判断多边形的核是否存在 19

㈣ 圆的基本运算 
1 .点是否在圆内 20 
2 .求不共线的三点所确定的圆 21

㈤ 矩形的基本运算 
1.已知矩形三点坐标,求第4点坐标 22

㈥ 常用算法的描述 22

㈦ 补充 
1.两圆关系: 24 
2.判断圆是否在矩形内: 24 
3.点到平面的距离: 25 
4.点是否在直线同侧: 25 
5.镜面反射线: 25 
6.矩形包含: 26 
7.两圆交点: 27 
8.两圆公共面积: 28 
9. 圆和直线关系: 29 
10. 内切圆: 30 
11. 求切点: 31 
12. 线段的左右旋: 31 
13.公式: 32


/* 需要包含的头文件 */ 
#include

/* 常用的常量定义 */ 
const double INF = 1E200 
const double EP = 1E-10 
const int MAXV = 300 
const double PI = 3.14159265

/* 基本几何结构 */ 
struct POINT 

double x; 
double y; POINT(double a=0, double b=0) { x=a; y=b;} file://constructor/ 
}; 
struct LINESEG 

POINT s; 
POINT e; LINESEG(POINT a, POINT b) { s=a; e=b;} 
LINESEG() { } 
}; 
struct LINE // 直线的解析方程 a*x+b*y+c=0 为统一表示,约定 a >= 0 

double a; 
double b; 
double c; LINE(double d1=1, double d2=-1, double d3=0) {a=d1; b=d2; c=d3;} 
};

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

double dist(POINT p1,POINT p2) // 返回两点之间欧氏距离 

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

bool equal_point(POINT p1,POINT p2) // 判断两个点是否重合 

return ( (abs(p1.x-p2.x)}

/*****************************************************************************

r=multiply(sp,ep,op),得到(sp-op)*(ep-op)的叉积 
r>0:ep在矢量opsp的逆时针方向; 
r=0:opspep三点共线; 
r<0:ep在矢量opsp的顺时针方向 
******************************************************************************
*/

double multiply(POINT sp,POINT ep,POINT op) 

return((sp.x-op.x)*(ep.y-op.y)-(ep.x-op.x)*(sp.y-op.y)); 
}

/*****************************************************************************
** 
r=dotmultiply(p1,p2,op),得到矢量(p1-op)和(p2-op)的点积,如果两个矢量都非零矢量
 
r<0:两矢量夹角为锐角;r=0:两矢量夹角为直角;r>0:两矢量夹角为钝角 
******************************************************************************
*/ 
double dotmultiply(POINT p1,POINT p2,POINT p0) 

return ((p1.x-p0.x)*(p2.x-p0.x)+(p1.y-p0.y)*(p2.y-p0.y)); 
}

/* 判断点p是否在线段l上,条件:(p在线段l所在的直线上)&& (点p在以线段l为对角线的
矩形内) */ 
bool online(LINESEG l,POINT p) 

return((multiply(l.e,p,l.s)==0) 
&&( ( (p.x-l.s.x)*(p.x-l.e.x)<=0 )&&( (p.y-l.s.y)*(p.y-l.e.y)<=0 ) ) ); 
}

// 返回点p以点o为圆心逆时针旋转alpha(单位:弧度)后所在的位置 
POINT rotate(POINT o,double alpha,POINT p) 

POINT tp; 
p.x-=o.x; 
p.y-=o.y; 
tp.x=p.x*cos(alpha)-p.y*sin(alpha)+o.x; 
tp.y=p.y*cos(alpha)+p.x*sin(alpha)+o.y; 
return tp; 
}

/* 返回顶角在o点,起始边为os,终止边为oe的夹角(单位:弧度) 
角度小于pi,返回正值 
角度大于pi,返回负值 
可以用于求线段之间的夹角 
*/ 
double angle(POINT o,POINT s,POINT e) 

double cosfi,fi,norm; 
double dsx = s.x - o.x; 
double dsy = s.y - o.y; 
double dex = e.x - o.x; 
double dey = e.y - o.y;

cosfi=dsx*dex+dsy*dey; 
norm=(dsx*dsx+dey*dey)*(dex*dex+dey*dey); 
cosfi /= sqrt( norm );

if (cosfi >= 1.0 ) return 0; 
if (cosfi <= -1.0 ) return -3.1415926;

fi=acos(cosfi); 
if (dsx*dey-dsy*dex>0) return fi; // 说明矢量os 在矢量 oe的顺时针方向 
return -fi; 
}


/*****************************\ 
* * 
* 线段及直线的基本运算 * 
* * 
\*****************************/

/* 判断点与线段的关系,用途很广泛 
本函数是根据下面的公式写的,P是点C到线段AB所在直线的垂足

AC dot AB 
r = --------- 
||AB||^2 
(Cx-Ax)(Bx-Ax) + (Cy-Ay)(By-Ay) 
= ------------------------------- 
L^2

r has the following meaning:

r=0 P = A 
r=1 P = B 
r<0 P is on the backward extension of AB 
r>1 P is on the forward extension of AB 
0*/ 
double relation(POINT p,LINESEG l) 

LINESEG tl; 
tl.s=l.s; 
tl.e=p; 
return dotmultiply(tl.e,l.e,l.s)/(dist(l.s,l.e)*dist(l.s,l.e)); 
}

// 求点C到线段AB所在直线的垂足 P 
POINT perpendicular(POINT p,LINESEG l) 

double r=relation(p,l); 
POINT tp; 
tp.x=l.s.x+r*(l.e.x-l.s.x); 
tp.y=l.s.y+r*(l.e.y-l.s.y); 
return tp; 

/* 求点p到线段l的最短距离,并返回线段上距该点最近的点np 
注意:np是线段l上到点p最近的点,不一定是垂足 */ 
double ptolinesegdist(POINT p,LINESEG l,POINT &np) 

double r=relation(p,l); 
if(r<0) 

np=l.s; 
return dist(p,l.s); 

if(r>1) 

np=l.e; 
return dist(p,l.e); 

np=perpendicular(p,l); 
return dist(p,np); 
}

// 求点p到线段l所在直线的距离,请注意本函数与上个函数的区别 
double ptoldist(POINT p,LINESEG l) 

return abs(multiply(p,l.e,l.s))/dist(l.s,l.e); 
}

/* 计算点到折线集的最近距离,并返回最近点. 
注意:调用的是ptolineseg()函数 */ 
double ptopointset(int vcount,POINT pointset[],POINT p,POINT &q) 

int i; 
double cd=double(INF),td; 
LINESEG l; 
POINT tq,cq;

for(i=0;i{ 
l.s=pointset[i]; 
l.e=pointset[i+1]; 
td=ptolinesegdist(p,l,tq); 
if(td{ 
cd=td; 
cq=tq; 


q=cq; 
return cd; 

/* 判断圆是否在多边形内.ptolineseg()函数的应用2 */ 
bool CircleInsidePolygon(int vcount,POINT center,double radius,POINT polygon[]


POINT q; 
double d; 
q.x=0; 
q.y=0; 
d=ptopointset(vcount,polygon,center,q); 
if(dreturn true; 
else 
return false; 
}

/* 返回两个矢量l1和l2的夹角的余弦(-1 --- 1)注意:如果想从余弦求夹角的话,注意反
余弦函数的定义域是从 0到pi */ 
double cosine(LINESEG l1,LINESEG l2) 

return (((l1.e.x-l1.s.x)*(l2.e.x-l2.s.x) + 
(l1.e.y-l1.s.y)*(l2.e.y-l2.s.y))/(dist(l1.e,l1.s)*dist(l2.e,l2.s))) ); 

// 返回线段l1与l2之间的夹角 单位:弧度 范围(-pi,pi) 
double lsangle(LINESEG l1,LINESEG l2) 

POINT o,s,e; 
o.x=o.y=0; 
s.x=l1.e.x-l1.s.x; 
s.y=l1.e.y-l1.s.y; 
e.x=l2.e.x-l2.s.x; 
e.y=l2.e.y-l2.s.y; 
return angle(o,s,e); 

// 如果线段u和v相交(包括相交在端点处)时,返回true 
bool intersect(LINESEG u,LINESEG v) 

return( (max(u.s.x,u.e.x)>=min(v.s.x,v.e.x))&& file://排斥实验 
(max(v.s.x,v.e.x)>=min(u.s.x,u.e.x))&& 
(max(u.s.y,u.e.y)>=min(v.s.y,v.e.y))&& 
(max(v.s.y,v.e.y)>=min(u.s.y,u.e.y))&& 
(multiply(v.s,u.e,u.s)*multiply(u.e,v.e,u.s)>=0)&& file://跨立实验 
(multiply(u.s,v.e,v.s)*multiply(v.e,u.e,v.s)>=0)); 
}


// (线段u和v相交)&&(交点不是双方的端点) 时返回true 
bool intersect_A(LINESEG u,LINESEG v) 

return((intersect(u,v))&& 
(!online(u,v.s))&& 
(!online(u,v.e))&& 
(!online(v,u.e))&& 
(!online(v,u.s))); 
}


// 线段v所在直线与线段u相交时返回true;方法:判断线段u是否跨立线段v 
bool intersect_l(LINESEG u,LINESEG v) 

return multiply(u.s,v.e,v.s)*multiply(v.e,u.e,v.s)>=0; 
}

 

Given a polygon, your task is to find the centre of gravity for the given polygon. 

#include<stdio.h>
#include<math.h>
struct point
{
	double x,y;
}pp[1000008];
point becenter(point pnt[],int m)//借的模板;
{
	point p,s;
	double tp,area=0,tpx=0,tpy=0;
	int i;
	p.x=pnt[0].x,p.y=pnt[0].y;
	for(i=1;i<=m;i++)
	{
		s.x=pnt[i%m].x;
		s.y=pnt[i%m].y;
		tp=(p.x*s.y-s.x*p.y);
		area+=tp/2;
		tpx+=(p.x+s.x)*tp;
		tpy+=(p.y+s.y)*tp;
		p.x=s.x;p.y=s.y;
	}
	s.x=tpx/(6*area);
	s.y=tpy/(6*area);
	//printf("tpx=%.2lf tpy=%.2lf area=%.2lf\n",tpx,tpy,area);
	//printf("%.2lf %.2lf\n",s.x,s.y);
	return s;
}
int main()
{
	int i,j;
	double l,k,d,n;
	int ncase,m;
	scanf("%d",&ncase);
	while(ncase--)
	{
		scanf("%d",&m);
		for(i=0;i<m;i++)
		scanf("%lf%lf",&pp[i].x,&pp[i].y);
		point c=becenter(pp,m);
		printf("%.2lf %.2lf\n",c.x,c.y);
	}
	return 0;
}

多边形与圆交面积

double sign( double x ){
      if( x < -eps ) return -1;
      return x > eps;
}

struct point{
      double x, y ;
      point(){}
      point(double x, double y):x(x),y(y){}
      double operator *(const point &b)const{
            return x * b.y - y * b.x;
      }
      point operator -(const point &b)const{
            return point( x - b.x, y - b.y );
      }
      void in(){
              scanf("%lf%lf",&x,&y);
      }
}p[maxn] , q[maxn], o;
double r ,ans ; // 已原点为圆心的圆半径

圆o与线段x-y的交点,,如果存在顺序存入p[0]和p[1],,返回交点的个数
int xianduan_jiao_yuan( point o, double r,point u, point v, point &o1, point &o2){
      u = u - o;   v = v - o;
      point   w = v - u ;
      double a = w.x * w.x + w.y * w.y ;
      double b = w.x * u.x * 2 + w.y * u.y * 2;
      double c = u.x * u.x + u.y * u.y - r * r;
      double d = b * b - 4 * a * c;
      if( sign( d ) < 0) return 0;
      d = sqrt( d );
      double t1 = ( - b + d )/2/a, t2 = ( - b - d )/2/a;
      if( t1 > t2 )   swap( t1, t2 );
      int ans = 0;
      注意这里端点在圆上没有算作相交,有需要改成等号既可
      if( sign( t1 ) > 0 && sign( t1 - 1 ) < 0 ){
            ans ++ ;
            o1 = point( u.x + w.x * t1 + o.x, u.y + w.y*t1+o.y);
      }注意重根这里算作了一个交点,,
      if(sign( t1 - t2) != 0 && sign(t2)>0&&sign(t2-1)<0){
            ans ++;
            if( ans == 1 ) o1 = point( u.x+w.x*t2+o.x,u.y+w.y*t2+o.y );
            else o2 = point( u.x+w.x*t2+o.x,u.y+w.y*t2+o.y );
      }
      return ans ;
}

double get( point p[], int n ){
        double ans = 0;
        for( int i = 1; i < n; i ++ ){
               point a, b;
               a = p[ i - 1]; b = p[ i ];//只要有一个人在圆外面 那就只需要计算扇形的面积即可
               if( sqrt(a.x*a.x+a.y*a.y) - r > eps || sqrt( b.x * b.x + b.y*b.y) - r > eps ){
                     double t = atan2(b.y , b.x ) - atan2(a.y , a.x ) ;
                     while( t - pi > eps   ) t -= pi * 2 ;//t的范围是( -pi , pi )
                     while( t + pi < -eps ) t += pi * 2 ;
                     ans += t * r * r ;
               }else ans += a * b ;   //叉乘   这个在圆里面
        }
        return ans ;
}

void add( point a , point b ){
        point p1, p2;
        int end = 0;
        p[end++] = a;
        point o1, o2;
        int t = xianduan_jiao_yuan(point(0,0),r,a,b,o1,o2);
        if( t >= 1) p[end++] = o1;
        if( t >= 2) p[end++] = o2;
        p[end++] = b;
        ans += get( p, end );
}

int main(){
        //o:圆心 
        ans=0;
             scanf("%d",&n);
            for(   i = 0; i < n; i ++ ){
                  scanf("%lf%lf",&q[i].x, &q[i].y );
                  q[i] = q[i] - o;
            }
            q[n] = q[0];
            for( i   =0; i < n; i++ ){
                  add( q[i], q[i+1] );
            }
            printf("%.2lf\n",fabs(ans)/2);

      return 0;
}

struct point1{
    double x, y ;
    point1(){}
    point1(double x, double y):x(x),y(y){}
    double operator *(const point1 &b)const{
        return x * b.y - y * b.x;
    }
    point1 operator -(const point1 &b)const{
        return point1( x - b.x, y - b.y );
    }
    void in(){
         scanf("%lf%lf",&x,&y);
    }
}p[maxn] , q[maxn], o;

const double eps = 1e-10;
inline double max (double a, double b) { if (a > b) return a; else return b; }
inline double min (double a, double b) { if (a < b) return a; else return b; }
inline int fi (double a)
{
    if (a > eps) return 1;
    else if (a >= -eps) return 0;
    else return -1;
}
class vector
{
public:
    double x, y;
    vector (void) {}
    vector (double x0, double y0) : x(x0), y(y0) {}
    double operator * (const vector& a) const { return x * a.y - y * a.x; }
    double operator % (const vector& a) const { return x * a.x + y * a.y; }
    vector verti (void) const { return vector(-y, x); }
    double length (void) const { return sqrt(x * x + y * y); }
    vector adjust (double len)
    {
        double ol = len / length();
        return vector(x * ol, y * ol);
    }
    vector oppose (void) { return vector(-x, -y); }
};
class point
{
public:
    double x, y;
    point (void) {}
    point (double x0, double y0) : x(x0), y(y0) {}
    vector operator - (const point& a) const { return vector(x - a.x, y - a.y); }
    point operator + (const vector& a) const { return point(x + a.x, y + a.y); }
};
class segment
{
public:
    point a, b;
    segment (void) {}
    segment (point a0, point b0) : a(a0), b(b0) {}
    point intersect (const segment& s) const
    {
        vector v1 = s.a - a, v2 = s.b - a, v3 = s.b - b, v4 = s.a - b;
        double s1 = v1 * v2, s2 = v3 * v4;
        double se = s1 + s2;
        s1 /= se, s2 /= se;
        return point(a.x * s2 + b.x * s1, a.y * s2 + b.y * s1);
    }
    point pverti (const point& p) const
    {
        vector t = (b - a).verti();
        segment uv(p, p + t);
        return intersect(uv);
    }
    bool on_segment (const point& p) const
    {
        if (fi(min(a.x, b.x) - p.x) <= 0 && fi(p.x - max(a.x, b.x)) <= 0 &&
            fi(min(a.y, b.y) - p.y) <= 0 && fi(p.y - max(a.y, b.y)) <= 0) return true;
        else return false;
    }
};
double radius;
point polygon[550];
double kuras_area (point a, point b)
{
    point ori(0, 0);
    int sgn = fi((b - ori) * (a - ori));
    double da = (a - ori).length(), db = (b - ori).length();
    int ra = fi(da - radius), rb = fi(db - radius);
    double angle = acos(((b - ori) % (a - ori)) / (da * db));
    segment t(a, b); point h, u; vector seg;
    double ans, dlt, mov, tangle;

    if (fi(da) == 0 || fi(db) == 0) return 0;
    else if (sgn == 0) return 0;
    else if (ra <= 0 && rb <= 0) return fabs((b - ori) * (a - ori)) / 2 * sgn;
    else if (ra >= 0 && rb >= 0)
    {
        h = t.pverti(ori);
        dlt = (h - ori).length();
        if (!t.on_segment(h) || fi(dlt - radius) >= 0)
            return radius * radius * (angle / 2) * sgn;
        else
        {
            ans = radius * radius * (angle / 2);
            tangle = acos(dlt / radius);
            ans -= radius * radius * tangle;
            ans += radius * sin(tangle) * dlt;
            return ans * sgn;
        }
    }
    else
    {
        h = t.pverti(ori);
        dlt = (h - ori).length();
        seg = b - a;
        mov = sqrt(radius * radius - dlt * dlt);
        seg = seg.adjust(mov);
        if (t.on_segment(h + seg)) u = h + seg;
        else u = h + seg.oppose();
        if (ra == 1) swap(a, b);
        ans = fabs((a - ori) * (u - ori)) / 2;
        tangle = acos(((u - ori) % (b - ori)) / ((u - ori).length() * (b - ori).length()));
        ans += radius * radius * (tangle / 2);
        return ans * sgn;
    }
}

int main()
{
       //o:圆心
        double x,y;
        scanf("%d",&n);
        for (i = 0; i < n; i++)
        {
            scanf("%lf %lf", &x, &y);
            x-=o.x;
            y-=o.y;
            polygon[i] = point(x, y);
        }
        double area = 0;
        for (int i = 0; i < n; i++)
            area += kuras_area(polygon[i], polygon[(i + 1) % n]);
        printf("%.2f\n", fabs(area));

    return 0;
}


  • 0
    点赞
  • 0
    收藏
    觉得还不错? 一键收藏
  • 0
    评论
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值