Bentley-Ottman Segment Intersection O((n+k) log n)

bentley-ottman.cc

#include <iostream>
#include <vector>
#include <cstdio>
#include <iomanip>
#include <algorithm>
#include <cmath>
#include <map>
#include <cassert>
#include <queue>
#include <set>
#include <functional>
#include <ctime>
#include <numeric>

using namespace std;

#define fst first
#define snd second
#define all(c) ((c).begin()), ((c).end())
#define TEST(x,a) { auto y=(x); if (sign(y-a)) { cout << "line " << __LINE__  << #x << " = " << y << " != " << a << endl; exit(-1); } }

#define dout if(0)cout

double urand() { return rand() / (1.0 + RAND_MAX); }

const double PI = acos(-1.0);

// implementation note: use EPS only this function
// usage note: check sign(x) < 0, sign(x) > 0, or sign(x) == 0
// notice: should be normalize to O(1)
const double EPS = 1e-8;
int sign(double x) {
  if (x < -EPS) return -1;
  if (x > +EPS) return +1;
  return 0;
}
struct point {
  typedef double T;
  T x, y; 
  point &operator+=(point p) { x += p.x; y += p.y; return *this; }
  point &operator*=(T a)     { x *= a;   y *= a;   return *this; }
  point &operator-=(point p) { x -= p.x; y -= p.y; return *this; }
  point &operator/=(T a)     { return *this *= (1.0/a); }
  point operator-() const    { return {-x, -y}; }
  bool operator<(point p) const { 
    if (sign(x - p.x)) return sign(x - p.x) < 0;
    return sign(y - p.y) < 0;
  }
};
bool operator>(point p, point q) { return q<p; }
bool operator<=(point p, point q) { return !(q<p); }
bool operator>=(point p, point q) { return !(p<q); }
bool operator==(point p, point q) { return !(p<q) && !(q<p); }
bool operator!=(point p, point q) { return (p<q) || !(q<p); }
point operator+(point p, point q) { return p += q; }
point operator-(point p, point q) { return p -= q; }
point operator*(point::T a, point p) { return p *= a; }
point operator*(point p, point::T a) { return p *= a; }
point operator/(point p, point::T a) { return p /= a; }
point::T dot(point p, point q) { return p.x*q.x+p.y*q.y; }
point::T cross(point p, point q) { return p.x*q.y-p.y*q.x; } // left turn > 0
point normalize(point p) { 
  auto s = sqrt(dot(p, p));
  return sign(s) ? p / s : point({0,0});
}
point::T norm2(point p) { return dot(p,p); }
point orth(point p) { return {-p.y, p.x}; }
point::T norm(point p) { return sqrt(dot(p,p)); }

ostream &operator<<(ostream &os, const point &p) { os<<"("<<p.x<<","<<p.y<<")"; return os; }

struct segment { point p, q; };

vector<point> intersect(segment s, segment t) {
  auto a = cross(s.q - s.p, t.q - t.p);
  auto b = cross(t.p - s.p, t.q - t.p); 
  auto c = cross(s.q - s.p, s.p - t.p);
  if (a < 0) { a = -a; b = -b; c = -c; }
  if (sign(b) < 0 || sign(a-b) < 0 ||
      sign(c) < 0 || sign(a-c) < 0) return {};      // disjoint
  if (sign(a) != 0) return {s.p + b/a*(s.q - s.p)}; // properly crossing
  vector<point> ps;                                 // same line
  auto insert_if_possible = [&](point p) {
    for (auto q: ps) if (sign(dot(p-q, p-q)) == 0) return;
    ps.push_back(p);
  };
  if (sign(dot(s.p-t.p, s.q-t.p)) <= 0) insert_if_possible(t.p);
  if (sign(dot(s.p-t.q, s.q-t.q)) <= 0) insert_if_possible(t.q);
  if (sign(dot(t.p-s.p, t.q-s.p)) <= 0) insert_if_possible(s.p);
  if (sign(dot(t.p-s.q, t.q-s.q)) <= 0) insert_if_possible(s.q);
  return ps;
}

//
// Verified: AOJ1226, AOJ2448
//
struct arrangement {
  struct Vertex; struct Edge; // Doubly Connected Edge List
  struct Vertex {
    point p;    // 
    Edge *edge; // any edge incident to this vertex
  };
  struct Edge { 
    Vertex *vertex;    // origin vertex of the edge
    Edge *twin;        // reverse edge of e
    Edge *prev, *next; // surrounding edges of face (CCW)
  };
  map<Vertex*, map<Vertex*, Edge*>> adj;
  vector<Vertex*> vertices;
  vector<Edge*> edges;

