标签：cnt return Point double ...... 更新 Vector 几何 向量

// 2022/01/22更新

更新内容主要为算法竞赛入门经典——训练指南（升级版）（刘汝佳、陈锋编著）第4章几何问题

1 二维几何基础

在平面坐标系下，向量和点一样，也用两个数x,y表示。第6章介绍齐次坐标的概念，从而在程序上区别开点和向量。以下点和向量都用两个数x,y表示。

不要把点和向量弄混。有如下 点-点=向量，向量+向量=向量，向量-向量=向量，点+向量=点，而点+点是没有意义的。

struct Point{
	double x,y;
	Point(double x=0, double y=0):x(x),y(y) { } // 构造函数，方便代码编写 
};

typedef Point Vector; // 从程序实现上，Vector只是Point的别名

//向量+向量=向量，点+向量=点
Vector operator + (Vector A, Vector 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 * (Vector A, double p) { return Vector(A.x*p, A.y*p); }
//向量/数=向量
Vector operator / (Vector A, double p) { return Vector(A.x/p, A.y/p); }

bool operator < (const Point& a, const Point& b) {
	return a.x<b.x || ( a.x == b.x && a.y < b.y);
} 

const double eps = 1e-10;
int dcmp(double x){
	if(fabs(x) < eps) return 0; else return x < 0 ? -1 : 1;
}

bool operator == (const Point& a, const Point& b){
	return dcmp(a.x-b.x) == 0 && dcmp(a.y-b.y) == 0;
}

注意，上面的“相等”函数用到了“三态函数”dcmp，减少了精度问题。

极角：

向量的极角：从x轴正半轴旋转到该向量方向所需要的角度。
C标准库里的函数就是用来求极角的，如向量(x,y)的极角就是（单位：弧度）。

//极角函数 atan2()函数 的例子
#include<bits/stdc++.h>
using namespace std;
struct Point{
	double x,y;
	Point(double x=0, double y=0):x(x),y(y) { } // 构造函数，方便代码编写 
};
typedef Point Vector; // 从程序实现上，Vector只是Point的别名
int main(){
	ios::sync_with_stdio(false);
	double x,y,Polar_angle; Vector vec; 
	cin>>x>>y;
	vec=Vector(x,y);
	Polar_angle=atan2(y,x);
	cout<<Polar_angle<<endl;
	return 0;
}
/*
#1        #2        #3
0 1       -1 0      13 34
1.5708    3.14159   1.20559
*/

1.1 基本运算

点积

叉积

两个向量的位置关系

向量旋转

基于复数的几何计算

点积：两个向量v和w的点积等于二者长度的乘积再乘上它们夹角的余弦。其中的是指从v到w逆时针旋转的角，因此当夹角大于90度时点积为负。

余弦是偶函数，因此点积满足交换律。如果两向量垂直，点积等于0。不难推导出：在平面坐标系下，两个向量和的点积等于。这里给出点积计算方法，以及利用点积计算向量长度和夹角的函数。代码如下：

struct Point{
	double x,y;
	Point(double x=0, double y=0):x(x),y(y) { } // 构造函数，方便代码编写 
};

typedef Point Vector; // 从程序实现上，Vector只是Point的别名

double Dot(Vector A, Vector B) { return A.x*B.x + A.y*B.y; }
double Length(Vector A) { return sqrt(Dot(A, A)); }
double Angle(Vector A, Vector B) { return acos(Dot(A, B) / Length(A) / Length(B)); }

//2022/01/23

叉积：简单来说，两个向量v和w的叉积等于v和w组成的三角形的有向面积的两倍。不难发现，叉积不满足交换律。事实上，。在坐标系下，两个向量 和的叉积等于。代码如下：

struct Point{
	double x,y;
	Point(double x=0, double y=0):x(x),y(y) { } // 构造函数，方便代码编写 
};

typedef Point Vector; // 从程序实现上，Vector只是Point的别名

//点-点=向量
Vector operator - (Point A, Point B) { return Vector(A.x-B.x, A.y-B.y); }

double Cross(Vector A, Vector B){ return  A.x*B.y-A.y*B.x; }
double Area2(Point A, Point B, Piont C){ return Cross(B-A, C-A); }

两个向量的位置关系：把叉积和点积组合在一起，我们可以更细致地判断两个向量的位置关系。

如图，括号里的第一个数是点积的符号，第二个数是叉积的符号，第一个向量v总是水平向右，另一个向量w的各种情况都包含在了上图中。比如，当w的终点在左上方的第二象限时，点积为负但叉积均为正，用（-，+）表示。

向量旋转： 向量可以绕起点旋转，公式为，其中为逆时针旋转的角。代码如下：

//rad是弧度
Vector Rotate(Vector A, double rad){
	return Vector(A.x*cos(rad)-A.y*sin(rad), A.x*sin(rad)+A.y*cos(rad));
} 

特殊情况：下面函数用来计算向量的单位法线，即左转90度以后把长度归一化。

//调用前请确保A不是零向量
Vector Normal(Vector A){
	double L=Length(A);
	return Vector(-A.y/L, A.x/L);
} 

基于复数的几何计算：使用C++里的复数实现点和向量的方法。如下定义后，我们自动拥有了构造函数、加减法和数量积。用real(p)和imag(p)访问实部和虚部，conj(p)返回共轭复数，即。相关代码如下：

#include<complex>
using namespace std;
typedef complex<double> Point;
typedef Point Vector;

double Dot(Vector A, Vector B){ return real(conj(A)*B); }
double Cross(Vector A, Vector B){ return imag(conj(A)*B); }
Vector Rotate(Vector A, double rad){ return A*exp(Point(0, rad)); }

上述函数的效率不是很高，但是相当方便、好记。

1.2 点和直线

直线的参数表示

直线交点

点到直线的距离

点到线段的距离

点在直线上的投影

线段相交判定

直线的参数表示：

公式：(P0是直线上一点，v是方向向量，t为参数)，如果已知直线上的两个不同点A和B，则方向向量为B-A，所以参数方程为。参数方程最方便的地方在于直线、射线和线段的方程形式是一样的，区别仅仅在于参数t。

类型	t范围
直线	没有范围限制
射线
线段

直线交点：设直线分别为和，设向量，设交点在第一条直线上的参数为，第二条直线上的参数为，则x和y坐标可以各列出一个方程，解得：

， 。代码如下：

//调用前请确保两条直线 P+tv 和 Q+tw 有唯一交点，当且仅当 Cross(v,w)非0
Point GetLineIntersection(Point P, Vector v, Point Q, Vector w){
	Vector u = P-Q;
	double t = Cross(w, u) / Cross(v, w);
	return P+v*t; 
} 

需要注意的是，从上述公式可以看出，如果P,v,Q,w的各个分量均为有理数，则交点坐标也是有理数。在精度要求极高的情况下，可以考虑自定义分数类。（？？？）

点到线段的距离：点到线段的距离有两种可能。设投影点为Q。

①如果Q在线段AB上，则所求距离就是P点到直线AB的距离。

②如果Q不在线段AB上，则分两种情况：

1）如果Q在射线BA上，则所求距离为PA距离；

2）如果Q在射线AB上，则所求距离为PB距离。

判断Q的位置可以用点积进行。代码如下：

double DistanceToSegment(Point P, Point A, Point B){
	if(A == B) return Length(P-A);
	Vector v1 = B - A, v2 = P - A, v3 = P - B;
	if(dcmp(Dot(v1, v2)) < 0) return Length(v2);
	else if(dcmp(Dot(v1, v3)) > 0) return Length(v3);
	else return fabs(Cross(v1, v2)) / Length(v1);
}

点在直线上的投影：为求上述的投影点Q，我们把直线AB写成参数式（v为向量）,并且设Q的参数为，那么。根据PQ垂直于AB，两个向量的点积应该为0，因此。根据分配律有，这样就可以解出，从而得到Q点。代码如下：

Point GetLineProjection(Point P, Point A, Point B){
	Vector v = B-A;
	return A+v*(Dot(v, P-A) / Dot(v, v) ); 
}

线段相交判定：即给定两条线段，判断是否相交。

我们定义“规范相交”为两线段恰好有一个公共点，且公共点不是任何一条线段的端点。（以可以理解成两条线段均不含端点）

线段规范相交的充要条件是：每条线段的两个端点都在另一条线段的两侧（这里的“两侧”是指叉积的符号不同）。代码如下：

bool SegmentProperIntersection(Point a1, Point a2, Point b1, Point b2){
	double c1 = Cross(a2-a1,b1-a1), c2 = Cross(a2-a1,b2-b1);
		   c2 = Cross(b2-b1,a1-b1), c4 = Cross(b2-b1,a2-b1);
	return dcmp(c1)*dcmp(c2)<0 && dcmp(c3)*dcmp(c4)<0;
}

如果允许在端点处相交，情况就比较复杂了：

如果c1和c2都是0，表示两线段共线，这时可能会有部分重叠的情况；

如果c1和c2不都是0，则只有一种相交方法，即某个端点在另外一条线段上。

为了判断上述情况是否发生，还需要判断一个点是否在一条线段上（不包含端点）。代码如下：

bool OnSegment(Point p, Point a1, Point a2){
	return dcmp(Cross(a1-p, a2-p)) == 0 && dcmp(Dot(a1-p, a2-p) < 0);
}

1.3 多边形

计算多边形的有向面积 

①凸多边形，可以从第一个顶点出发把凸多边形分成n-2个三角形，然后把面积加起来。代码如下

double ConvexPolygonArea(Point* p, int n){
	double area = 0;
	for(int i = 1; i < n-1; i++)
		area += Cross(p[i]-p[0], p[i+1]-p[0]);
	return area/2;
}

②非凸多边形，上述方法也对非凸多边形适用。由于三角形面积是有向的，在外面的部分可以正负抵消掉。实际上，可以从任意点出发进行划分。

可以取p[0]点为划分顶点，一方面可以少算两个叉积（0和任意向量的叉积都等于0），另一方面也减少了乘法溢出的可能性，还不用特殊处理（i=n-1的时候，下一个顶点是p[0]而不是p[n]，因为p[n]不存在）。代码如下：

//多边形的有向面积
double PolygonArea(Point* p, int n){
	double area = 0;
	for(int i = 1; i < n-1; i++)
		area += Cross(p[i]-p[0], p[i+1]-p[0]);
	return area/2;
} 

也可以取坐标原点为划分点，乘法次数减少，代码待自行编写。

// 2022/01/25

1.4 例题选讲

1. UVA11178 Morley's Theorem(Morley定理)

#include<bits/stdc++.h>
using namespace std;
struct Point{
	double x,y;
	Point(double x=0, double y=0):x(x),y(y) { }
};
typedef Point Vector;
Vector operator + (Vector A, Vector 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 * (Vector A, double p) { return Vector(A.x*p, A.y*p); }
Vector operator / (Vector A, double p) { return Vector(A.x/p, A.y/p); }
double Dot(Vector A, Vector B) { return A.x*B.x + A.y*B.y; }
double Length(Vector A) { return sqrt(Dot(A, A)); }
double Angle(Vector A, Vector B) { return acos(Dot(A, B) / Length(A) / Length(B)); }
double Cross(Vector A, Vector B){ return  A.x*B.y-A.y*B.x; }
Vector Rotate(Vector A, double rad){
	return Vector(A.x*cos(rad)-A.y*sin(rad), A.x*sin(rad)+A.y*cos(rad));
} 
Point GetLineIntersection(Point P, Vector v, Point Q, Vector w){
	Vector u = P-Q;
	double t = Cross(w, u) / Cross(v, w);
	return P+v*t; 
} 
Point read_point(){
	Point p;
	scanf("%lf %lf",&p.x,&p.y);
	return p;
}
Point getD(Point A, Point B, Point C){
	Vector v1 = C-B;
	double a1 = Angle(A-B, v1);
	v1 = Rotate(v1, a1/3);
	
	Vector v2 = B-C;
	double a2 = Angle(A-C, v2);
	v2 = Rotate(v2, -a2/3); // 负号表示顺时针旋转
	
	return GetLineIntersection(B, v1, C, v2); 
}
int main(){
	int T;
	Point A, B, C, D, E, F;
	scanf("%d",&T);
	while(T--){
		A = read_point();
		B = read_point();
		C = read_point();
		D = getD(A, B, C);
		E = getD(B, C, A);
		F = getD(C, A, B);
		printf("%.6lf %.6lf %.6lf %.6lf %.6lf %.6lf\n", D.x, D.y, E.x, E.y, F.x, F.y); 
	}
	return 0;
} 

2 好看的一笔画 That Nice Euler Circuit，上海 2004，LA 3263

 分析：(还需要思考)

        若是要直接找出所有区域，会非常麻烦，而且容易出错。但用欧拉定理可以将问题进行转化，使解法更容易。

        欧拉定理：设平面图的顶点数、边数和面数分别为V，E和F，则。

        这样，只需求出顶点数V和边数E，就可以求出。

        该平面图的结点由两部分组成，即原来的结点和新增的结点。由于可能出现三线共点，需要删除重复的点。代码如下：

/*********************************************************************************
欧拉定理：设平面图的顶点数、边数和面数分别为V、E和F，则V+F-E=2
so...F=E+2-V;
该平面图的结点有原来的和新增结点构成，由于可能出现三线共点，需要删除重复点
*********************************************************************************/
#include<bits/stdc++.h>
using namespace std;

struct Point{
    double x,y;
    Point(double x=0,double y=0):x(x),y(y){}
};

typedef Point Vector;

//向量+向量=向量;    向量+点=点
Vector operator + (Vector A,Vector 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 * (Vector A,double p){return Vector(A.x*p,A.y*p);}

//向量/数=向量
Vector operator / (Vector A,double p){return Vector(A.x/p,A.y/p);}


bool operator < (const Point& a,const Point& b){return a.x<b.x||(a.x==b.x && a.y<b.y);}


const double eps = 1e-10;

int dcmp(double x){if(fabs(x)<eps)return 0;else return x < 0 ? -1 : 1;}

bool operator == (const Point& a,const Point& b){return dcmp(a.x-b.x)==0 && dcmp(a.y-b.y)==0;}


double Dot(Vector A,Vector B){return A.x*B.x+A.y*B.y;}
double length(Vector A){return sqrt(Dot(A,A));}
double Angle(Vector A,Vector B){return acos(Dot(A,B)/length(A)/length(B));}

double Cross(Vector A,Vector B){return A.x*B.y-B.x*A.y;}
double Area2(Point A,Point B,Point C){return Cross(B-A,C-A);}

/*******两直线交点*******/
//调用前确保两条直线P+tv和Q+tv有唯一交点，当且仅当Cross(v,w)非0;
Point GetLineIntersection(Point P,Vector v,Point Q,Vector w){
    Vector u=P-Q;
    if(Cross(v,w))
    {
        double t=Cross(w,u)/Cross(v,w);//精度高的时候，考虑自定义分数类
        return P+v*t;
    }
//    else
 //       return ;
}

/************************
线段相交判定（规范相交）
************************/
bool SegmentProperIntersection(Point a1,Point a2,Point b1,Point b2){
    double c1=Cross(a2-a1,b1-a1);
    double c2=Cross(a2-a1,b2-a1);
    double c3=Cross(b2-b1,a1-b1);
    double c4=Cross(b2-b1,a2-b1);
    return dcmp(c1)*dcmp(c2)<0&&dcmp(c3)*dcmp(c4)<0;
}
/**如果允许在端点处相交：如果c1和c2都是0，表示共线，如果c1和c2不都是0，则表示某个端点在另一条直线上**/
bool OnSegment(Point p,Point a1,Point a2){
    return dcmp(Cross(a1-p,a2-p))==0&&dcmp(Dot(a1-p,a2-p))<0;
}

const int mmax=310;
Point P[mmax],V[mmax*mmax];

Point read_point(Point &P){
    scanf("%lf%lf",&P.x,&P.y);
    return P;
}

int main(){
    int n,kase=0;
    while(scanf("%d",&n)==1 && n){
        for(int i=0;i<n;i++){
            P[i]=read_point(P[i]);
            V[i]=P[i];
        }
        n--;
        int c=n,e=n;
        for(int i=0;i<n;i++){
            for(int j=i+1;j<n;j++){
                if(SegmentProperIntersection(P[i],P[i+1],P[j],P[j+1])){//严格相交
                    V[c++]=GetLineIntersection(P[i],P[i+1]-P[i],P[j],P[j+1]-P[j]);//交点
                }
            }
        }
 //       printf("c=%d\n",c);
        sort(V,V+c);
        c=unique(V,V+c)-V;
  //      printf("%d=%d-%d\n",c,unique(V,V+c),V);
        for(int i=0;i<c;i++){
            for(int j=0;j<n;j++){
                if(OnSegment(V[i],P[j],P[j+1])) e++;//边数
            }
        }
        printf("Case %d: There are %d pieces.\n",++kase,e+2-c);
    }
    return 0;
}

/*
5
0 0 0 1 1 1 1 0 0 0
7
1 1 1 5 2 1 2 5 5 1 3 5 1 1
7
1 1 1 5 2 1 2 5 5 1 3 9 1 1
0
*/

  

3 UVA11796 狗的距离 Dog Distance

分析：

        先来看一种简单的情况：甲狗和乙狗的路线都是一条线段。因为运动是相对的，因此也可以认为甲狗静止不动，乙狗自己沿着直线走，因此问题转化为求点到线段的最小或最大距离。

        有了简化版的分析，只需模拟整个过程。设现在甲狗的位置在，刚经过编号为的拐点；乙的位置是，刚经过编号为的拐点，则我们只需要计算两狗谁先到拐点，那么在这个时间点之前的问题就是我们刚才讨论过的“简化版”、求解完毕之后，需要更新甲狗和乙狗的位置，如果正好到达下一个拐点，还要更新和，然后继续模拟。因为每次至少有一条狗到达拐点，所以子问题的求解次数不超过A+B。代码如下：

#include<bits/stdc++.h>
using namespace std;
struct Point{
	double x,y;
	Point(double x=0, double y=0):x(x),y(y) { }
};
typedef Point Vector;
Vector operator + (Vector A, Vector 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 * (Vector A, double p) { return Vector(A.x*p, A.y*p); }
Vector operator / (Vector A, double p) { return Vector(A.x/p, A.y/p); }
bool operator < (const Point& a, const Point& b) {
	return a.x<b.x || ( a.x == b.x && a.y < b.y);
} 
const double eps = 1e-10;
int dcmp(double x){
	if(fabs(x) < eps) return 0; else return x < 0 ? -1 : 1;
}
bool operator == (const Point& a, const Point& b){
	return dcmp(a.x-b.x) == 0 && dcmp(a.y-b.y) == 0;
}
double Dot(Vector A, Vector B) { return A.x*B.x + A.y*B.y; }
double Length(Vector A) { return sqrt(Dot(A, A)); }
double Angle(Vector A, Vector B) { return acos(Dot(A, B) / Length(A) / Length(B)); }
double Cross(Vector A, Vector B){ return  A.x*B.y-A.y*B.x; }
Vector Rotate(Vector A, double rad){
	return Vector(A.x*cos(rad)-A.y*sin(rad), A.x*sin(rad)+A.y*cos(rad));
} 
Point GetLineIntersection(Point P, Vector v, Point Q, Vector w){
	Vector u = P-Q;
	double t = Cross(w, u) / Cross(v, w);
	return P+v*t; 
} 
Point read_point(){
	Point p;
	scanf("%lf %lf",&p.x,&p.y);
	return p;
}
double DistanceToSegment(Point P, Point A, Point B){
	if(A == B) return Length(P-A);
	Vector v1 = B - A, v2 = P - A, v3 = P - B;
	if(dcmp(Dot(v1, v2) < 0)) return Length(v2);
	else if(dcmp(Dot(v1, v3)) > 0) return Length(v3);
	else return fabs(Cross(v1, v2)) / Length(v1); 
}
const int maxn = 60;
int T, A, B;
Point P[maxn], Q[maxn];
double Min, Max;
void update(Point P, Point A, Point B){
	Min = min(Min, DistanceToSegment(P, A, B));
	Max = max(Max, Length(P-A));
	Max = max(Max, Length(P-B));
}
int main(){
	scanf("%d",&T);
	for(int kase = 1; kase <= T; kase++){
		scanf("%d%d",&A,&B);
		for(int i=0;i<A;i++) P[i]=read_point();
		for(int i=0;i<B;i++) Q[i]=read_point();
		
		double LenA=0,LenB=0;
		for(int i=0;i<A-1;i++) LenA+=Length(P[i+1]-P[i]);
		for(int i=0;i<B-1;i++) LenB+=Length(Q[i+1]-Q[i]);
		
		int Sa=0,Sb=0;
		Point Pa=P[0],Pb=Q[0];
		Min=1e9,Max=-1e9;
		while(Sa<A-1&&Sb<B-1){
			double La=Length(P[Sa+1]-Pa); //甲到下一拐点的距离 
			double Lb=Length(Q[Sb+1]-Pb); //乙到下一拐点的距离 
			double T=min(La/LenA,Lb/LenB);
			//取合适的单位，可以让甲和乙的速度分别是LenA和LenB 
			Vector Va=(P[Sa+1]-Pa)/La*T*LenA; //甲的位移向量 
			Vector Vb=(Q[Sb+1]-Pb)/Lb*T*LenB; //乙的位移向量 
			update(Pa, Pb, Pb+Vb-Va); //求解“简化版”，更新最小最大距离 
			Pa=Pa+Va;
			Pb=Pb+Vb;
			if(Pa==P[Sa+1]) Sa++;
			if(Pb==Q[Sb+1]) Sb++;
		}
		printf("Case %d: %.0lf\n",kase,Max-Min); 
	}
	return 0;
} 

2 与圆和球有关的计算问题

2.1 圆的相关计算

        圆上任意一点拥有唯一的圆心角，所以在定义圆的时候，可以加一个通过圆心角求坐标的函数。代码如下：

struct Point{
	double x,y;
	Point(double x, double y):x(x), y(y) { }
};
struct Circle {
	Point c; //圆心
	double r; //半径
	Circle(Point c, double r):c(c), r(r) { }
	Point point(double a){
		return Point(c.x + cos(a)*r, c.y + sin(a)*r);
	} 
};

直线和圆的交点： 假设直线为AB，圆的圆心为C，半径为r。

①第一种方法是解方程组。设交点为，代入圆方程后整理得到，进一步整理后得到一元二次方程。根据判别式的值可以分成3种情况，即无交点（相离）、一个交点（相切）和两个交点（相交）。代码如下（变量a,b,c,d,e,f,g对应于上述方程中的字母）：

int getLineCircleIntersection(Line L, Circle C, double& t1, double& t2, Vector<Point>& sol){
	double a = L.v.x, b = L.p.x-C.c.x, c = L.v.y, d = L.p.y-C.c.y;
	double e = a*a+c*c, f = 2*(a*b+c*d), g = b*b+d*d-C.r*C.r;
	double delta = f*f-4*e*g; //判别式
	if(dcmp(delta) < 0) return 0; //相离
	if(dcmp(delta) == 0){ //相切 
		t1=t2= -f / (2*e); sol.push_back(L.point(t1));
		return 1;
	}
	//相交
	t1 = (-f - sqrt(delta)) / (2*e); sol.push_back(L.point(t1)); 
	t2 = (-f + sqrt(delta)) / (2*e); sol.push_back(L.point(t2));
	return 2;
}

        函数返回的是交点的个数，参数sol存放的是交点本身。注意，上述代码并没有清空sol，这给很多题目带来方便：可以反复调用这个函数，把所有交点放在一个sol里。

②第二种方法是几何法。先求圆心C在AB上的投影P，再求向量对应的单位向量，则两个交点分别为和，其中L为P到交点的距离（P与两个交点等距），可以由勾股定理算出。代码略，代补。

两圆相交：假设圆心分别为和，半径分别为和，圆心距为，根据余弦定理可以算出到的角，由向量的极角，加减就可以得到和的极角。有了极角，就可以很方便地计算出和的坐标了（和是两圆交点），代码如下：

// 计算向量极角
double angle(Vector v){ return atan2(v.y, v.x); }
// 两圆相交 
int getCircleCircleIntersection(Circle C1, Circle C2, vector<Point>& sol){
	double d = Length(C1.c - C2.c); //圆心距
	if(dcmp(d) == 0){
		if(dcmp(C1.r - C2.r) == 0) return -1; //两圆重合
		return 0; //两个同心圆 
	} 
	if(dcmp(C1.r + C2.r - d) < 0) return 0; //两圆相离
	if(dcmp(fabs(C1.r-C2.r)-d) > 0) return 0; 
	
	double a = angle(C2.c - C1.c); //向量C1C2的极角
	double da = acos((C1.r*C1.r + d*d - C2.r*C2.r) / (2*C1.r*d));
	//C1C2到C1P1的角
	Point p1 = C1.point(a-da), p2 = C1.point(a+da);
	
	sol.push_back(p1); 
	if(p1 == p2) return 1;
	sol.push_back(p2);
	return 2;
} 

过定点作圆的切线：先求出向量的距离和向量的夹角ang，则向量的极角加减ang就是两条切线的极角，注意切线不存在和只有一条的情况。代码如下：

//过点p到圆C的切线，v[i]是第i条切线的向量，返回切线条数
int getTangents(Point p, Circle C, Vector* v){
	Vector u = C.c - p;
	double dist = Length(u);
	if(dist < C.r) return 0;
	else if(dcmp(dist-C.r) == 0){ //p在圆上，只有一条切线
		v[0] = Rotate(u, PI/2);
		return 1; 
	}else{
		double ang = asin(C.r / dist);
		v[0] = Rotate(u, -ang);
		v[1] = Rotate(u, +ang);
		return 2;
	}
} 

 两圆的公切线：根据两圆的圆心距，从小到大排列一共有6种情况：

①情况一：两圆完全重合，有无数条公切线。

②情况二：两圆内含，没有公共点，没有公切线。

③情况三：两圆内切，有1条外公切线。

④情况四：两圆相交，有2条外公切线。

⑤情况五：两圆外切，有3条公切线（1条内公切线，2条外公切线）。

⑥情况六：两圆外离，有4条公切线（2条内公切线，2条外公切线）。

        可以根据圆心距和半径的关系辨别出这6种情况，然后逐一求解。情况一和情况二没什么需要求解的；情况三和情况五中的内公切线都对应于“过圆上一点求圆的切线”，只需连接圆心和切点，旋转90度后即可知道切线的方向向量。这样，问题的关键是求出情况四、五中的外公切线和情况六中的内外公切线。

        先考虑情况六中的内公切线，它对应于两圆相离的情况。

        根据三角函数定义不难求出角度，然后和前面一样通过极角进行计算即可。        

        外公切线类似。假定，不管两圆是相离、相切还是相交，都是。剩下的过程又和前面一样了。代码如下：

//返回切线的条数，-1表示有无穷多条切线
//a[i]和b[i]分别是第i条切线在圆A和圆B上的切点
int getTangents(Circle A, Circle B, Point* a, Point* b){
	int cnt = 0;
	if(A.r < B.r) { swap(A, B); swap(a, b); }
	int d2 = (A.x-B.x)*(A.x-B.x) + (A.y-B.y)*(A.y-B.y);
	int rdiff = A.r-B.r;
	int rsum = A.r+B.r;
	if(d2 < rdiff*rdiff) return 0; //内含
	
	double base = atan2(B.y-A.y, B.x-A.x);
	if(d2 == 0 && A.r == B.r) return -1; //无限多条切线
	if(d2 == rdiff*rdiff){ //内切，1条切线 
		a[cnt] = A.getPoint(base); b[cnt] = B.getPoint(base); cnt++;
		return 1; 
	} 
	//有外公切线
	double ang = acos((A.r-B.r) / sqrt(d2));
	a[cnt] = A.getPoint(base+ang); b[cnt] = B.getPoint(base+ang); cnt++;
	a[cnt] = A.getPoint(base-ang); b[cnt] = B.getPoint(base-ang); cnt++;
	if(d2 == rsum*rsum){ //一条内公切线
		a[cnt] = A.getPoint(base); b[cnt] = B.getPoint(PI+base); cnt++; 
	}else if(d2 > rsum*rsum){ //两条内公切线
		double ang = acos((A.r+B.r) / sqrt(d2));
		a[cnt] = A.getPoint(base+ang); b[cnt] = B.getPoint(PI+base+ang); cnt++;
		a[cnt] = A.getPoint(base-ang); b[cnt] = B.getPoint(PI+base-ang); cnt++; 
	}
	return cnt;
} 

 2.2 球面相关问题

经纬度转换为空间坐标：公式如下：

代码如下：

//角度转换成弧度
double torad(double deg){
	return deg/180 *acos(-1); //acos(-1)就是PI 
} 

//经纬度（角度）转化为空间坐标
void get_coord(double R, double lat, double lng, double& x, double& y, double& z){
	lat = torad(lat);
	lng = torad(lng);
	x = R*cos(lat)*cos(lng);
	y = R*cos(lat)*sin(lng);
	z = R*sin(lat);
} 

球面距离：问题：已知球面两点，如何求出它们的最短路?注意，只能沿着球面走，不能穿过球的内部。从表面走的话，最近的路径是走圆弧，具体来说是走大圆(Great Circle)圆弧。用一个穿过球心的平面截这个球，截面就是一个大圆。

        怎么求大圆弧长呢？你无须想象那个大圆的空间位置，而可以把它们想象成一个平面问题：求半径为r，弦长为d的圆弧长度。圆心角为，因此弧长为。

2.3 例题选讲

1 UVA12304 2D Geometry 110 in 1! 二维几何110合一！

2 UVA1308 Viva Confetti 圆盘问题

3 二维几何常用算法

点在多边形内的判定

凸包

半平面交

平面区域

标签：cnt,return,Point,double,......,更新,Vector,几何,向量
来源： https://blog.csdn.net/qq_57351574/article/details/122644183

本站声明： 1. iCode9 技术分享网（下文简称本站）提供的所有内容，仅供技术学习、探讨和分享；
2. 关于本站的所有留言、评论、转载及引用，纯属内容发起人的个人观点，与本站观点和立场无关；
3. 关于本站的所有言论和文字，纯属内容发起人的个人观点，与本站观点和立场无关；
4. 本站文章均是网友提供，不完全保证技术分享内容的完整性、准确性、时效性、风险性和版权归属；如您发现该文章侵犯了您的权益，可联系我们第一时间进行删除；
5. 本站为非盈利性的个人网站，所有内容不会用来进行牟利，也不会利用任何形式的广告来间接获益，纯粹是为了广大技术爱好者提供技术内容和技术思想的分享性交流网站。