计算几何模板

标签：

34.过圆外一点的直线与圆的两个切点

<span style="font-size:14px;">#include<iostream>
#include<cmath>
#include<cstdio>
#define INF 1E200
using namespace std;
const double PI  = 3.14159265;
const double eps=0.0000000001;
struct Point //点
{
    double x,y;
    Point (double a=0,double b=0):x(a),y(b) {}
};

struct Line_segment //线段
{
    Point s,e;
    Line_segment() {}
    Line_segment(Point a,Point b):s(a),e(b) {}
};

struct Line //直线
{
    double A,B,C;
    Line(double A=1,double B=-1,double C=0):A(A),B(B),C(C) {}
};

inline double Max(double a,double b)
{
    return a>b?a:b;
}

inline double Min(double a,double b)
{
    return a<b?a:b;
}

//计算几何 开始!
double Dist(Point a,Point b) //1.两点之间距离
{
    return sqrt((a.x-b.x)*(a.x-b.x)+(a.y-b.y)*(a.y-b.y));
}

bool equal_Point(Point a,Point b) //2.判断两点是否重合
{
    return fabs(a.x-b.x)<eps&&fabs(a.y-b.y)<eps;
}

double multiply(Point sp,Point  ep,Point op) //3.叉积
{
    /******************************************************************************
    r=multiply(sp,ep,op),得到(sp-op) 和 ( ep-op)的叉积
    r>0; ep在矢量op sp的逆时针方向；
    r=0;op sp ep 三点共线；
    r<0;ep在矢量op sp的顺时针方向
    *******************************************************************************/
    return((sp.x-op.x)*(ep.y-op.y)-(ep.x-op.x)*(sp.y-op.y));
}

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

bool online(Line_segment l,Point p)//5.判断点p是否在线段l上 +3
{
    /******************************************************************************
    判断点p是否在线段l上
    条件：(p在线段l所在的直线上) && (点p在以线段l为对角线的矩形内)
    *******************************************************************************/
    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 ) ) );
}

Point Rotate(Point o,double alpha,Point p) // 6.返回点p以点o为圆心逆时针旋转alpha(单位：弧度)后所在的位置
{
    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;
}

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

    原理：
    r = dotmultiply(s,e,o) / (dist(o,s)*dist(o,e))
    r'= multiply(s,e,o)
    r >= 1 angle = 0;
    r <= -1 angle = -PI
    -1<r<1 && r'>0 angle = arccos(r)
    -1<r<1 && r'<=0 angle = -arccos(r)
    ********************************************************************************/
    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+dsy*dsy)*(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;
}


double relation(Point p,Line_segment l) //8.判断点与线段的关系
{
    /******************************************************************************
    判断点与线段的关系,用途很广泛
    本函数是根据下面的公式写的，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<r<1  P is interior to AB
    ********************************************************************************/
    Line_segment 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));
}




Point perpendicular(Point p,Line_segment l) //9.求点C到线段AB所在直线的垂足 P
{
    /******************************************************************************
       求点C到线段AB所在直线的垂足 P
       *******************************************************************************/
    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;
}

double ptolinesegdist(Point p,Line_segment l,Point &np) //10.求点p到线段l的最短距离,并返回线段上距该点最近的点np
{
    /******************************************************************************
       求点p到线段l的最短距离,并返回线段上距该点最近的点np
    注意：np是线段l上到点p最近的点，不一定是垂足

    *******************************************************************************/
    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);
}

double ptoldist(Point p,Line_segment l) // 11.求点p到线段l所在直线的距离,请注意本函数与上个函数的区别
{
    return abs(multiply(p,l.e,l.s))/Dist(l.s,l.e);
}


double ptopointset(int vcount,Point pointset[],Point p,Point &q) //12.计算点到折线集的最近距离,并返回最近点.
{
    /******************************************************************************
     计算点到折线集的最近距离,并返回最近点.
    注意：调用的是ptolineseg()函数
    *******************************************************************************/
    int i;
    double cd=double(INF),td;
    Line_segment l;
    Point tq,cq;

    for(i=0; i<vcount-1; i++)
    {
        l.s=pointset[i];

        l.e=pointset[i+1];
        td=ptolinesegdist(p,l,tq);
        if(td<cd)
        {
            cd=td;
            cq=tq;
        }
    }
    q=cq;
    return cd;
}

