Area of a polygon (Recursively using Python)


The essence of recursion is as follows:

  1. Identify a simple base case and solve for that.
  2. Conceive of a process that, when repeated, reduces a more complex case to that base case.
  3. Apply the process in #2 to the problem until you reach the base case.

In your case, the first step is easy. The smallest polygon is a triangle. The area of a triangle is (x1y2 + x2y3 + x3y1 – y1x2 –y2x3 – y3x1) / 2. (It looks like they misstated it in the problem though...)

The second step is also easy, because the problem statement gives it to you: given an n-vertex polygon, lop off a triangle, determine its area, and add it to the area of the resulting (n-1)-vertex polygon.

We'll break it down into parts. First, a function to solve #1:

def area_of_triangle(points):
    (x1, y1), (x2, y2), (x3, y3) = points
    return abs(x1 * y2 + x2 * y3 + x3 * y1 - y1 * x2 - y2 * x3 - y3 * x1) / 2

Easy. Now a function to solve #2. All we need is a function that lops off a triangle and returns both it and the resulting smaller polygon:

def lop_triangle(points):
    triangle = [points[0], points[-1], points[-2]]
    polygon = points[:-1]
    return triangle, polygon

If it's not obvious, this simply creates a triangle using the first and the last two points of the polygon. Then it removes the last point of the polygon, which is now equivalent to chopping off the triangle. (Draw a n-polygon and label its vertices from 0 to n to see how it works.) So now we have a triangle and a simpler polygon.

Now, let's put it all together. This third step is in some ways the hardest, but because we solved the first two problems already, the third is easier to understand.

def area_of_polygon(points):
    if len(points) == 3:
        return area_of_triangle(points)
        triangle, polygon = lop_triangle(points)
        return area_of_triangle(triangle) + area_of_polygon(polygon)

All the magic happens in that last line. Every time area_of_polygon gets a triangle, it just returns the area of a triangle. But when it gets a larger polygon, it lops off a triangle, takes the area of that triangle, and adds it to... the area of a smaller polygon. So say the polygon has 5 vertices. The first time area_of_polygon is called (c1), it lops off a triangle, takes its area, and then calls area_of_polygon (c2) again, but this time with a 4-vertex polygon. Then area_of_polygon lops of a triangle, and calls area_of_polygon (c3) again, but this time with a 3-vertex polygon. And then it doesn't have to call area_of_polygon again. It just returns the area of a triangle to the previous call (c2). That sums the result with the triangle in (c2) and returns that value to (c1). And then you have your answer.

Also, for what it's worth, the wolfram formula can be written with great clarity in three lines:

def area_of_polygon(vertices):
    pairs = zip(vertices, vertices[1:] + vertices[0:1])
    return sum(x1 * y2 - y1 * x2 for (x1, y1), (x2, y2) in pairs) / 2

The implementation of your formula is flawed. It looks ahead to values in your x, y lists that haven't been set yet with (x[i]*y[i+1] - x[i+1]*y[i])

If you put a print statement inside your try-except block you will see that you are simply multiplying by zero and getting zero area:

    x[i], y[i] = polygon[i]
    area += (x[i]*y[i+1] - x[i+1]*y[i])
    print x[i], y[i+1], x[i+1], y[i]
except IndexError, e:
    return area/2

#1 0 0 2
#2 0 0 5

Also, you are not returning the results of your recursive call to areaofpolygon, so you will never get that area/2. You want: return areaofpolygon(polygon, i+1). And make sure you actually divide by 2.0 so that you get float precision since your points are ints.

Try just using the formula you found or that was suggested in another question.


Here is a fixed version of your code:

#!/usr/bin/env python

from random import randint
from shapely.geometry import Polygon

area = 0

def areaofpolygon(polygon, i):
    global area
    if i == 0: 
        area = 0

        x1, y1 = polygon[i]
        x2, y2 = polygon[i+1]
        area += (x1*y2) - (x2*y1)

    except IndexError, e:
        x1, y1 = polygon[0]
        x2, y2 = polygon[-1]
        area += (x2*y1) - (x1*y2)
        return abs(area/2.0)

    return areaofpolygon(polygon, i+1)   

def main():
    mypolygon = [(randint(0, 100), randint(0, 100)) for _ in xrange(10)]
    print mypolygon

    area = areaofpolygon(mypolygon, 0)
    assert area == Polygon(mypolygon).area

    print "Test passed."
    return 0

if __name__ == '__main__':

### Output ###
$ ./ 
[(29, 76), (85, 49), (27, 80), (94, 98), (19, 1), (75, 6), (55, 38), (74, 62), (0, 25), (93, 94)]
Test passed.
$ ./ 
[(13, 37), (98, 74), (42, 58), (32, 64), (95, 97), (34, 62), (34, 59), (21, 76), (55, 32), (76, 31)]
Test passed.
$ ./ 
[(38, 67), (66, 59), (16, 71), (53, 100), (64, 52), (69, 31), (45, 23), (52, 37), (27, 21), (42, 74)]
Test passed.

Notice that you didn't need global x,y lists. And you also missed the last part of the equation where you use the last point and the first point.

Use this formula.

You accomplish your task in one for loop.