计算几何模板笔记
计算几何
四元数
纯旋转时，采用四元数比仿射矩阵计算步骤少，思路简单，主要是快。
自己写的代码，没交过题，自以为正确，不保证正确。
注意四元数中旋转角度要除2。
证明跟二维复平面计算方法一模一样，把所有的\(i\)变为\((i,j,k)\)即可，但四维损失交换律，不损失结合律。
表示方式：想要记录一种旋转变换（不包括被旋转点），它直接把过原点的旋转轴\((x,y,z)\)和旋转角度\(\theta\)，方向恒为逆时针，记录下来共4维，就包含使用\(4\times4\)其次矩阵或\(3\times3\)旋转矩阵的全部信息。
记录为\(q=[q_v\ q_s]=[q_x\ q_y\ q_z\ q_s]=[(x,y,z)\sin(\frac\theta 2)\ \cos(\frac\theta2)]\) 为单位四元数，对比二维复平面上表示\((\cos(\theta),\sin(\theta)i)\)
乘法：
\(pq=(p_xi+p_yj+p_zk+p_s)\times(q_xi+q_yj+q_zk+q_s)\)
\(=(p_xi+p_yj+p_zk)\times(q_xi+q_yj+q_zk)+p_v*q_s+q_v*p_s+p_s*q_s\)
\(=-p_xq_x+p_xq_yk-p_xq_zj-p_yq_y+p_yq_zi-p_yq_xk-p_zq_z+p_zq_xj-p_zq_yi+p_v*q_s+q_v*p_s+p_s*q_s\)
\(=-(p_xq_x+p_yq_y+p_zq_z)+(p_xq_yk-p_xq_zj+p_yq_zi-p_yq_xk+p_zq_xj-p_zq_yi)+p_v*q_s+q_v*p_s+p_s*q_s\)
\(=-p_v\cdot q_v+p_v\times q_v+p_v*q_s+q_v*p_s+p_s*q_s\)
\(=[p_v*q_s+q_v*p_s+p_v\times q_v\quad p_sq_s-p_v\cdot q_v]\)
共轭：
\(\theta_2=-\theta_1\)
\(\therefore q*=[q_v\frac{\sin(-\theta)}{\sin(\theta)}\ q_s\frac{\cos(-\theta)}{\cos(\theta)}]=[-q_v\ q_s]\)
逆：
单位四元数：\(I=[0\ 0\ 0\ 1]\)，因为乘它后位置不变即旋转\(\theta=0\)。
对比二维复平面乘法意义为旋转+缩放，\(q\)为旋转\(\theta\)+缩放\(|q|\)，\(qq^*\)则为旋转\(\theta+(-\theta)=0\)，缩放为\(|q|^2\)。为防止每次缩放后再乘\(\frac1{|q|^2}\)恢复模长，在计算前就将所有旋转变换规范化为单位四元数而其物理意义不变，之后再作多次乘法模长不变。
\(\because qq^-=I,\ qq^*=|q|^2=|q|^2I\quad I=\frac{qq^*}{|q|^2}\)
\(\therefore q^-=\frac{q^*}{|q|^2}\)
使用时先规范化，令\(|q|=1\)，\(q^-=q^*\)
或按上面推出的乘法公式计算结果也相同：
\(qq^*=[q_v\ q_s][-q_v\ q_s]=[q_v*q_s-q_v*q_s+q_v\times(-q_v)\quad q_sq_s-q_v\cdot(-q_v)]\)
\(=[0\ 1]=I\)
\(\therefore q^*=q^-\)
#define db double
#define eps 1e-10
#define pi acos(-1)
struct vec{//点(这里更好理解成从原点出发的向量)-三维
db x, y, z;
vec(){}
vec(db _x, db _y, db _z){
x=_x, y=_y, z=_z;
}
vec operator ^(const vec &n)const{
return vec(y*n.z-z*n.y, z*n.x-x*n.z, x*n.y-y*n.x);
}
db operator *(const vec &n)const{
return x*n.x+y*n.y+z*n.z;
}
vec operator *(const db t)const{
return vec(x*t, y*t, z*t);
}
vec operator +(const vec &n)const{
return vec(x+n.x, y+n.y, z+n.z);
}
};
struct sy{//四元素-包含旋转轴(从原点出发的向量)+旋转角度(方向只能向右手定则转)
vec v;
db s;
sy(){}
sy(vec _v, db th, bool it){
th/=2;
s=cos(th);
db t=sin(th), len=_v.x*_v.x+_v.y*_v.y+_v.z*_v.z;
len=sqrt(len);
v.x=_v.x/len*t, v.y=_v.y/len*t, v.z=_v.z/len*t;
}
sy(vec _v, db _s){
v=_v, s=_s;
}
sy(vec _v){
v=_v, s=0;
}
sy operator ^(const sy &n)const{
return sy(n.v*s+v*n.s+(v^n.v), s*n.s-v*n.v);
}
sy inv(){
return sy(v*-1.0, s);
}
sy rot(vec _v){
sy tv=sy(_v);
return (*this)^tv^(*this).inv();
}
void p(){
printf("%.6f %.6f %.6f\n%.6f\n\n", v.x, v.y, v.z, s);
}
};
db x, y, z, px, py, pz, th;
int main() {
th=pi/2;
while(cin>>x>>y>>z>>px>>py>>pz){
sy n=sy(vec(x, y, z), th, 1);
sy p=n.rot(vec(px, py, pz));
p.p();
}
return 0;
}
凸包
二维凸包
凸多边形
所有内角大小都在\([0,\pi]\)范围内的 简单多边形。
凸包
在平面上能包含所有给定点的最小凸多边形叫做凸包。
对于给定集合\(X\)，所有包含\(X\)的凸集的交集S被称为\(X\)的 凸包 。
凸包用最小的周长围住了给定的所有点。如果一个凹多边形围住了所有的点，它的周长一定不是最小。根据三角不等式，凸多边形在周长上一定是最优的。
凸包的求法
常用的求法有 Graham 扫描法和 Andrew 算法。
Andrew
把所有点以横坐标为第一关键字，纵坐标为第二关键字排序。
排序后最小的元素和最大的元素一定在凸包上。因为是凸多边形，从一个点出发逆时针走，轨迹总是“左拐”的，一旦出现右拐，就说明这一段不在凸包上。用一个单调栈来维护上下凸壳。
注：判断两向量\(\overrightarrow a\)逆时针方向到\(\overrightarrow b\)转过角度是否大于\(\pi\)，用叉乘而不是点乘。
\(\mid \overrightarrow a\times\overrightarrow b\mid=a\times b\times\sin\theta=\begin{vmatrix}i&j&k\\a_x&a_y&a_z\\b_x&b_y&b_z\end{vmatrix}=(a_yb_z-a_zb_y,a_zb_x-a_xb_z,a_xb_y-a_yb_x)\)
本题中方向\(j=a_xb_y-a_yb_x\)。
ll andrew(){
int stk[maxn], tp=0;//为了好取2个
bool used[maxn]={};//为了重复计算减少重复计算下凸壳,不加也行就慢一点
tp=0;
sort(p+1, p+1+n);
for(int i=1; i<=n; i++){
while(tp>=2 && (p[stk[tp]]-p[stk[tp-1]])*(p[i]-p[stk[tp]])<=0)//*是重载的x运算符，-也是
used[stk[tp--]]=0;
used[i]=1, stk[++tp]=i;
}
int tmp=tp;
used[1]=0; //使得1在最后封闭凸包时也对单调栈更新
for (int i=n-1; i>0; i--)
if(!used[i]){
while(tp>tmp &&
(p[stk[tp]]-p[stk[tp-1]])*(p[i]-p[stk[tp]])<=0)
used[stk[tp--]]=0;
used[i]=1, stk[++tp]=i;
}
for(int i=1; i<=tp-1; i++)
h[i]=p[stk[i]];
return tp-1;
}
Graham 扫描
除入队顺序外，与 Andrew 相同。
入队顺序：找最左下角的点s。（可使以s为基点的\(\theta\)在\(0-\frac{\pi}{2}\)中，从而得到\(sin\theta\)排序即为\(\theta\)排序，虽然在这里没用）。以s为基点求出其他所有点的极角，按照极角从小到大排序。（不必计算，重载运算符即可，当\(\overrightarrow a\times\overrightarrow b\)大于零时，\(\overrightarrow a\)的极角小于\(\overrightarrow b\)的极角）
平面几何
重复用sin，cos 500次要预先存好，而不是每次计算，会t。
基本结构
int sgn(double x)
{
if(fabs(x) < eps)return 0;
if(x < 0) return -1;
else return 1;
}
struct Point
{
double x,y;
Point(){}
Point(double _x,double _y)
{
x = _x;y = _y;
}
Point operator -(const Point &b)const
{
return Point(x - b.x,y - b.y);
}
//叉乘
double operator ^(const Point &b)const
{
return x*b.y - y*b.x;
}
//点乘
double operator *(const Point &b)const
{
return x*b.x + y*b.y;
}
};
struct Line
{
Point s,e;
double k;//极角
Line(){}
Line(Point _s,Point _e)
{
s = _s;e = _e;
k = atan2(e.y-s.y,e.x-s.x);
}
Point operator &(const Line &b)const
{ //求两直线交点
Point res=s;
double t=((s-b.s)^(b.s-b.e))/((s-e)^(b.s-b.e));
//double t=((b.s-s)^(b.e-b-s))/((e-s)^(b.e-b.s));等价
res.x+=(e.x-s.x)*t;
res.y+=(e.y-s.y)*t;
return res;
}
//lxl2
double cc(Line l2){
return ((e-s)^(l2.e-l2.s));
}
//线段长
double Len(){
return sqrt((e.x-s.x)*(e.x-s.x)
+(e.y-s.y)*(e.y-s.y));
}
//点到直线距离
double PtoL(Point p){
if(Len()<=0) return Line(s, p).Len();
return fabs((p-s)^(e-s))/Len();
}
};
直线相交
void f(Line l1, Line l2) {
if(sgn((l1.e-l1.s)^(l2.e-l2.s))==0) {
if(sgn((l1.e-l2.s)^(l2.e-l2.s))==0); //重叠
else; //平行
} else {
double t=((l1.s-l2.s)^(l2.e-l2.s))
/(-((l1.e-l1.s)^(l2.e-l2.s)));
double tx=l1.s.x+(l1.e.x-l1.s.x)*t;
double ty=l1.s.y+(l1.e.y-l1.s.y)*t;
//交于(tx, ty)
}
}
线段交
快速排斥实验+跨立实验
当两条线段共线但不相交时也可以通过跨立实验，需要与快速排斥实验结合。
bool inter(Line l1,Line l2)
{
return
max(l1.s.x,l1.e.x) >= min(l2.s.x,l2.e.x) &&
max(l2.s.x,l2.e.x) >= min(l1.s.x,l1.e.x) &&
max(l1.s.y,l1.e.y) >= min(l2.s.y,l2.e.y) &&
max(l2.s.y,l2.e.y) >= min(l1.s.y,l1.e.y) &&
sgn((l2.s-l1.s)^(l1.e-l1.s))*sgn((l2.e-l1.s)^(l1.e-l1.s)) <= 0 &&
sgn((l1.s-l2.s)^(l2.e-l2.s))*sgn((l1.e-l2.s)^(l2.e-l2.s)) <= 0;
}
多边形
多边形表示方法，逆时针方向的的点集。
PIP问题
判断一点是否在任意多边形内部。
光线投射算法
判断矩阵+找射线焦点，射线斜率取无理数，避免与边重合。奇个交点在内部。
回转数算法
回转数 ： 平面内闭合曲线逆时针绕过该点的总次数 。
非零规则 ： 当回转数等于0，点在曲线外部 。
把该点与多边形的所有顶点连接起来，计算相邻两边夹角的和 。如果夹角和为0，则点在外部。
多边形面积
找一个随便找个参考点，叠加三角形面积，面积用\(\frac{1}{2}\)叉乘求。
逆时针存储多边形的各点，从0开始。
double calS(){
double ans=0;
for(int i=0; i<n; i++)
ans+=(p[i]^p[(i+1)%n]);
return ans/2;
}
平面交
Line Q[maxn];//辅助
Point p[maxn]; //记录最初给的点集，读入时用，从0开始，main中
Line line[maxn]; //由最初的点集生成直线的集合，辅助，我不用手动写
Point pp[maxn]; //记录半平面交的结果的点集
bool HPIcmp(Line a,Line b)
{
if(fabs(a.k-b.k)>eps)return a.k<b.k;
return ((a.s-b.s)^(b.e-b.s))<0;
}
void HPI(Line line[], int n, Point res[], int &resn)
{
//res是结果
int tot=n;
sort(line,line+n,HPIcmp);
tot=1;
for(int i=1;i<n;i++)
if(fabs(line[i].k-line[i-1].k)>eps) //去掉斜率重复的
line[tot++]=line[i];
int head=0, tail=1;
Q[0]=line[0];
Q[1]=line[1];
resn=0;
for(int i=2; i<tot; i++)
{
if(fabs((Q[tail].e-Q[tail].s)^(Q[tail-1].e-Q[tail-1].s))<eps ||
fabs((Q[head].e-Q[head].s)^(Q[head+1].e-Q[head+1].s))<eps)
return;
while(head<tail && (((Q[tail]&Q[tail-1])-line[i].s)^
(line[i].e-line[i].s))>eps) tail--;
while(head<tail && (((Q[head]&Q[head+1])-line[i].s)^
(line[i].e-line[i].s))>eps) head++;
Q[++tail]=line[i];
}
while(head<tail && (((Q[tail]&Q[tail-1])-Q[head].s)^
(Q[head].e-Q[head].s))>eps) tail--;
while(head<tail && (((Q[head]&Q[head-1])-Q[tail].s)^
(Q[tail].e-Q[tail].e))>eps) head++;
if(tail<=head+1) return;
for(int i=head; i<tail; i++)
res[resn++]=Q[i]&Q[i+1];
if(head<tail-1)
res[resn++]=Q[head]&Q[tail];
return;
}
void change(Point a, Point b, Point &c, Point &d, double p)
{ //将线段ab往左移动距离p,修改得到线段cd
double len=a.Len(b);
/*三角形相似推出下面公式*/
double dx=(a.y-b.y)*p/len;
double dy=(b.x-a.x)*p/len;
c.x=a.x+dx; c.y=a.y+dy;
d.x=b.x+dx; d.y=b.y+dy;
}
//调用
for(int i=0; i<n; i++){
Point t1,t2;
change(p[i], p[(i+1)%n], t1, t2, d);
line[i]=Line(t1, t2);
}
int si;
HPI(line, n, pp, si);
//最终形成凸多边形的点逆时针存入pp，大小为si，从0开始
kuangbin模板
int sgn(double x) {
if(fabs(x) < eps)return 0;
if(x < 0) return -1;
else return 1;
}
db hypot(db x, db y){
return sqrt(x*x+y*y);
}
点
struct Point {
db x, y;
Point() {}
Point(db x, db y):x(x), y(y) {}
void input() {
scanf("%lf%lf",&x, &y);
}
// 1.比较两点是否相同
bool operator == (Point b)const {
return sgn(x-b.x) == 0 && sgn(y-b.y) == 0;
}
// 2.按照(x,y)的优先级比较，因题目不同要修改
bool operator < (Point b)const {
return sgn(x-b.x) == 0 ? sgn(y - b.y) < 0 : x < b.x;
}
Point operator - (const Point &b)const {
return Point(x - b.x, y - b.y);
}
db operator ^ (const Point &b)const {
return x * b.y - y * b.x;
}
db operator * (const Point &b)const {
return x * b.x + y * b.y;
}
db len() {
return hypot(x, y);
}
db len2() {
return x * x + y * y;
}
// 8. 返回两点距离
db distance(Point p) {
return hypot(x - p.x, y - p.y);
}
// 9. 向量加
Point operator + (const Point &b)const {
return Point(x + b.x, y + b.y);
}
// 10. 标量乘
Point operator * (const db &k) const {
return Point(x * k, y * k);
}
// 11. 标量除
Point operator /(const db &k)const {
return Point(x / k, y / k);
}
// 12. 计算 pa 和 pb 的夹角
db rad(Point a, Point b) {
Point p = *this;
return fabs(atan2( fabs((a-p)^(b-p)), (a-p) * (b-p) ));
}
// 13. 化为长度为 r 的向量
Point trunc(db r) {
db l = len();
if(!sgn(l)) return *this;
r /= l;
return Point(x*r, y*r);
}
// 14. 逆时针旋转90度
Point rotleft() {
return Point(-y, x);
}
// 15. 顺时针转90度
Point rotright() {
return Point(y, -x);
}
// 16. 绕点 p 逆时针旋转 angle
Point rotate(Point p, db angle) {
Point v = (*this) - p;
db c = cos(angle), s = sin(angle);
return Point(p.x + v.x * c - v.y * s, p.y + v.x * s + v.y * c);
}
};
线
struct Line {
Point s, e;
Line() {}
Line(Point s, Point e):s(s),e(e) {}
void input() {
s.input();
e.input();
}
// 17. 判断有向线段是否相等
bool operator == (Line v) {
return (s == v.s) && (e == v.e);
}
// 20. 调整直线两点顺序（按Point重载的<调整）
void adjust() {
if(e < s) swap(s, e);
}
// 21. 求直线长度
db length() {
return s.distance(e);
}
// 18. 根据一个点和倾斜角angle确定直线， 0 <= angle < pi
Line (Point p, db angle) {
s = p;
if(sgn(angle - pi / 2) == 0) {
e = (s + Point(0, 1));
} else {
e = (s + Point(1, tan(angle)));
}
}
// 19. ax + by + c = 0
Line (db a, db b, db c) {
if(sgn(a) == 0) { // y = -c / b
s = Point(0, -c/b);
e = Point(1, -c/b);
} else if (sgn(b) == 0) { // x = -c / a
s = Point(-c/a, 0);
e = Point(-c/a, 1);
} else { //(0, -c/b), (1, (-c-a)/b)
s = Point(0, -c/b);
e = Point(1, (-c-a)/b);
}
}
// 22. 返回直线倾斜角 0 <= angle < pi
db angle() {
db k = atan2(e.y - s.y, e.x - s.x);
if(sgn(k) < 0) k += pi;
if(sgn(k - pi) == 0) k -= pi;
return k;
}
/*
23.
点和直线的关系
1 在左侧
2 在右侧
3 在直线上
*/
int relation(Point p) {
int c = sgn((p-s) ^ (e-s));
if(c < 0) return 1;
else if(c > 0) return 2;
else return 3;
}
// 24. 点在线段上的判断，包括端点 第二个判断改为小于表示不包括端点
bool pointonseg(Point p) {
return sgn((p-s)^(e-s)) == 0 && sgn((p-s) * (p-e)) <= 0;
}
// 25. 两向量平行（对应直线平行或重合）
bool parallel(Line v) {
return sgn((e - s) ^ (v.e - v.s)) == 0;
}
/*
26.
两线段相交判断
2 规范相交
1 非规范相交
0 不相交
*/
int segcrosseg(Line v) {
int d1 = sgn((e - s) ^ (v.s - s)); //v.s 在 线段的哪一侧
int d2 = sgn((e - s) ^ (v.e - s)); //v.e 在 线段的哪一侧
int d3 = sgn((v.e - v.s) ^ (s - v.s));
int d4 = sgn((v.e - v.s) ^ (e - v.s));
if((d1^d2) == -2 && (d3^d4) == -2) return 2; // 跨立实验满足 一个是-1一个是1
// 1. v.s在线段上 || v.e 在线段上 || s 在另外一条线段上 || e在另外一条线段上
return (d1 == 0 && sgn((v.s - s) * (v.s - e)) <= 0) ||
(d2 == 0 && sgn((v.e - s) * (v.e - e)) <= 0) ||
(d3 == 0 && sgn((s - v.s) * (s - v.e)) <= 0) ||
(d4 == 0 && sgn((e - v.s) * (e - v.e)) <= 0);
}
/*
27.
直线与线段相交判断
*this line -v seg
2 规范相交
1 非规范相交
0 不相交
*/
int linecorssseg(Line v) { // v是线段
int d1 = sgn((e - s) ^ (v.s - s));
int d2 = sgn((e - s) ^ (v.e - s));
if((d1 ^ d2) == -2) return 2;
return (d1 == 0 || d2 == 0);
}
/*
28.
两直线关系
0 平行
1 重合
2 相交
*/
int linecrossline(Line v) {
// 如果平行，则看点是否在直线上
if((*this).parallel(v)) return v.relation(s) == 3;
return 2;
}
/*
29.
求两直线的交点
要保证两直线不平行或重合
*/
Point crosspoint(Line v) {
db a1 = (v.e - v.s) ^ (s - v.s);
db a2 = (v.e - v.s) ^ (e - v.s);
return Point((s.x * a2 - e.x * a1) / (a2 - a1), (s.y * a2 - e.y * a1) / (a2 - a1));
}
// 30. 点到直线的距离
db dispointtoline(Point p) {
return fabs((p-s) ^ (e-s)) / length();
}
// 31. 点到线段的距离
db dispointtoseg(Point p) {
if(sgn((p-s)*(e-s)) < 0 || sgn((p-e) * (s-e)) < 0)
return min(p.distance(s), p.distance(e));
return dispointtoline(p);
}
/*
32.
返回线段到线段的距离
前提是两线段不相交，相交距离就是 0 了
*/
db dissegtoseg(Line v) {
return min(min(dispointtoseg(v.s), dispointtoseg(v.e)), min(v.dispointtoseg(s), v.dispointtoseg(e)));
}
/*
33. 返回 p 在直线上的投影
*/
Point lineprog(Point p) {
return s + ( ((e-s)*((e-s)*(p-s)))/((e-s).len2()) );
}
// 34. 返回点 p 关于直线的对称点
Point symmetrypoint(Point p) {
Point q = lineprog(p);
return Point(2*q.x - p.x, 2 * q.y - p.y);
}
};
圆
struct circle {
Point p;
db r;
circle() {}
circle(Point p, db r):p(p), r(r) {}
void input() {
p.input();
scanf("%lf", &r);
}
// 通过圆心角确定圆上的一个点
Point point(double a) {
return Point(p.x + cos(a) * r, p.y + sin(a) * r);
}
bool operator == (circle v) {
return (p == v.p) && sgn(r - v.r) == 0;
}
bool operator < (circle v)const {
return (((p<v.p) || (p == v.p)) && sgn(r - v.r) < 0);
}
// 面积
db area() {
return pi * r * r;
}
// 周长
db circumference() {
return 2 * pi * r;
}
/*
三角形的外接圆
需要Point 的 + / rotate() 以及 Line 的crosspoint()
利用两条边的中垂线得到圆心
*/
circle(Point a, Point b, Point c) {
Line u = Line((a+b)/2,((a+b)/2)+((b-a).rotleft()));
Line v = Line((b+c)/2,((b+c)/2)+((c-b).rotleft()));
p = u.crosspoint(v);
r = p.distance(a);
}
/*
三角形的内切圆
bool t 没有作用，只是为了和上面外接圆函数区别
*/
circle(Point a, Point b, Point c, bool t) {
Line u, v;
// u 为角 a 的平分线
db m = atan2(b.y-a.y, b.x-a.x), n = atan2(c.y - a.y, c.x - a.x);
u.s = a;
u.e = u.s + Point(cos((n+m)/2), sin((n+m)/2));
// v 为角 b 的平分线
m = atan2(a.y-b.y, a.x-b.x), n = atan2(c.y-b.y, c.x-b.x);
v.s = b;
v.e = v.s + Point(cos((n+m)/2), sin((n+m)/2));
p = u.crosspoint(v);
r = Line(a,b).dispointtoseg(p);
}
/*
点和圆的关系
0 圆外
1 圆上
2 圆内
*/
int relation(Point b) {
db dst = b.distance(p);
if(sgn(dst - r) < 0) return 2;
else if(sgn(dst - r) == 0) return 1;
return 0;
}
/*
线段和圆的关系
比较的是圆心到线段的距离和半径的关系
2 交
1 切
0 不交
*/
int relationseg(Line v) {
db dst = v.dispointtoseg(p);
if(sgn(dst - r) < 0) return 2;
else if(sgn(dst - r) == 0) return 1;
return 0;
}
int relationline(Line v) {
db dst = v.dispointtoline(p);
if(sgn(dst - r) < 0) return 2;
else if(sgn(dst - r) == 0) return 1;
return 0;
}
/*
两圆的关系
5 相离
4 外切
3 相交
2 内切
1 内含
*/
int relationcircle(circle v) {
db d = p.distance(v.p);
if(sgn(d - r - v.r) > 0) return 5;
if(sgn(d - r - v.r) == 0) return 4;
db l = fabs(r - v.r);
if(sgn(d - r - v.r) < 0 && sgn(d - l) > 0) return 3;
if(sgn(d - l) == 0) return 2;
if(sgn(d - l) < 0) return 1;
}
/*
求两个圆的交点，返回0表示没有交点，返回1是一个交点，2是两个交点
*/
int pointcrosscircle(circle v, Point &p1, Point &p2) {
int rel = relationcircle(v);
if(rel == 1 || rel == 5) return 0;
db d = p.distance(v.p);
db l = (d * d + r * r - v.r * v.r) / (2 * d);
db h = sqrt(r * r - l * l);
Point tmp = p + (v.p - p).trunc(l);
p1 = tmp + ((v.p - p).rotleft().trunc(h));
p2 = tmp + ((v.p - p).rotright().trunc(h));
if(rel == 2 || rel == 4)return 1;
return 2;
}
// 求直线与圆的交点，返回交点个数
int pointcrossline(Line v, Point &p1, Point &p2) {
if(!(*this).relationline(v)) return 0;
Point a = v.lineprog(p);
db d = v.dispointtoline(p);
d = sqrt(r * r - d * d);
if(sgn(d) == 0) {
p1 = a;
p2 = a;
return 1;
}
p1 = a + (v.e - v.s).trunc(d);
p2 = a - (v.e - v.s).trunc(d);
return 2;
}
// 得到 通过a，b两点，半径为r1的两个圆
int getcircle(Point a, Point b, db r1, circle &c1, circle &c2) {
circle x(a, r1), y(b, r1);
int t = x.pointcrosscircle(y, c1.p, c2.p);
if(!t) return 0;
c1.r = c2.r = r1;
return t;
}
// 得到与直线 u 相切，过点 q， 半径为 r1 的圆
int getcircle(Line u, Point q, db r1, circle &c1, circle &c2) {
db dis = u.dispointtoline(q);
if(sgn(dis - r1 * 2) > 0) return 0;
if(sgn(dis) == 0) {
c1.p = q + ((u.e - u.s).rotleft().trunc(r1));
c2.p = q + ((u.e - u.s).rotright().trunc(r1));
c1.r = c2.r = r1;
return 2;
}
Line u1 = Line((u.s + (u.e - u.s).rotleft().trunc(r1)), (u.e + (u.e - u.s).rotleft().trunc(r1)));
Line u2 = Line((u.s + (u.e - u.s).rotright().trunc(r1)), (u.e + (u.e - u.s).rotright().trunc(r1)));
circle cc = circle(q, r1);
Point p1, p2;
if(!cc.pointcrossline(u1, p1, p2)) cc.pointcrossline(u2, p1, p2);
c1 = circle(p1, r1);
if(p1 == p2) {
c2 = c1;
return 1;
}
c2 = circle(p2, r1);
return 2;
}
// 同时与直线u，v相切，半径为r1的圆 , u,v不平行
int getcircle(Line u, Line v, db r1, circle &c1, circle &c2, circle &c3, circle &c4) {
if(u.parallel(v)) return 0;
Line u1 = Line(u.s + (u.e - u.s).rotleft().trunc(r1), u.e + (u.e - u.s).rotleft().trunc(r1));
Line u2 = Line(u.s + (u.e - u.s).rotright().trunc(r1), u.e + (u.e - u.s).rotright().trunc(r1));
Line v1 = Line(v.s + (v.e - v.s).rotleft().trunc(r1), v.e + (v.e - v.s).rotright().trunc(r1));
Line v2 = Line(v.s + (v.e - v.s).rotright().trunc(r1), v.e + (v.e - v.s).rotright().trunc(r1));
c1.r = c2.r = c3.r = c4.r = r1;
c1.p = u1.crosspoint(v1);
c2.p = u1.crosspoint(v2);
c3.p = u2.crosspoint(v1);
c4.p = u2.crosspoint(v2);
return 4;
}
// 同时与不相交圆 cx, cy 相切，半径为r1的圆
int getcircle(circle cx, circle cy, db r1, circle &c1, circle &c2) {
circle x(cx.p, r1+cx.r), y(cy.p, r1+cy.r);
int t = x.pointcrosscircle(y, c1.p, c2.p);
if(!t) return 0;
c1.r = c2.r = r1;
return t;
}
// 过一点作圆的切线 (先判断点和圆的关系)
int tangentline(Point q, Line &u, Line &v) {
int x = relation(q);
if(x == 2) return 0; //圆内
if(x == 1) { //圆上
u = Line(q, q+(q-p).rotleft());
v = u;
return 1;
}
db d = p.distance(q);
db l = r * r / d;
db h = sqrt(r * r - l * l);
u = Line(q, p + ((q - p).trunc(l) + (q-p).rotleft().trunc(h)));
v = Line(q, p + (q - p).trunc(l) + (q - p).rotright().trunc(h));
return 2;
}
// 求两圆相交的面积
db areacircle(circle v) {
int rel = relationcircle(v);
if(rel >= 4) return 0.0;
if(rel <= 2) return min(area(), v.area()); //内部
db d = p.distance(v.p);
db hf = (r + v.r + d) / 2.0;
db ss = 2 * sqrt(hf*(hf - r) * (hf - v.r) * (hf - d));
db a1 = acos((r * r + d * d - v.r * v.r) / (2.0 * r * d));
a1 = a1 * r * r;
db a2 = acos((v.r * v.r + d * d - r * r) / (2.0 * v.r * d));
a2 = a2 * v.r * v.r;
return a1 + a2 - ss;
}
// 求圆和三角形 pab 的相交面积
// POJ3675 HDU3982 HDU2892
db areatriangle(Point a, Point b) {
if(sgn((p-a)^(p-b)) == 0)return 0.0;
Point q[5];
int len = 0;
q[len++] = a;
Line l(a, b);
Point p1, p2;
if(pointcrossline(l, q[1], q[2]) == 2) {
if(sgn((a-q[1]) * (b-q[1])) < 0) q[len++] = q[1];
if(sgn((a-q[2]) * (b-q[2])) < 0) q[len++] = q[2];
}
q[len++] = b;
if(len == 4 && sgn((q[0] - q[1]) * (q[2] - q[1])) > 0) swap(q[1], q[2]);
db res = 0;
for(int i=0; i<len-1; i++) {
if(relation(q[i]) == 0 || relation(q[i+1]) == 0) {
db arg = p.rad(q[i], q[i+1]);
res += r * r * arg / 2.0;
} else {
res += fabs((q[i] - p) ^ (q[i+1] - p)) / 2.0;
}
}
return res;
}
// a[i] 存放第 i 条公切线与 圆A 的交点
int getTangents(circle A, circle B, Point*a, Point *b) {
int cnt = 0;
// 以A为半径更大的那个圆进行计算
if(A.r < B.r) return getTangents(B, A, b, a);
db d2 = (A.p-B.p).len2(); // 圆心距平方
db rdiff = A.r - B.r; // 半径差
db rsum = A.r + B.r; //半径和
if(d2 < rdiff * rdiff) return 0; // 情况1，内含,没有公切线
Point AB = B.p - A.p; // 向量AB，其模对应圆心距
db base = atan2(AB.y, AB.x); // 求出向量AB对应的极角
if(d2 == 0 && A.r == B.r) return -1;// 情况3，两个圆重合，无限多切线
if(d2 == rdiff * rdiff) { // 情况2，内切，有一条公切线
a[cnt] = A.point(base);
b[cnt] = B.point(base);
cnt++;
return 1;
}
// 求外公切线
db ang = acos((A.r - B.r) / sqrt(d2)); //求阿尔法
// 两条外公切线
a[cnt] = A.point(base+ang);
b[cnt] = B.point(base+ang);
cnt++;
a[cnt] = A.point(base-ang);
b[cnt] = B.point(base-ang);
cnt++;
if(d2 == rsum * rsum) { // 情况5，外切，if里面求出内公切线
a[cnt] = A.point(base);
b[cnt] = B.point(pi+base);
cnt++;
} else if(d2 > rsum * rsum) { //情况6，相离，再求出内公切线
db ang = acos((A.r + B.r) / sqrt(d2));
a[cnt] = A.point(base + ang);
b[cnt] = B.point(pi+base+ang);
cnt++;
a[cnt] = A.point(base - ang);
b[cnt] = B.point(pi+base-ang);
cnt++;
}
// 此时，d2 < rsum * rsum 代表情况 4 只有两条外公切线
return cnt;
}
};
多圆交面积
const int N=1013;
struct circles {
circle c[N];
double ans[N];
double pre[N];
int n;
circles() {}
void add(circle cc) {
c[n++] = cc;
}
// x 包含在 y 中
bool inner(circle x, circle y) {
if(x.relationcircle(y) != 1) return 0;
return sgn(x.r - y.r) <= 0 ? 1 : 0;
}
double areaarc(double th, double r) {
return 0.5 * r * r * (th - sin(th));
}
// 圆的面积并，去掉内含的圆
void init_or() {
int i, j, k = 0;
bool mark[N] = {0};
for(i = 0; i < n; i++) {
for(j = 0; j < n; j ++) {
if(i != j && !mark[j]) {
if((c[i] == c[j]) || inner(c[i], c[j])) break;
}
}
if(j < n) mark[i] = 1;
}
for(i = 0; i < n; i++) {
if(!mark[i]) {
c[k++] = c[i];
}
}
n = k;
}
// 圆的面积交，去掉内含的圆
void init_and() {
int i, j, k = 0;
bool mark[N] = {0};
for(i = 0; i < n; i++) {
for(j = 0; j < n; j ++) {
if(i != j && !mark[j]) {
if((c[i] == c[j]) || inner(c[j], c[i])) break;
}
}
if(j < n) mark[i] = 1;
}
for(i = 0; i < n; i++) {
if(!mark[i]) {
c[k++] = c[i];
}
}
n = k;
}
void getarea() {
memset(ans, 0, sizeof ans);
vector<pair<double,int> > v;
for(int i=0; i<n; i++) {
v.clear();
v.push_back(make_pair(-pi, 1));
v.push_back(make_pair(pi, -1));
for(int j=0; j<n; j++) {
if(i != j) {
Point q = (c[j].p - c[i].p);
db ab = q.len(), ac = c[i].r, bc = c[j].r;
if(sgn(ab + ac - bc) <= 0) {
v.push_back(make_pair(-pi, 1));
v.push_back(make_pair(pi, -1));
continue;
}
if(sgn(ab + bc - ac) <= 0) continue;
if(sgn(ab - ac - bc) > 0) continue;
double th = atan2(q.y, q.x), fai = acos((ac*ac + ab * ab - bc * bc) / (2.0*ac*ab));
double a0 = th - fai;
if(sgn(a0 + pi) < 0) a0 += 2 * pi;
db a1 = th + fai;
if(sgn(a1 - pi) > 0) a1 -= 2 * pi;
if(sgn(a0 - a1) > 0) {
v.push_back(make_pair(a0, 1));
v.push_back(make_pair(pi, -1));
v.push_back(make_pair(-pi, 1));
v.push_back(make_pair(a1, -1));
} else {
v.push_back(make_pair(a0, 1));
v.push_back(make_pair(a1, -1));
}
}
}
sort(v.begin(),v.end());
int cur = 0;
for(int j=0; j<(int)v.size(); j++) {
if(cur && sgn(v[j].first - pre[cur])) {
ans[cur] += areaarc(v[j].first - pre[cur], c[i].r);
ans[cur] += 0.5 * (c[i].point(pre[cur])^c[i].point(v[j].first));
}
cur += v[j].second;
pre[cur] = v[j].first;
}
}
}
}cs;
多边形
const int maxp=1013;
struct polygon {
int n;
Point p[maxp];
Line l[maxp];
void input(int n) {
this->n=n;
for(int i=0; i<n; i++) p[i].input();
}
// 得到所有边
void getline() {
for(int i=0; i<n; i++) {
l[i] = Line(p[i],p[(i+1)%n]);
}
}
// 以 p0 为标准极角排序
struct cmp {
Point p;
cmp(const Point &p0) {
p = p0;
}
bool operator()(const Point &aa, const Point &bb) {
Point a = aa, b = bb;
int d = sgn((a-p)^(b-p));
if(d == 0) {
return sgn(a.distance(p) - b.distance(p)) < 0;
}
return d > 0;
}
};
/*
进行极角排序
首先找打最左下角的点
需要重载好Point的 < 操作符
*/
void norm() {
Point mi = p[0];
for(int i=1; i<n; i++) mi = min(mi, p[i]);
sort(p, p+n, cmp(mi));
}
// 得到周长
db getcircumference() {
db sum = 0;
for(int i=0; i<n; i++) {
sum += p[i].distance(p[(i+1)%n]);
}
return sum;
}
// 得到面积
db getarea() {
db sum = 0;
// 以原点为划分点
for(int i=0; i<n; i++) {
sum += (p[i]^p[(i+1)%n]);
}
return fabs(sum)/2;
}
// 得到方向，1 表示逆时针，0表示顺时针
bool getdir() {
db sum = 0;
for(int i=0; i<n; i++) {
sum += (p[i] ^ p[(i+1)%n]);
}
if(sgn(sum) > 0) return 1;
return 0;
}
// 得到重心
Point getbarycentre() {
Point ret(0,0);
db area = 0;
for(int i=1; i<n-1; i++) {
db tmp = (p[i]-p[0])^(p[i+1]-p[0]);
if(sgn(tmp) == 0) continue;
area += tmp;
ret.x += (p[0].x + p[i].x + p[i+1].x) / 3 * tmp;
ret.y += (p[0].y + p[i].y + p[i+1].y) / 3 * tmp;
}
if(sgn(area)) ret = ret / area;
return ret;
}
/*
得到凸包 凸包点编号0 ~ n-1
*/
void getconvex(polygon &convex) {
sort(p, p+n);
convex.n = n;
for(int i=0; i<min(n, 2); i++) {
convex.p[i] = p[i];
}
if(convex.n == 2 && (convex.p[0] == convex.p[1])) convex.n --;
if(n <= 2) return;
int &top = convex.n;
top = 1;
for(int i=2; i<n; i++) {
while(top && sgn((convex.p[top] - p[i]) ^ (convex.p[top-1] - p[i])) <= 0)
top--;
convex.p[++top] = p[i];
}
int temp = top;
convex.p[++top] = p[n-2];
for(int i=n-3; i>=0; i--) {
while(top != temp && sgn((convex.p[top] - p[i]) ^ (convex.p[top-1] - p[i])) <= 0)
top--;
convex.p[++top] = p[i];
}
if(convex.n == 2 && (convex.p[0] == convex.p[1])) convex.n --;
convex.norm();
}
void Graham(polygon &convex) {
norm();
int &top = convex.n;
top = 0;
if(n == 1) {
top = 1;
convex.p[0] = p[0];
return ;
}
if(n == 2) {
top = 2;
convex.p[0] = p[0];
convex.p[1] = p[1];
if(convex.p[0] == convex.p[1]) top--;
return;
}
convex.p[0] = p[0];
convex.p[1] = p[1];
top = 2;
for(int i=2; i<n; i++) {
while(top > 1 && sgn((convex.p[top-1] - convex.p[top-2]) ^ (p[i]-convex.p[top-2])) <= 0)
top--;
convex.p[top++] = p[i];
}
if(convex.n == 2 && (convex.p[0] == convex.p[1])) convex.n--;
}
// 判断是不是凸的
bool isconvex() {
bool s[3];
memset(s, false, sizeof s);
for(int i=0; i<n; i++) {
int j = (i + 1) % n;
int k = (j + 1) % n;
s[sgn((p[j] - p[i]) ^ (p[k] - p[i])) + 1] = true;
if(s[0] && s[2]) return false;
}
return true;
}
/*
判断点和任意多边形的关系
3 点上
2 边上
1 内部
0 外部
*/
int relationpoint(Point q) {
for(int i=0; i<n; i++) {
if(p[i] == q) return 3;
}
getline();
for(int i=0; i<n; i++) {
if(l[i].pointonseg(q)) return 2;
}
int cnt = 0;
for(int i=0; i<n; i++) {
int j = (i + 1) % n;
int k = sgn((q-p[j])^(p[i]-p[j]));
int u = sgn(p[i].y - q.y);
int v = sgn(p[j].y - q.y);
if(k > 0 && u < 0 && v >= 0) cnt ++;
if(k < 0 && v < 0 && u >= 0) cnt --;
}
return cnt != 0;
}
//直线 u 切割凸多边形左侧
//注意直线方向
//HDU3982
void convexcut(Line u, polygon &po) {
int &top = po.n;
top = 0;
for(int i=0; i<n; i++) {
int d1 = sgn((u.e-u.s)^(p[i]-u.s));
int d2 = sgn((u.e-u.s)^(p[(i+1)%n]-u.s));
if(d1 >= 0) po.p[top++] = p[i];
if(d1 * d2 < 0) po.p[top++] = u.crosspoint(Line(p[i], p[(i+1)%n]));
}
}
// 多边形和圆交的面积
db areacircle(circle c) {
db ans = 0;
for(int i=0; i<n; i++) {
int j = (i + 1) % n;
if(sgn((p[j]-c.p)^(p[i]-c.p)) >= 0)
ans += c.areatriangle(p[i], p[j]);
else ans -= c.areatriangle(p[i], p[j]);
}
return fabs(ans);
}
/*
多边形与圆的关系
2. 圆完全在多边形内
1. 圆在多边形里面，碰到了多边形边界
0. 其他
*/
db relationcircle(circle c) {
getline();
int x = 2;
if(relationpoint(c.p) != 1) return 0; // 圆心不在内部
for(int i=0; i<n; i++) {
if(c.relationseg(l[i]) == 2) return 0;
if(c.relationseg(l[i]) == 1) x = 1; // 相切
}
return x;
}
db minRectangleCover() {
if(n < 3) return 0.0;
p[n] = p[0];
db ans = -1;
int up = 1, r = 1, l;
for(int i=0; i<n; i++) {
Point vec = p[i+1] - p[i];
while(sgn( (vec^(p[up+1]-p[i])) - (vec^(p[up]-p[i])) ) >= 0)
up = (up + 1) % n;
while(sgn( (vec*(p[r+1]-p[i])) - (vec * (p[r]-p[i])) ) >= 0)
r = (r + 1) % n;
if(i == 0) l = r;
while(sgn( (vec*(p[l+1]-p[i])) - (vec * (p[l]-p[i])) ) <= 0)
l = (l + 1) % n;
db d = (p[i] - p[i+1]).len2();
db tmp = (vec^(p[up]-p[i])) * ( (vec*(p[r]-p[i])) - (vec *(p[l]-p[i]))) / d;
if(ans < 0 || ans > tmp) ans = tmp;
}
return ans;
}
// 最远的一对点的距离
db diameter() {
if(n == 2) return p[0].distance(p[1]);
int i = 0, j = 0;
for(int k=0; k<n; k++) {
if(p[i] < p[k]) i = k;
if(!(p[j] < p[k])) j = k;
}
int si = i, sj = j;
db res = 0;
while(i != sj || j != si) {
res = max(res, p[i].distance(p[j]));
int ni = (i+1)%n;
int nj = (j+1)%n;
if(sgn((p[ni]-p[i])^(p[nj]-p[j])) <= 0) {
i = ni;
} else j = nj;
}
return res;
}
};
浙公网安备 33010602011771号