bool CircleInsidePolygon(int vcount,Point center,double radius,Point polygon[]) //13.判断圆是否在多边形内.ptolineseg()函数的应用
{
    Point q;
    double d;
    q.x=0;
    q.y=0;
    d=ptopointset(vcount,polygon,center,q);
    if(d<radius||fabs(d-radius)<eps)
        return true;
    else
        return false;
}

double cosine(Line_segment l1,Line_segment l2)   //14.返回两个 矢量 l1和l2的夹角的余弦
{
    /******************************************************************************
    返回两个矢量l1和l2的夹角的余弦(-1 --- 1)注意：如果想从余弦求夹角的话，注意反余弦函数的定义域是从 0到pi
    *******************************************************************************/
    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)));
}

double lsangle(Line_segment l1,Line_segment l2)  // 15.返回线段l1与l2之间的夹角 单位:弧度 范围(-pi，pi)
{
    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);
}


bool intersect(Line_segment u,Line_segment v) // 16.如果线段u和v相交(包括相交在端点处)时,返回true
{
    /******************************************************************************
    如果线段u和v相交(包括相交在端点处)时,返回true
    判断P1P2跨立Q1Q2的依据是: ( P1 - Q1 ) x ( Q2 - Q1 ) x ( Q2 - Q1 ) x ( P2 - Q1 ) >= 0
    判断Q1Q2跨立P1P2的依据是: ( Q1 - P1 ) x ( P2 -  P1 ) x ( P2 - P1 ) x ( Q2 - P1 ) >= 0
    *******************************************************************************/
    return((Max(u.s.x,u.e.x)>=Min(v.s.x,v.e.x))&&                     //排斥实验
           (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)&&         //跨立实验
           (multiply(u.s,v.e,v.s)*multiply(v.e,u.e,v.s)>=0));
}


bool intersect_A(Line_segment u,Line_segment v) // 17.(线段u和v相交)&&(交点不是端点)时返回true
{
    return ((intersect(u,v))&&
            (!online(u,v.s))&&
            (!online(u,v.e))&&
            (!online(v,u.e))&&
            (!online(v,u.s)));
}

bool intersect_l(Line_segment u,Line_segment v)// 18.线段v所在直线与线段u相交时返回true
{
    /******************************************************************************
    线段v所在直线与线段u相交时返回true;方法;判断线段u是否跨立线段v
    *******************************************************************************/
    return multiply(u.s,v.e,v.s)*multiply(v.e,u.e,v.s)>=0;
}


Line makeline(Point p1,Point p2)   // 19.根据已知两点坐标，求过这两点的直线解析方程
{
    /******************************************************************************
    根据已知两点坐标，求过这两点的直线解析方程Ax+By+C=0 (A>=0)
    *******************************************************************************/
    Line tl;
    int sign = 1;
    tl.A=p2.y-p1.y;
    if(tl.A<0)
    {
        sign = -1;
        tl.A=sign*tl.A;
    }
    tl.B=sign*(p1.x-p2.x);
    tl.C=sign*(p1.y*p2.x-p1.x*p2.y);
    return tl;
}


double slope(Line l)   // 20.根据直线解析方程返回直线的斜率k,  水平线返回 0, 竖直线返回 1e200
{
    if(abs(l.A) < 1e-20)
        return 0;
    if(abs(l.B) < 1e-20)
        return INF;
    return -(l.A/l.B);
}


double alpha(Line l)   //21. 返回直线的倾斜角 alpha ( 0 - pi)
{
    if(abs(l.A)< eps)
        return 0;
    if(abs(l.B)< eps)
        return PI/2;
    double k=slope(l);
    if(k>0)
        return atan(k);
    else
        return PI+atan(k);
}


Point symmetry(Line l,Point p) //22. 求点p关于直线l的对称点
{
    Point tp;
    tp.x=((l.B*l.B-l.A*l.A)*p.x-2*l.A*l.B*p.y-2*l.A*l.C)/(l.A*l.A+l.B*l.B);
    tp.y=((l.A*l.A-l.B*l.B)*p.y-2*l.A*l.B*p.x-2*l.B*l.C)/(l.A*l.A+l.B*l.B);
    return tp;
}


