计算几何基础小结
已解决问题:
判断点是否在线段上
判断两线段是否相交
判断点是否在多边形内
判断线段、折线、多边形是否在多变形内
判断上述是否在圆内
计算上述与线段及直线的交点
凸包
待解决:
半平面交
#include<ctime>
#include<cstdio>
#include<cstring>
#include<cstdlib>
#include<algorithm>
#include<cmath>
#include<vector>
#include<stack>
#include<queue>
using namespace std;
//================================ struct =================================
struct V{
double px,py;
V operator - (const V c){
V a;
a.px = px - c.px,a.py = py - c.py;
return a;
}
V operator + (const V c){
V a;
a.px = px + c.px,a.py = py + c.py;
return a;
}
V operator * (const double c){
V a = *this;
a.px *= c,a.py *= c;
return a;
}
V operator / (const double c){
V a = *this;
a.px /= c,a.py /= c;
return a;
}
bool operator <= (const V c){
//return px <= c.px && py <= c.py;
return px*px + py*py <= c.px*c.px + c.py*c.py;
}
bool operator == (const V c){
return px == c.px && py == c.py;
}
double MO2(){
return px * px + py * py;
}
double DIS(V p1,V p2,int pa){
if(pa == 1){
return abs(py - p1.py);
}
if(pa == 2){
return abs(px - p1.px);
}
if(pa == 0){
double k = (p2.py-p1.py)/(p2.px-p1.px),x0,y0;
x0 = (k*k * p1.px + k * (py - p1.py) + px) / (k*k + 1);
y0 = k * (px - p1.px) + p1.py;
return sqrt((px-x0)*(px-x0) + (py-y0)*(py-y0));
}
}
double CP(V c){
return px*c.py - py*c.px;
}//cross product
double DP(V c){
return px*c.px + py*c.py;
}//dot product
bool min(V a,V b){
if(b <= a) swap(a,b);
return a <= *this;
}
bool max(V a,V b){
if(a <= b) swap(a,b);
return *this <= a;
}
};//vector
struct L{
V p1,p2;
int p0;
/*
switch(p0){
case 1: beeline
case 2: segment
case 3: ray
}
*/
int pa;
/*
switch(pa){
case 1: parallel-y;
case 2: parallel-x;
case 0: k < oo && k != 0;
}
*/
double k;
void getk(){
if(p1.px == p2.px) pa = 1;
else if(p1.py == p2.py) pa = 2;
else pa = 0,k = (p1.py-p2.py) / (p1.px-p2.px);
}
bool IP(V q){
return q.min(p1,p2) && q.max(p1,p2);
}//in the ploygon of the line
bool PA(L c){
if(pa != c.pa) return false;
if(pa) return true;
else if(k == c.k) return true;
else return false;
}
bool CR(L c){
if((*this).PA(c)) return false;
V v1 = p2 - p1,v2 = c.p1 - p1,v3 = c.p2 - p1;
double exp = v1.CP(v2) * v1.CP(v3);
//if(v1.CP(v2) * v1.CP(v3) < 0) return true;
switch(p0){
case 1://B_S
return exp <= 0;
break;
case 2://S_S
return exp < 0;
break;
case 3://R_S
if(exp > 0) return false;
else if(!exp){
V v = v1.CP(v2)? c.p2:c.p1;
if(p1-v <= p2-v) return false;
else return true;
}
else{
V v4 = c.p1 - p2,v5 = c.p2 - p2;
double s1 = abs(v2.CP(v3)),s2 = abs(v4.CP(v5));
return s1 > s2;
}
break;
}
}//line cross
};//line
struct BL{
vector<L> E;
};//broken line
struct PL{
vector<L> E;
};//ploygon
struct O{
V o;
double r;
};//circle
//================================ data ===================================
double CP(V a,V b){
return a.px*b.py - a.py*b.px;
}//cross product
double DP(V a,V b){
return a.px*b.px + a.py*b.py;
}//dot product
double P_P_DIST(V a,V b){
return sqrt((a.px-b.px)*(a.px-b.px) + (a.py-b.py)*(a.py-b.py));
}//disdance for points
double P_L_DIST(V a,L b){
double ret = a.DIS(b.p1,b.p2,b.pa);
if(b.p0 != 1){
ret = min(ret,P_P_DIST(a,b.p1));
}
if(b.p0 == 2){
ret = min(ret,P_P_DIST(a,b.p2));
}
return ret;
}//least distance for point to line(segment/beeline/ray)
double P_BL_DIST(V p,BL bl){
double ret = oo;
for(int i = 0;i < bl.E.size();++ i){
ret = min(ret,P_L_DIST(p,bl.E[i]));
}
return ret;
}//least distance for point to broken line
double P_PL_DIST(V p,PL pl){
double ret = oo;
for(int i = 0;i < pl.E.size();++ i){
ret = min(ret,P_L_DIST(p,pl.E[i]));
}
return ret;
}//least distance for point to ploygon
double P_O_DIST(V p,O o){
return o.r - P_P_DIST(p,o.o);
}//least distance for point to circle
double AREA(PL pl){
V v0 = pl.E[0].p1;
double ret = 0;
for(int i = 0;i < pl.E.size();++ i){
ret += v0.CP(pl.E[i].p2-pl.E[i].p1);
v0 = v0 + (pl.E[i].p2-pl.E[i].p1);
}
return ret;
}//the measure of area
//================================= judge =================================
bool P_L_ON(V q,L p){
/*double x1 = p.p1.px,x2 = p.p2.px,y1 = p.p1.py,y2 = p.p2.py,
x0 = q.px,y0 = q.py;
if((!p.p0) && (!(x0 <= max(x1,x2) && x0 >= min (x1,x2) && y0 <= max(y1,y2) && y0 >= min(y1,y2))))
return false;*/
//if((!p.p0) && (!(q.min(p.p1,p.p2) && q.max(p.p1,p.p2)))) return false;
if((p.p0 == 2) && (!p.IP(q))) return false;
int pd = CP(p.p1-q,p.p2-q);
if(pd != 0) return false;
pd = DP(p.p1-q,p.p2-q);
if(pd <= 0) return true;
else{
if(p.p0 == 1) return true;
if(p.p0 == 2) return false;
if(p.p0 == 3){
if(p.p1-q <= p.p2-q) return false;
else return true;
}
}
//if(!CP(q-p.p1,p.p2-p.p1)) return true;
//else return false;
}//wheather a point is on a line(beeline/segment/ray)
bool S_S_CR(L a,L b){
return (a.CR(b) && b.CR(a)) || ((P_L_ON(a.p1,b) || P_L_ON(a.p2,b) || P_L_ON(b.p1,a) || P_L_ON(b.p2,a))/* && ((a.p2-a.p1).DP(b.p1-a.p1) * (a.p2-a.p1).DP(b.p2-a.p1) == 0)*/);
}//segments' cross
bool S_L_CR(L a,L b){
return b.CR(a);
}//segment and beeline's cross
bool S_R_CR(L a,L b){
return b.CR(a);
}//segment and ray's cross
bool P_PL_IN(V p,PL pl){
L o;
o.p1 = p,o.p2.px = p.px-1,o.p2.py = p.py,o.p0 = 3;
int cnt = 0;
for(int i = 0;i < pl.E.size();++ i){
if(P_L_ON(p,pl.E[i])) return true;
if(pl.E[i].pa != 1){
if(P_L_ON(pl.E[i].p1,o) || (P_L_ON(pl.E[i].p2,o))){
V v1 = P_L_ON(pl.E[i].p1,o)? pl.E[i].p1:pl.E[i].p2,
v2 = P_L_ON(pl.E[i].p1,o)? pl.E[i].p2:pl.E[i].p1;
if(v1.py > v2.py) cnt ++;
}
}
else if(S_R_CR(pl.E[i],o)) cnt ++;
}
return cnt & 1;
}//point in ploygon
bool S_PL_IN(L p,PL pl){
if(!(P_PL_IN(p.p1,pl) && P_PL_IN(p.p2,pl))) return false;
vector<V> v;
for(int i = 0;i < pl.E.size();++ i){
if(P_L_ON(p.p1,pl.E[i]) || P_L_ON(p.p2,pl.E[i])){
if(P_L_ON(p.p1,pl.E[i])) v.push_back(p.p1);
if(P_L_ON(p.p2,pl.E[i])) v.push_back(p.p2);
}
else if(P_L_ON(pl.E[i].p1,p) || P_L_ON(pl.E[i].p2,p)){
if(P_L_ON(pl.E[i].p1,p)) v.push_back(pl.E[i].p1);
if(P_L_ON(pl.E[i].p2,p)) v.push_back(pl.E[i].p2);
}
else if(p.CR(pl.E[i])) return false;
}
for(int i = 0;i < v.size() - 1;++ i){
if(!P_PL_IN((v[i] + v[i+1])/2,pl)) return false;
}
return true;
}//segment in ploygon
bool BL_PL_IN(BL bl,PL pl){
for(int i = 0;i < bl.E.size();++ i){
if(!S_PL_IN(bl.E[i],pl)) return false;
}
return true;
}//broken line in ploygon
bool PL_PL_IN(PL pl1,PL pl2){
for(int i = 0;i < pl1.E.size();++ i){
if(!S_PL_IN(pl1.E[i],pl2)) return false;
}
return true;
}//ploygon in ploygon
bool O_PL_IN(O o,PL pl){
double dis = oo;
for(int i = 0;i < pl.E.size();++ i){
dis = min(dis,P_L_DIST(o.o,pl.E[i]));
}
return o.r <= dis;
}//circle in ploygon
bool P_O_IN(V p,O o){
return (p-o.o).MO2() <= o.r * o.r;
}//point in circle
bool S_O_IN(L s,O o){
return P_O_IN(s.p1,o) && P_O_IN(s.p2,o);
}//segment in circle
bool BL_O_IN(BL bl,O o){
for(int i = 0;i < bl.E.size();++ i){
if(!P_O_IN(bl.E[i].p1,o)) return false;
}
return P_O_IN(bl.E[bl.E.size()-1].p2,o);
}//broken line in circle
bool PL_O_IN(PL pl,O o){
for(int i = 0;i < pl.E.size();++ i){
if(!P_O_IN(pl.E[i].p1,o)) return false;
}
return true;
}//ploygon in circle
bool O_O_IN(O o1,O o2){
if(o1.r > o2.r) return false;
else return o2.r-o1.r <= P_P_DIST(o1.o,o2.o);
}//circle in circle
//================================= node ==================================
void P_O_ND(V p,O o,double &x,double &y){
if(p == o.o){
x = y = oo;
return;
}
L a;
double lth = P_O_DIST(p,o);
a.p1 = o.o,a.p2 = p;
a.getk();
if(a.pa == 1){
if(p.py > o.o.py) x = p.px,y = o.o.py + o.r;
else x = p.px,y = o.o.py - o.r;
}
else if(a.pa == 2){
if(p.px > o.o.px) y = p.py,x = o.o.px + o.r;
else y = p.py,x = o.o.px - o.r;
}
else{
double xp = o.r / sqrt(a.k*a.k + 1),yp = abs(a.k * xp);
x = o.o.px,y = o.o.py;
if(p.px > o.o.px) x += xp;
else x -= xp;
if(p.py > o.o.py) y += yp;
else y -= yp;
}
}//closest node of point and circle
void sameL_S_S_ND(L a,L b,double &x,double &y){
V h = a.p1,v1,v2,v3,v4;
if(a.pa == 1) h.px += 100;
else if(a.pa == 2) h.py += 100;
else a.getk(),h.px += 100,h.py += 100 * a.k;
if(a.p1.MO2() > a.p2.MO2()) swap(a.p1,a.p2);
if(b.p1.MO2() > b.p2.MO2()) swap(b.p1,b.p2);
if(a.p2.MO2() > b.p2.MO2()) swap(a,b);
v1 = a.p1 - h,v2 = a.p2 - h,v3 = b.p1 - h,v4 = b.p2 - h;
double s1 = abs(CP(v1,v2)),s2 = abs(CP(v3,v4)),s3 = abs(CP(v1,v4));
if(s1 + s2 < s3) x = y = - oo;//no node
else if(s1 + s2 > s3) x = y = oo;//oo nodes
else x = a.p2.px,y = a.p2.py;
}//node of segments which are on the same beeline
void S_SorB_ND(L a,L b,double &x,double &y){
if(!((b.p0 == 1 && S_L_CR(a,b))||(b.p0 == 2 && S_S_CR(a,b)))){
x = y = - oo;//no node
return;
}
if(a.pa == 1){
if(b.pa == 1){
if(b.p0 == 1) x = y = oo;//oo nodes
else sameL_S_S_ND(a,b,x,y);
}
else{
x = a.p1.px;
b.getk();
y = b.k * (x - b.p1.px) + b.p1.py;
}
}
else if(b.pa == 1){
x = b.p1.px;
a.getk();
y = a.k * (x - a.p1.px) + a.p1.py;
}
else if(a.pa == 2){
if(b.pa == 2){
if(b.p0 == 1) x = y = oo;//oo nodes
else sameL_S_S_ND(a,b,x,y);
}
else{
y = a.p1.py;
b.getk();
x = 1.0/a.k * (y - b.p1.py) + b.p1.px;
}
}
else if(b.pa == 2){
y = b.p1.py;
a.getk();
x = 1.0/b.k * (y - a.p1.py) + a.p1.px;
}
else{
a.getk(),b.getk();
if(a.k == b.k) sameL_S_S_ND(a,b,x,y);
else{
x = (a.k*a.p1.px - b.k*b.p1.px + b.p1.py - a.p1.px) / (a.k - b.k);
y = a.k * (x - a.p1.px) + a.p1.py;
}
}
}//node of segment and segment(or beeline)
void SorB_BL_ND(L l,BL bl,V v[]){
int j = 1;
for(int i = 0;i < bl.E.size();++ i){
S_SorB_ND(l,bl.E[i],v[j].px,v[j].py);
if(v[j].px != oo && v[j].px != -oo)
j ++;
}
}//node(s) of segment(or beeline) and broken line
void SorB_PL_ND(L l,PL pl,V v[]){
int j = 1;
for(int i = 0;i < pl.E.size();++ i){
S_SorB_ND(l,pl.E[i],v[j].px,v[j].py);
if(v[j].px != oo && v[j].px != -oo)
j ++;
}
}//node(s) of segment(or beeline) and ploygon
void SorB_O_ND(L l,O o,V v[]){
if(l.pa == 2 && P_O_IN(l.p1,o) && P_O_IN(l.p2,o)) return;
//L l0 = l;l0.pa = 1;
double dis = P_L_DIST(o.o,l);
V v1,v2;
if(dis > o.r) return;
if(l.p0 == 1){
double h = sqrt(o.r * o.r - dis * dis);
if(h == 0){
v[1].px = l.p1.px,v[1].py = o.o.py;
}
else{
v1.px = v2.px = l.p1.px;
v1.py = v2.py = o.o.py;
v1.py += h,v2.py -= h;
if(l.pa == 1){
v[1] = v1,v[2] = v2;
}
else{
if((P_O_IN(l.p1,o) && P_O_IN(l.p2,o))||(!P_O_IN(l.p1,o) && !P_O_IN(l.p2,o))) return;
if(P_L_ON(v1,l) && P_L_ON(v2,l)){
v[1] = v1,v[2] = v2;
}
else if(P_L_ON(v1,l)){
v[1] = v1;
}
else if(P_L_ON(v2,l)){
v[1] = v2;
}
}
}
}
else if(l.p0 == 2){
double w = sqrt(o.r * o.r - dis * dis);
if(w == 0){
v[1].py = l.p1.py,v[1].px = o.o.px;
}
else{
v1.py = v2.py = l.p1.py;
v1.px = v2.px = o.o.px;
v1.px += w,v2.px -= w;
if(l.pa == 1){
v[1] = v1,v[2] = v2;
}
else{
if((P_O_IN(l.p1,o) && P_O_IN(l.p2,o))||(!P_O_IN(l.p1,o) && !P_O_IN(l.p2,o))) return;
if(P_L_ON(v1,l) && P_L_ON(v2,l)){
v[1] = v1,v[2] = v2;
}
else if(P_L_ON(v1,l)){
v[1] = v1;
}
else if(P_L_ON(v2,l)){
v[1] = v2;
}
}
}
}
else{
l.getk();
double A = l.k * l.k + 1,B = 2 * (-l.k * l.p1.px + l.p1.py - o.o.py - o.o.px),C = (-l.k * l.p1.px + l.p1.py - o.o.py) * (-l.k * l.p1.px + l.p1.py - o.o.py) - o.r * o.r,D;// = 4 * ((-l.k * l.p1.px + l.p1.py - o.o.py - o.o.px) * (-l.k * l.p1.px + l.p1.py - o.o.py - o.o.px) - (l.k * l.k + 1) * ((-l.k * l.p1.px + l.p1.py - o.o.py) * (-l.k * l.p1.px + l.p1.py - o.o.py) - o.r * o.r));
D = B * B - 4 * A * C;
v1.px = (-B + sqrt(D))/(2*A),v1.py = l.k * (v1.px - l.p1.px) + l.p1.py;
v2.px = (-B - sqrt(D))/(2*A),v2.py = l.k * (v2.px - l.p1.px) + l.p1.py;
if(D == 0){
v[1] = v1;
}
else{
if(l.pa == 1){
v[1] = v1,v[2] = v2;
}
else{
if((P_O_IN(l.p1,o) && P_O_IN(l.p2,o))||(!P_O_IN(l.p1,o) && !P_O_IN(l.p2,o))) return;
if(P_L_ON(v1,l) && P_L_ON(v2,l)){
v[1] = v1,v[2] = v2;
}
else if(P_L_ON(v1,l)){
v[1] = v1;
}
else if(P_L_ON(v2,l)){
v[1] = v2;
}
}
}
}
}//node(s) of segment(or beeline) and circle
//============================ convex closure =============================
bool cmp(V a,V b){
return CP(a,b) > 0;
}
stack<V,vector<V> > GCC(V src[],int n){
V p0 = src[1];
int o0 = 1;
for(int i = 2;i <= n;++ i){
if(p0.py != src[i].py){
p0 = p0.py < src[i].py? p0:(o0 = i,src[i]);
}
else p0 = p0.px < src[i].py? p0:(o0 = i,src[i]);
}
V tmp[maxn];
int j = 0;
for(int i = 1;i <= n;++ i){
if(i != o0){
tmp[++j] = src[i];
}
}
sort(tmp+1,tmp+n,cmp);
stack<V,vector<V> > S;
S.push(p0),S.push(tmp[1]),S.push(tmp[2]);
V ut = tmp[1];//ut: undertop
for(int i = 3;i < n;++ i){
while(CP(S.top() - ut,tmp[i] - S.top()) <= 0){
S.pop(),S.pop();
V tp = ut;ut = S.top();
S.push(tp);
}
S.push(tmp[i]);
}
return S;
}//get convex closure
//================================= main ==================================
const double OO = 1e-8;
const double oo = 1e9;
const int maxn = 10000;
void init(){
}
void work(){
}
int main(){
init();
work();
return 0;
}