  vector<segment> segs;

  struct node { // Sweep-Line Structure (RBST)
    int index;
    int size = 1;
    node *left = 0, *right = 0;
  } *root = 0;
  vector<node> ns;

  node *update(node *x) {
    if (x) {
      x->size = 1;
      if (x->left)  x->size += x->left->size;
      if (x->right) x->size += x->right->size;
    }
    return x;
  }
  node *merge(node *x, node *y) {
    if (!x) return y;
    if (!y) return x;
    if (rand() % (x->size + y->size) < x->size) {
      x->right = merge(x->right, y);
      return update(x);
    } else {
      y->left = merge(x, y->left);
      return update(y);
    }
  }
  template <class C> // 3-way split: cond(x) < 0, cond(x) == 0, cond(x) > 0
  tuple<node*, node*, node*> split(node *x, C cond) {
    if (!x) return make_tuple(x,x,x);
    if (cond(x) == 0) {
      auto a = split(x->left, cond);
      auto b = split(x->right, cond);
      x->left = x->right = 0; update(x);
      get<1>(a) = merge(merge(get<1>(a), x), get<1>(b));
      get<2>(a) = get<2>(b);
      return a;
    }
    if (cond(x) < 0) {
      auto a = split(x->right, cond);
      x->right = 0; update(x);
      get<0>(a) = merge(x, get<0>(a));
      return a;
    }
    if (cond(x) > 0) {
      auto a = split(x->left, cond);
      x->left = 0; update(x);
      get<2>(a) = merge(get<2>(a), x);
      return a;
    }
  }
  node *leftmost(node *x) { while (x && x->left) x = x->left; return x; }
  node *rightmost(node *x) { while (x && x->right) x = x->right; return x; }
  template <class F>
  void process(node *x, F func) {
    if (!x) return;
    process(x->left, func);
    func(x);
    process(x->right, func);
  }

  arrangement(vector<segment> segs_) : segs(segs_) {
    ns.resize(segs.size());

    set<point> events;
    map<point, set<int>> L, R;
    for (int i = 0; i < segs.size(); ++i) {
      if (segs[i].q < segs[i].p) swap(segs[i].p, segs[i].q);
      events.insert(segs[i].p);
      events.insert(segs[i].q);
      L[segs[i].p].insert(i);
      R[segs[i].q].insert(i);
      ns[i].index = i;
    }
    vector<Vertex*> last(segs.size());

    while (!events.empty()) {
      const point p = *events.begin();
      events.erase(events.begin());

      Vertex *u = new Vertex({p});
      vertices.push_back(u);
      
      auto cond = [&](node *x) {
        const segment &s = segs[x->index];
        if (sign(s.q.x - s.p.x) == 0) {
          if (sign(p.y - s.p.y) < 0) return -1;
          if (sign(s.q.y - p.y) < 0) return +1;
          return 0;
        }
        return -sign(cross(normalize(s.p - p), normalize(s.q - p)));
      };
      auto z = split(root, cond);
      vector<node*> inserter;
      process(get<1>(z), [&](node *x) {
        Vertex *v = last[x->index];
        if (!adj[u][v]) {
          adj[u][v] = new Edge({u});
          adj[v][u] = new Edge({v});
          adj[u][v]->twin = adj[v][u];
          adj[v][u]->twin = adj[u][v];
          edges.push_back(adj[u][v]);
          edges.push_back(adj[v][u]);
        }
        if (!R[p].count(x->index)) 
          inserter.push_back(x);
      });
      for (int i: L[p]) 
        if (!R[p].count(i))
          inserter.push_back(&ns[i]);

      sort(all(inserter), [&](node *x, node *y) {
        const segment &s = segs[x->index], &t = segs[y->index];
        return sign(cross(normalize(s.q - s.p), normalize(t.q - t.p))) >= 0;
      });
      auto add_event = [&](node *x, node *y) {
        if (!x || !y) return;
        vector<point> ps = intersect(segs[x->index], segs[y->index]);
        for (point q: ps) 
          if (p < q) events.insert(q);
      };
      if (inserter.empty()) {
        add_event(rightmost(get<0>(z)), leftmost(get<2>(z)));
      } else {
        add_event(rightmost(get<0>(z)), inserter[0]);
        add_event(leftmost(get<2>(z)), inserter.back());
      }
      root = 0;
      for (int i = 0; i < inserter.size(); ++i) {
        last[inserter[i]->index] = u;
        inserter[i]->left = inserter[i]->right = 0;
        root = merge(root, update(inserter[i]));
      }
      root = merge(merge(get<0>(z), root), get<2>(z));
    }
    for (auto &pp: adj) {
      Vertex *u = pp.fst;
      vector<Edge*> es;
      for (auto z: pp.snd) es.push_back(z.snd);
      sort(all(es), [&](Edge *e, Edge *f) {
        auto quad = [](point p) {
          for (int i = 1; i <= 4; ++i, swap(p.x = -p.x, p.y))
            if (p.x > 0 && p.y >= 0) return i;
          return 0;
        };
        const point p = e->twin->vertex->p - e->vertex->p;
        const point q = f->twin->vertex->p - f->vertex->p;
        if (quad(p) != quad(q)) return quad(p) < quad(q);
        return cross(normalize(p), normalize(q)) > 0;
      });
      u->edge = es.back();
      for (Edge *e: es) {
        u->edge->next = e;
        u->edge->next->prev = u->edge;
        u->edge = u->edge->next;
      }
    }
  }
};

