Fastest way to Shoelace formula

Question

I have made a function who calculate area polygon with Shoelace way.

That's works perfectly but right now I wonder if there is not a faster way to have the same result. I want to know that because this function must work faster with polygon with a lot of coordinates.

My function :

def shoelace_formula(polygonBoundary, absoluteValue = True):
    nbCoordinates = len(polygonBoundary)
    nbSegment = nbCoordinates - 1

    l = [(polygonBoundary[i+1][0] - polygonBoundary[i][0]) * (polygonBoundary[i+1][1] + polygonBoundary[i][1]) for i in xrange(nbSegment)]

    if absoluteValue:
        return abs(sum(l) / 2.)
    else:
        return sum(l) / 2.

My polygon :

polygonBoundary = ((5, 0), (6, 4), (4, 5), (1, 5), (1, 0))

Result :

22.

Any ideas?

I try with Numpy : It's speedest but you have to convert your coordinates first.

import numpy as np
x, y = zip(*polygonBoundary)

def shoelace_formula_3(x, y, absoluteValue = True):

    result = 0.5 * np.array(np.dot(x, np.roll(y, 1)) - np.dot(y, np.roll(x, 1)))
    if absoluteValue:
        return abs(result)
    else:
        return result

Bavardage

Using numpy should help. https://mcmap.net/q/102332/-calculate-area-of-polygon-given-x-y-coordinates – Bonaventure

@Jakub. I try with numpy, same performance. – Bavardage

Try it with many coordinates. When I tried both of your functions with 500 coordinate pairs, shoelace_formula_3 was twice as fast (115 microseconds) as shoelace_formula (321 microseconds). – Bonaventure

And if you do x, y = zip(*polygonBoundary) outside of the function and include x and y as function parameters, it runs in 93.7 microseconds. And import numpy outside of the function. – Bonaventure

Faster: 0.5 * np.abs(np.dot(x[:-1], y[1:]) + x[-1]*y[0] - np.dot(y[:-1], x[1:]) - y[-1]*x[0]) – Merth

MikeD in John Cook's blog suggested that the coordinates should be translated to be relative to one of the points in the polygon, in order to minimise precision loss when the absolute coordinate values are much larger than the area to be calculated. – Rickrickard

Answer 1

For me the fastest way would be using numpy where you have to send a numpy array of (x,y) cordinates as an argument in shoelace method:

import numpy as np
def shoelace(x_y):
    x_y = np.array(x_y)
    x_y = x_y.reshape(-1,2)

    x = x_y[:,0]
    y = x_y[:,1]

    S1 = np.sum(x*np.roll(y,-1))
    S2 = np.sum(y*np.roll(x,-1))

    area = .5*np.absolute(S1 - S2)

    return area

Aidaaidan

Welcome to SO! When you are replying to a question and there are some other answers, please provide the Cons and Pros on your answer vs remaining. – Staffordshire

Answer 2

Take a look at the Rosetta Code example using Numpy. Numpy gives a fast solution.

For your convenience, I paste below the Rosetta Code snippet:

import numpy as np
# x,y are arrays containing coordinates of the polygon vertices
x=np.array([3,5,12,9,5]) 
y=np.array([4,11,8,5,6]) 
i=np.arange(len(x))
#Area=np.sum(x[i-1]*y[i]-x[i]*y[i-1])*0.5 # signed area, positive if the vertex sequence is counterclockwise
Area=np.abs(np.sum(x[i-1]*y[i]-x[i]*y[i-1])*0.5) # one line of code for the shoelace formula

EDIT

You can now find the Shoelace formula implemented in the Scikit-Spatial library.

Alcahest

Your answer could be improved with additional supporting information. Please edit to add further details, such as citations or documentation, so that others can confirm that your answer is correct. You can find more information on how to write good answers in the help center. – Inhalant

I tested this answer, the one from @Mustaqim Khan and another implementation. I have not investigate why but this one gives a very different from the 2 other method I will suggest using Mustaqim Khan's answer – Randallrandan

Answer 3

Here's a version that uses 1/2 as many multiplications: https://mcmap.net/q/102333/-how-do-i-calculate-the-area-of-a-2d-polygon-duplicate

If you need even greater performance, you could consider doing this in a Python C extension. C can be dramatically faster than Python, especially for math operations -- sometimes 100-1000x.

Westernize

Thank you so much. It's better, but the numpy solution is speedest ! – Bavardage

Answer 4

Another interesting approach (although slower)

m = np.vstack([x,y])
result = 0.5 * np.abs(np.linalg.det(as_strided(m, (m.shape[1]-1, 2, 2), (m.itemsize, m.itemsize*m.shape[1], m.itemsize))).sum())

Merth

Answer 5

This is a very simple implementation of shoelace formula in python

class Polygon:
  def __init__(self,arr):
    self.arr = arr
  def area(self):
    total=0
    i = 0
    while i != len(self.arr)-1:
      total+=self.arr[i][0]*self.arr[i+1][1]
      total-=self.arr[i+1][0]*self.arr[i][1]
      i+=1
    return abs(total +self.arr[-1][0]*self.arr[0][1] -self.arr[-1][-1]*self.arr[0][0])/2

Enthusiast

Answer 6

Try simplest way, raw shoelace formula for triangles and polygons:

def shoelace_formula(x1, y1, x2, y2, x3, y3, x4, y4, x5, y5):
      return abs(0.5 * (x1*y2 + x2*y3 + x3*y4 + x4*y5 + x5*y1
                        - x2*y1 - x3*y2 - x4*y3 - x5*y4 - y1*y5))

print(shoelace_formula(5, 0, 6, 4, 4, 5, 1, 5, 1, 0))

Doordie

Answer 7

class Point //a new class for an any point a(X,Y), b(X,Y), c(X,Y), d(X,Y)
{
    //private int x, y;
    public int X { get; set; }
    public int Y { get; set; }

}

static class Geometry
{       

    public static void GerArea(Point a, Point b, Point c)
    {

        double area = 0.5 * ( (a.X * b.Y) + (b.X * c.Y) + (c.X * a.Y) - (b.X * a.Y) - (c.X * b.Y) - (a.X * c.Y) );

        Console.WriteLine(Math.Abs(area));
    }
    public static void GerArea(Point a, Point b, Point c, Point d)
    {
        double area = 0.5 * ((a.X * b.Y) + (b.X * c.Y) + (c.X * d.Y) + (d.X * a.Y) - (b.X * a.Y) - (c.X * b.Y) - (d.X * c.Y) - (a.X * d.Y ) );

        Console.WriteLine(Math.Abs(area));
    }
}
class Program
{
    static void Main(string[] args)
    {

        Point a = new Point() { X = -12, Y = 12 }; 
        Point b = new Point() { X = 15, Y = 15 };
        Point c = new Point() { X = -15, Y = -16 };
        Point d = new Point() { X = 16, Y = -15 };

        Console.WriteLine("****Shoelace formula****\n");


        Console.Write("Area of tringle: ");
        Geometry.GerArea(a, b, c);
        Console.Write("Area of quad: ");
        Geometry.GerArea(a, b, c, d);


        Console.ReadLine();

    }
}

Easley