add HPI edge version

geometry/halfPlaneIntersection.cpp

@@ -1,74 +1,31 @@
-#define N 100010
-#define EPS 1e-8
-#define SIDE 10000000
-struct PO{ double x , y ; } p[ N ], o ;
-struct LI{
-  PO a, b;
-  double angle;
-  void in( double x1 , double y1 , double x2 , double y2 ){
-    a.x = x1 ; a.y = y1 ; b.x = x2 ; b.y = y2;
-  }
-}li[ N ] , deq[ N ];
-int n , m , cnt;
-inline int dc( double x ){
-  if ( x > EPS ) return 1;
-  else if ( x < -EPS ) return -1;
-  return 0;
-}
-inline PO operator-( PO a, PO b ){
-  PO c;
-  c.x = a.x - b.x ; c.y = a.y - b.y;
-  return c;
-}
-inline double cross( PO a , PO b , PO c ){
-  return ( b.x - a.x ) * ( c.y - a.y ) - ( b.y - a.y ) * ( c.x - a.x );
-}
-inline bool cmp( const LI &a , const LI &b ){
-  if( dc( a.angle - b.angle ) == 0 ) return dc( cross( a.a , a.b , b.a ) ) < 0;
-  return a.angle > b.angle;
-}
-inline PO getpoint( LI &a , LI &b ){
-  double k1 = cross( a.a , b.b , b.a );
-  double k2 = cross( a.b , b.a , b.b );
-  PO tmp = a.b - a.a , ans;
-  ans.x = a.a.x + tmp.x * k1 / ( k1 + k2 );
-  ans.y = a.a.y + tmp.y * k1 / ( k1 + k2 );
-  return ans;
-}
-inline void getcut(){
-  sort( li + 1 , li + 1 + n , cmp ); m = 1;
-  for( int i = 2 ; i <= n ; i ++ )
-    if( dc( li[ i ].angle - li[ m ].angle ) != 0 )
-      li[ ++ m ] = li[ i ];
-  deq[ 1 ] = li[ 1 ]; deq[ 2 ] = li[ 2 ];
-  int bot = 1 , top = 2;
-  for( int i = 3 ; i <= m ; i ++ ){
-    while( bot < top && dc( cross( li[ i ].a , li[ i ].b , getpoint( deq[ top ] , deq[ top - 1 ] ) ) ) < 0 ) top -- ;
-    while( bot < top && dc( cross( li[ i ].a , li[ i ].b , getpoint( deq[ bot ] , deq[ bot + 1 ] ) ) ) < 0 ) bot ++ ;
-    deq[ ++ top ] = li[ i ] ;
-  }
-  while( bot < top && dc( cross( deq[ bot ].a , deq[ bot ].b , getpoint( deq[ top ] , deq[ top - 1 ] ) ) ) < 0 ) top --;
-  while( bot < top && dc( cross( deq[ top ].a , deq[ top ].b , getpoint( deq[ bot ] , deq[ bot + 1 ] ) ) ) < 0 ) bot ++;
-  cnt = 0;
-  if( bot == top ) return;
-  for( int i = bot ; i < top ; i ++ ) p[ ++ cnt ] = getpoint( deq[ i ] , deq[ i + 1 ] );
-  if( top - 1 > bot ) p[ ++ cnt ] = getpoint( deq[ bot ] , deq[ top ] );
-}
-double px[ N ] , py[ N ];
-void read( int rm ) {
-  for( int i = 1 ; i <= n ; i ++ ) px[ i + n ] = px[ i ] , py[ i + n ] = py[ i ];
-  for( int i = 1 ; i <= n ; i ++ ){
-    // half-plane from li[ i ].a -> li[ i ].b
-    li[ i ].a.x = px[ i + rm + 1 ]; li[ i ].a.y = py[ i + rm + 1 ];
-    li[ i ].b.x = px[ i ]; li[ i ].b.y = py[ i ];
-    li[ i ].angle = atan2( li[ i ].b.y - li[ i ].a.y , li[ i ].b.x - li[ i ].a.x ) ;
-  }
-}
-inline double getarea( int rm ){
-  read( rm ); getcut();
-  double res = 0.0;
-  p[ cnt + 1 ] = p[ 1 ];
-  for( int i = 1 ; i <= cnt ; i ++ ) res += cross( o , p[ i ] , p[ i + 1 ] ) ;
-  if( res < 0.0 ) res *= -1.0;
-  return res;
-}
+int PtSide(Pt p, Line L) {
+    return dcmp((L.e - L.s)^(p - L.s));
+}
+bool argcmp(const Pt &a, const Pt &b) { // arg(a) < arg(b)
+    int f = (Pt{a.y, -a.x} > Pt{} ? 1 : -1) * (a != Pt{});
+    int g = (Pt{b.y, -b.x} > Pt{} ? 1 : -1) * (b != Pt{});
+    return f == g ? (a ^ b) > 0 : f < g;
+}
+vector<Line> HPI(vector<Line> P) {
+    sort(P.begin(), P.end(), [&](Line l, Line m) {
+        if (argcmp(l.v, m.v)) return true;
+        if (argcmp(m.v, l.v)) return false;
+        return PtSide(l.s, m) > 0;
+    });
+    int n = P.size(), l = 0, r = -1;
+    for (int i = 0; i < n; i++) {
+        if (i and !argcmp(P[i - 1].v, P[i].v)) continue;
+        while (l < r and PtSide(LLIntersect(P[r-1], P[r]), P[i]) <= 0) r--;
+        while (l < r and PtSide(LLIntersect(P[l], P[l+1]), P[i]) <= 0) l++;
+        P[++r] = P[i];
+    }
+    while (l < r and PtSide(LLIntersect(P[r-1], P[r]), P[l]) <= 0) r--;
+    while (l < r and PtSide(LLIntersect(P[l], P[l+1]), P[r]) <= 0) l++;
+    if (r - l <= 1 or !argcmp(P[l].v, P[r].v))
+        return {}; // empty
+    if (PtSide(LLIntersect(P[l], P[r]), P[l+1]) <= 0) {
+        assert(0);
+        return {}; // infinity
+    }
+    return vector(P.begin() + l, P.begin() + r + 1);
+}