void AOJ1226() {
  for (int n; ~scanf("%d", &n) && n; ) {
    //cout << n << endl;
    vector<vector<double>> a(4, vector<double>(n));
    for (int k = 0; k < 4; ++k)
      for (int i = 0; i < n; ++i)
        scanf("%lf", &a[k][i]);

    vector<segment> ss = {
      {{0,0},{0,1}},
      {{0,1},{1,1}},
      {{1,1},{1,0}},
      {{1,0},{0,0}}
    };
    for (int i = 0; i < n; ++i) {
      ss.push_back({{a[0][i],0},{a[1][i],1}});
      ss.push_back({{0,a[2][i]},{1,a[3][i]}});
    }
    arrangement arr(ss);

    double result = 0;
    set<arrangement::Edge*> seen;
    for (auto *e: arr.edges) {
      if (seen.count(e)) continue;
      auto *f = e;
      double area = 0;
      do {
        seen.insert(f);
        area += cross(f->vertex->p, f->twin->vertex->p);
        f = f->twin->prev;
      } while (f != e);
      result = max(result, area);
    }
    printf("%.6f\n", result/2);
  }
}
// 
// s.p + a (s.q - s.p) == t.p + b (t.q - t.p)
// を解いて左辺を返す．未知変数は２つ a の近似値 a* をもつ
//

void AOJ2448() {
  int n; scanf("%d", &n);
  vector<point> ps(n);
  for (int i = 0; i < n; ++i)
    scanf("%lf %lf", &ps[i].x, &ps[i].y);
  vector<segment> ss;
  for (int i = 0; i+1 < n; ++i)
    ss.push_back({ps[i], ps[i+1]});
  arrangement arr(ss);

  double result = 0;
  set<arrangement::Edge*> seen;
  for (auto *e: arr.edges) {
    if (seen.count(e)) continue;
    auto *f = e;
    double area = 0;
    do {
      seen.insert(f);
      area += cross(f->vertex->p, f->twin->vertex->p);
      f = f->twin->prev;
    } while (f != e);
    if (area > 0) result += area;
  }
  printf("%.12f\n", result/2);
}

int main() {
  //AOJ1226();
  AOJ2448();
  return 0;
  /*
  for (int seed =0; seed < 10000; ++seed) {
    cout << "seed = " << seed << endl;
    srand(seed);
    int n = 100;
    vector<segment> ss;
    int x = 0, y = 0;
    for (int i = 0; i < n; ++i) {
      int z = rand() % n, w = rand() % n;
      ss.push_back({{x,y},{z,w}});
      x = z; y = w;
      //cout << ss.back().p << " " << ss.back().q << endl;
    }
    arrangement arr(ss);
  }
  */
  /*
  ss.push_back({{0,0},{3,2}});
  ss.push_back({{0,1},{3,1}});
  ss.push_back({{2,0},{2,4}});
  //ss.push_back({{0,3},{3,2}});
  //ss.push_back({{1.5,0},{1.5,4}});
  //
  */

  /*
  auto *v = arr.get({1,1});
  auto *e = v->edge;
  do { 
    cout << e->vertex->p << " -- " << e->twin->vertex->p << endl;
    e = e->next;
  } while (e != v->edge);
  */
}