dmkz's blog

I implemented a triangulation of polygon in C++ in time O(n log(n))
By dmkz, history, 14 months ago, translation, In English

Hello, codeforces!

Today, based on the article Triangulation of polygon from a book "Computational Geometry: Algorithms and applications", I implemented an algorithm of triangulation of any polygon without self intersections (non-convex polygons too) in C++.

$$$\text{ }$$$ $$$\text{ }$$$

Why it was needed for me? Probably, you know about a competition RuCode 7.0. In this year, qualification round contained next problem H: you are given an endless plane, painted in black and white stripes with height 1 and infinity width. These stripes alternate: black, white, black, white... You are given a polygon on this plane as list of $$$n\text{ }\left(3\le n \le 100000\right)$$$ 2D points with integer coordinates $$$x_i, y_i \left(-10^9 \le x_i, y_i \le 10^9\right)$$$ in some traversal order. Need to calculate a total square inside this polygon which is painted in black color.

Firstly, I solved for any triangle ($$$n=3$$$). Secondly, I thought about the triangulation of polygon and accumulating answers for each triangle.

So, the heaviest part of solution is the triangulation of polygon. I implemented this, stress-tested and fixed all of incomprehensible and dubious places of pseudo code that can be found in above book. For example, the author of this book states that vertex which has a type "merge" should be connected with another "merge" vertex, and this is used in pseudo-code, but there is an example, where it is not true. Stress-testing shows all of corner cases, which are not covered by the pseudo code.

I fixed all of mistakes which the stress-testing has found. You can just call this function:

std::vector<Triangle> triangulate(vpii points) {
    makeCounterClockwise(points);
    return triangulateFast(points);
}

The output of this function is a vector of triangles. Be careful, because I'm erasing non-meaningful vertices like if $$$P_{i-1}$$$, $$$P_{i}$$$ and $$$P_{i+1}$$$ lie on the same line than erase $$$P_i$$$.

Source code of solution for given problem

This solution works in $$$296$$$ ms in Yandex.Contest when $$$n=10^5$$$.

UPD. updated the source code by removing all unused things and commented lines

UPD 2. orz supposed $$$O(N)$$$ solution of original problem H in russian comments for this blog. Lets use same trick which is used in an algorithm of calculation of square of polygon. We can fix a point $$$P=(0,0)$$$ and calculate the answer for each triangle $$$P A_i A_{i+1}$$$ (for each side of polygon): if the order of traversak of current triangle was clockwise, then we can take this answer with a sign "minus", otherwise we can take this answer with a sign "plus".

Linear solution of problem H

UPD 3. k1r1t0 supposed to use trapezoids instead of triangles (in russian comments for this blog). Seems that for trapezoids a square which is painted in black can be calculated simplier.

  • Vote: I like it
  • +114
  • Vote: I do not like it

»
14 months ago, # |
  Vote: I like it +3 Vote: I do not like it

Wow I got chills, literally about 10 hours ago I was looking for a valid Subquadratic Polygon Triangulation after seeing the lecture 2 days ago.

But I found this code by Takanori Maehara which does the same thing in $$$O(N)$$$

real triangulate(vector<point> ps) {
  int n = ps.size();
  vector<int> next(n);
  for (int i = 0; i < n-1; ++i) next[i] = i+1;
  auto is_ear = [&](int i, int j, int k) {
    if (sign(cross(ps[j]-ps[i], ps[k]-ps[i])) <= 0) return false;
    for (int l = next[k]; l != i; l = next[l])
      if (sign(cross(ps[i]-ps[l], ps[j]-ps[l])) >= 0
       && sign(cross(ps[j]-ps[l], ps[k]-ps[l])) >= 0
       && sign(cross(ps[k]-ps[l], ps[i]-ps[l])) >= 0) return false;
    return true;
  };
  real area = 0;
  for (int i = 0; next[next[i]] != i; ) {
    if (is_ear(i, next[i], next[next[i]])) {
      area  += abs(cross(ps[next[i]]-ps[i], ps[next[next[i]]] - ps[i])) / 2;
      next[i] = next[next[i]];
    } else i = next[i];
  }
  return area;
}
→ Reply
  • »
    »
    14 months ago, # ^ |
      Vote: I like it 0 Vote: I do not like it

    How is this $$$O(N)$$$?

    → Reply
  • »
    »
    14 months ago, # ^ |
      Vote: I like it +11 Vote: I do not like it

    Without even looking at the code, $$$O(n)$$$ triangulation is known to be notoriously complicated algorithm. Possible, but definitely not as simple

    → Reply
    • »
      »
      »
      14 months ago, # ^ |
        Vote: I like it 0 Vote: I do not like it

      my bad. I think it's $$$N^2$$$

      → Reply
  • »
    »
    14 months ago, # ^ |
      Vote: I like it -8 Vote: I do not like it

    https://github.com/spaghetti-source/algorithm/blob/master/geometry/_geom.cc#L729

    → Reply
  • »
    »
    14 months ago, # ^ |
    ← Rev. 2   Vote: I like it 0 Vote: I do not like it

    This approach can be used in $$$O(N)$$$ solution of original problem H without triangulation. I updated the blog with a source code: if order of traversal of triangle $$$O A_i A_{i+1}$$$ is clockwise, then we can take the answer for this triangle with a sign "minus", otherwise — with a sign "plus".

    Also, seems that there is an easier solution with trapezoids instead of triangles.

    → Reply