half_plane.cpp

半平面求交。zzy大牛论文改写的。有一定的实现。但有些功能没有尝试过

#include <iostream>
#include <algorithm>
#include <cmath>
using namespace std;
double const mini = 1e-6;
double const pi = 3.141592653589793;
int const maxn = 20000; 
struct point{ double x, y; };
double arctan(double t){return atan2(t,1);}
double cross(point x, point y){return x.x * y.y - x.y * y.x;}
double cross(point a, point b, point c){
  point x, y;
  x.x = b.x - a.x; x.y = b.y - a.y;
  y.x = c.x - a.x; y.y = c.y - a.y;
  return cross(x, y);
}
struct Hplane{
  point a, b;
  double angle;
  void calcangle();
  bool included(point h){return cross(a, b, h) > mini;};
}; //half-planes
void Hplane::calcangle(){ //calculate the angle of a half - plane
  if (abs(b.x - a.x) < mini)
    if (b.y > a.y) angle = pi/2; else angle = pi * 3 / 2;
  else 
    if (b.x > a.x) angle = arctan((b.y - a.y) / (b.x - a.x));
    else angle = arctan((b.y - a.y) / (b.x - a.x)) + pi;
}
point calcPoint(point a, point b,point c, point d){
  double s1,s2;
  point p;
  s1=cross(a,b,c); s2=cross(a,b,d);
  p.x = (c.x*s2-d.x*s1) / (s2 - s1);
  p.y = (c.y*s2-d.y*s1) / (s2 - s1);
  return p;
}
int m, deque[maxn];
Hplane a[maxn];
point P0, ar[maxn];  //record the new points
int f = 0, r = 0;
bool lessthan(Hplane a, Hplane b){
  return a.angle < b.angle;
}
bool half_plane_ins(){
  for (int i = 0; i < m; i++)
      a[i].calcangle();
  sort(a, a+m, lessthan);
  deque[0]=0; r=1; f=0;
  for (int i = 1; i < m;){
    int j = i + 1;
    while (abs(a[i].angle - a[j].angle) < mini) {
          if (!a[j].included(a[i].a)) i = j;
          j++;
    }
    while (r > f + 1 && !a[i].included(ar[r-1])) r--;
    while (r > f + 1 && !a[i].included(ar[f+1])) f++;
    deque[r++] = i;
    ar[r-1] = calcPoint(a[deque[r-2]].a,a[deque[r-2]].b,
   a[deque[r-1]].a,a[deque[r-1]].b);
    i = j;
  }
  
  bool p;
  do{
    p = false;
    while (r > f + 1 && !a[deque[f]].included(ar[r-1])) r--, p =true;
    while (r > f + 1 && !a[deque[r-1]].included(ar[f+1])) f++, p =true;
  }while (p);
  //finish the polygon
  ar[r] = calcPoint(a[deque[r-1]].a,a[deque[r-1]].b,
      a[deque[f]].a,a[deque[f]].b);
  ar[f] = ar[r];
  if (f+1==r) return true; else 
  return false;
}