Friday, May 31, 2013

Computing polygon area in C#

If you go to this page you'll see an algorithm for calculating the area of a polygon by considering the area of sections with respect to the X axis:

Each section consists of a rectangle and a triangle. The area of one section can be calculated based on two points A and B as follows:
        (area of rectangle)   (area of triangle on top)
 area = ((B.X - A.X) * A.Y) + (B.X-A.X) * (B.Y-A.Y) / 2
      = B.X*A.Y - A.X*A.Y + (B.X*B.Y - B.X*A.Y - A.X*B.Y + A.X*A.Y) / 2
      = B.X*A.Y/2 - A.X*A.Y/2 + (B.X*B.Y - A.X*B.Y) / 2
      = (B.X*A.Y - A.X*A.Y + B.X*B.Y - A.X*B.Y) / 2
Note that this formula works even if B is higher than A, and after you combine all the sections together, the total is translation invariant (still works if the polygon is not entirely in quadrant 1, the positive quadrant.)

We can simplify this formula further. First of all, instead of dividing by two when computing each section, we can skip that step and divide the total by 2 at the end. Secondly, let's consider two adjacent sections based on points A-B-C instead of just two points A-B:
 Area of A-B section * 2 = B.X*A.Y - A.X*A.Y + B.X*B.Y - A.X*B.Y
 Area of B-C section * 2 = C.X*B.Y - B.X*B.Y + C.X*C.Y - B.X*C.Y
 Area both sections  * 2 = B.X*A.Y - A.X*A.Y + B.X*B.Y - A.X*B.Y + C.X*B.Y - B.X*B.Y + C.X*C.Y - B.X*C.Y
 Simplified              = B.X*A.Y - A.X*A.Y - A.X*B.Y + C.X*B.Y + C.X*C.Y - B.X*C.Y
Notice that B.X*B.Y is positive in one section, and negative in the following section, so the two terms cancel out. In fact, if you extend this logic across the entire polygon, you'll find that all the P.X*P.Y terms cancel out for any point P. So really, for each pair of points A-B we only need to compute B.X*A.Y - A.X*B.Y! So that's what this algorithm does.
public static double PolygonArea(IEnumerable<PointD> polygon)
  var e = polygon.GetEnumerator();
  if (!e.MoveNext()) return 0;
  PointD first = e.Current, last = first;

  double area = 0;
  while (e.MoveNext()) {
    PointD next = e.Current;
    area += next.X * last.Y - last.X * next.Y;
    last = next;
  area += first.X * last.Y - last.X * first.Y;
  return area / 2;
The result is positive if the polygon is clockwise (assuming a coordinate system in which increasing Y goes upward), or negative if the polygon is counterclockwise.


Post a Comment

<< Home