bool lineintersect(Line l1,Line l2,Point &p) // 23.两直线相交返回true并返回交点p,不相交则返回false
{
    double d=l1.A*l2.B-l2.A*l1.B;
    if(abs(d)<eps) // 不相交
        return false;
    p.x = (l2.C*l1.B-l1.C*l2.B)/d;
    p.y = (l2.A*l1.C-l1.A*l2.C)/d;
    return true;
}


bool intersection(Line_segment l1,Line_segment l2,Point &inter)   // 24.如果线段l1和l2相交，返回true且交点由(inter)返回，否则返回false
{
    Line ll1,ll2;
    ll1=makeline(l1.s,l1.e);
    ll2=makeline(l2.s,l2.e);
    if(lineintersect(ll1,ll2,inter))
        return online(l1,inter) && online(l2,inter);
    else
        return false;
}

/*******************************************************************************/
bool point_in_circle(Point o,double r,Point p)   //25. 返回值:点p在圆内(包括边界)时,返回true
{
    /******************************************************************************
    参数o为圆心,r为半径,p为判断的点
    返回值:点p在圆内(包括边界)时,返回true
    *******************************************************************************/
    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)<eps;
}



bool cocircle(Point p1,Point p2,Point p3,Point &q,double &r) //26.三点确定一个圆,不能构成圆返回false
{
    /******************************************************************************
    用 途 :求不共线的三点确定一个圆
    输 入 :三个点p1,p2,p3
    返回值 :如果三点共线，返回false；反之，返回true。圆心由q返回，半径由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)<eps) //共线，圆不存在
        return false;
    q.x=(y13*z2-y12*z3)/d;
    q.y=(x12*z3-x13*z2)/d;
    r=Dist(p1,q);
    return true;
}

int CircleRelation(Point p1, double r1, Point p2, double r2) //27.两圆位置关系
{
    /******************************************************************************
    相离:return 1
    外切:return 2
    相交:return 3
    内切:return 4
    内含:return 5
    *******************************************************************************/
    double d = sqrt( (p1.x-p2.x)*(p1.x-p2.x)+(p1.y-p2.y)*(p1.y-p2.y) );

    if( fabs(d-r1-r2) < eps ) // 必须保证前两个if先被判定！
        return 2;
    if( fabs(d-fabs(r1-r2)) < eps )
        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!未知错误
}

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

bool SameSide(Point p1, Point p2, Line line) //29.点在直线同侧返回true
{
    return (line.A * p1.x + line.B * p1.y + line.C) *
           (line.A * p2.x + line.B * p2.y + line.C) > 0;
}

void  c2point(Point p1,double r1,Point p2,double r2,Point &rp1,Point &rp2) //30.两个圆(已判断为相交或相切)的交点rp1,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) //31.两相交圆公共面积 +30
{
    double TEMP;
    Point rp1,rp2,rp;
    c2point(p1,r1,p2,r2,rp1,rp2);
    if(r1>r2) //保证r2>r1
    {
        rp=p1;
        p1=p2;
        p2=rp;
        TEMP=r1;
        r1=r2;
        r2=TEMP;
    }
    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;
}

int clpoint(Point p,double r,double a,double b,double c,Point &rp1,Point &rp2) //32.圆和直线(ax+by+c=0,a>=0)关系
{
    /******************************************************************************
    相离 return 0
    相切 return 1
    相交 return 2
    *******************************************************************************/
    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) // 33.三角形内切圆
{
    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;
}

void cutpoint(Point p,double r,Point sp,Point &rp1,Point &rp2) //34.过圆外一点的直线与圆的两个切点(p为圆心,r为圆半径,点sp为圆外一点)
{
    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);
}

int main()
{
    double a,b,c,d;
    Point m(1,1),n(2,2);
    Line_segment w(m,n);
    while(~scanf("%lf %lf",&a,&b))
    {
        Point A(a,b);
        if(online(w,A)) printf("yes\n");
        else printf("n0\n");
    }
    return 0;
}</span>

计算几何模板

标签：

原文地址：http://blog.csdn.net/k_l_c_/article/details/51992298