Shoelace formula for polygonal area: Difference between revisions

Content added Content deleted

Inline

Revision as of 12:32, 4 October 2017

Task

Write a function/method/routine to use the the Shoelace formula to calculate the area of the polygon described by the ordered points:

     (3,4), (5,11), (12,8), (9,5), and (5,6)

Show the answer here, on this page.

ALGOL 68

<lang algol68>BEGIN

   # returns the area of the polygon defined by the points p using the Shoelace formula #
   OP  AREA = ( [,]REAL p )REAL:
       BEGIN
           [,]REAL points = p[ AT 1, AT 1 ]; # normalise array bounds to start at 1 #
           IF 2 UPB points /= 2 THEN
               # the points do not have 2 coordinates #
               -1
           ELSE
               REAL   result := 0;
               INT    n       = 1 UPB points;
               IF n > 1 THEN
                   # there at least two points #
                   []REAL x   = points[ :, 1 ];
                   []REAL y   = points[ :, 2 ];
                   FOR i TO 1 UPB points - 1 DO
                       result +:= x[ i     ] * y[ i + 1 ];
                       result -:= x[ i + 1 ] * y[ i     ]
                   OD;
                   result     +:= x[ n ] * y[ 1 ];
                   result     -:= x[ 1 ] * y[ n ]
               FI;
               ( ABS result ) / 2
           FI
       END # AREA # ;

   # test case as per the task #
   print( ( fixed( AREA [,]REAL( ( 3.0, 4.0 ), ( 5.0, 11.0 ), ( 12.0, 8.0 ), ( 9.0, 5.0 ), ( 5.0, 6.0 ) ), -6, 2 ), newline ) )

END </lang>

Output:

 30.00

C

Reads the points from a file whose name is supplied via the command line, prints out usage if invoked incorrectly. <lang C> /*Abhishek Ghosh, 25th September 2017*/

include<stdlib.h>
include<stdio.h>
include<math.h>

typedef struct{ double x,y; }point;

double shoelace(char* inputFile){ int i,numPoints; double leftSum = 0,rightSum = 0;

point* pointSet; FILE* fp = fopen(inputFile,"r");

fscanf(fp,"%d",&numPoints);

pointSet = (point*)malloc((numPoints + 1)*sizeof(point));

for(i=0;i<numPoints;i++){ fscanf(fp,"%lf %lf",&pointSet[i].x,&pointSet[i].y); }

fclose(fp);

pointSet[numPoints] = pointSet[0];

for(i=0;i<numPoints;i++){ leftSum += pointSet[i].x*pointSet[i+1].y; rightSum += pointSet[i+1].x*pointSet[i].y; }

free(pointSet);

return 0.5*fabs(leftSum - rightSum); }

int main(int argC,char* argV[]) { if(argC==1) printf("\nUsage : %s <full path of polygon vertices file>",argV[0]);

else printf("The polygon area is %lf square units.",shoelace(argV[1]));

return 0; } </lang> Input file, first line specifies number of points followed by the ordered vertices set with one vertex on each line.

Invocation and output :

C:\rosettaCode>shoelace.exe polyData.txt
The polygon area is 30.000000 square units.

Except for the use of "END FUNCTION name instead of just END, and the convenient function SUM with array span expressions (so SUM(P) rather than a DO-loop to sum the elements of array P), both standardised with F90, this would be acceptable to F66, which introduced complex number arithmetic. Otherwise, separate X and Y arrays would be needed, but complex numbers seemed convenient seeing as (x,y) pairs are involved. But because the MODULE facility of F90 has not been used, routines invoking functions must declare the type of the function names, especially if the defult types are unsuitable, as here. <lang Fortran> DOUBLE PRECISION FUNCTION AREA(N,P) !Calculates the area enclosed by the polygon P. C Uses the mid-point rule for integration. Consider the line joining (x1,y1) to (x2,y2) C The area under that line (down to the x-axis) is the midpoint (y1 + y2)/2 times the width (x2 - x1). C Now consider a sequence of such points heading in the +x direction: their areas are all positive. C Follow with a sequence of points heading in the -x direction, back to the first point: their areas are all negative. C The resulting sum is the area below the +x sequence and above the -x sequence: the area of the polygon. C The point sequence can wobble as it wishes and can meet the other side, but it must not cross itself c as would be done in a figure 8 drawn with a crossover instead of a meeting. C A clockwise traversal (as for an island) gives a positive area; use anti-clockwise for a lake.

      INTEGER N		!The number of points.
      DOUBLE COMPLEX P(N)	!The points.
      DOUBLE COMPLEX PP,PC	!Point Previous and Point Current.
      DOUBLE COMPLEX W		!Polygon centre. Map coordinates usually have large offsets.
      DOUBLE PRECISION A	!The area accumulator.
      INTEGER I		!A stepper.
       IF (N.LT.3) STOP "Area: at least three points are needed!"	!Good grief.
       W = (P(1) + P(N/3) + P(2*N/3))/3	!An initial working average.
       W = SUM(P(1:N) - W)/N + W	!A good working average is the average itself.
       A = 0			!The area enclosed by the point sequence.
       PC = P(N) - W		!The last point is implicitly joined to the first.
       DO I = 1,N		!Step through the positions.
         PP = PC			!Previous position.
         PC = P(I) - W			!Current position.
         A = (DIMAG(PC) + DIMAG(PP))*DBLE(PC - PP) + A	!Area integral component.
       END DO			!On to the next position.
       AREA = A/2		!Divide by two once.
     END FUNCTION AREA		!The units are those of the points.

     DOUBLE PRECISION FUNCTION AREASL(N,P)	!Area enclosed by polygon P, by the "shoestring" method.
      INTEGER N		!The number of points.
      DOUBLE COMPLEX P(N)	!The points.
      DOUBLE PRECISION A	!A scratchpad.
       A = SUM(DBLE(P(1:N - 1)*DIMAG(P(2:N)))) + DBLE(P(N))*DIMAG(P(1))
    1    - SUM(DBLE(P(2:N)*DIMAG(P(1:N - 1)))) - DBLE(P(1))*DIMAG(P(N))
       AREASL = A/2		!The midpoint formula requires a halving.
     END FUNCTION AREASL	!Negative for clockwise, positive for anti-clockwise.

     INTEGER ENUFF
     DOUBLE PRECISION AREA,AREASL	!The default types are not correct.
     DOUBLE PRECISION A1,A2		!Scratchpads, in case of a debugging WRITE within the functions.
     PARAMETER (ENUFF = 5)		!The specification.
     DOUBLE COMPLEX POINT(ENUFF)	!Could use X and Y arrays instead.
     DATA POINT/(3D0,4D0),(5D0,11D0),(12D0,8D0),(9D0,5D0),(5D0,6D0)/	!"D" for double precision.

     WRITE (6,*) POINT
     A1 = AREA(5,POINT)
     A2 = AREASL(5,POINT)
     WRITE (6,*) "A=",A1,A2
     END</lang>

Output: WRITE (6,*) means write to output unit six (standard output) with free-format (the *). Note the different sign convention.

 (3.00000000000000,4.00000000000000) (5.00000000000000,11.0000000000000)
 (12.0000000000000,8.00000000000000) (9.00000000000000,5.00000000000000)
 (5.00000000000000,6.00000000000000)
 A=   30.0000000000000       -30.0000000000000

The "shoestring" method came as a surprise to me, as I've always used what I had thought the "obvious" method. Note that function AREA makes one pass through the point data not two, and because map coordinate values often have large offsets a working average is used to reduce the loss of precision. This requires faith that SUM(P(1:N) - W) will be evaluated as written, not as SUM(P(1:N)) - N*W with even greater optimisation opportunity awaiting in cancelling further components of the expression. For example, the New Zealand metric grid has (2510000,6023150) as (Easting,Northing) or (x,y) at its central point of 41°S 173°E rather than (0,0) so seven digits of precision are used up. If anyone wants a copy of a set of point sequences for NZ (30,000 positions, 570KB) with lots of islands and lakes, even a pond in an island in a lake in the North Island...

FreeBASIC

<lang freebasic>' version 18-08-2017 ' compile with: fbc -s console

Type _point_

   As Double x, y

End Type

Function shoelace_formula(p() As _point_ ) As Double

   Dim As UInteger i
   Dim As Double sum

   For i = 1 To UBound(p) -1
       sum += p(i   ).x * p(i +1).y
       sum -= p(i +1).x * p(i   ).y
   Next
   sum += p(i).x * p(1).y
   sum -= p(1).x * p(i).y

   Return Abs(sum) / 2

End Function

' ------=< MAIN >=------

Dim As _point_ p_array(1 To ...) = {(3,4), (5,11), (12,8), (9,5), (5,6)}

Print "The area of the polygon ="; shoelace_formula(p_array())

' empty keyboard buffer While Inkey <> "" : Wend Print : Print "hit any key to end program" Sleep End</lang>

Output:

The area of the polygon = 30

J

Implementation:

<lang J>shoelace=:verb define

 0.5*|+/((* 1&|.)/ - (* _1&|.)/)|:y

)</lang>

Task example:

<lang J> shoelace 3 4,5 11,12 8,9 5,:5 6 30</lang>

Exposition:

We start with our list of coordinate pairs

<lang J> 3 4,5 11,12 8,9 5,:5 6

3  4
5 11

12 8

9  5
5  6</lang>

But the first thing we do is transpose them so that x coordinates and y coordinates are the two items we are working with:

We want to rotate the y list by one (in each direction) and multiply the x list items by the corresponding y list items. Something like this, for example:

Or, rephrased:

We'll be subtracting what we get when we rotate in the other direction, which looks like this:

Finally, we add up that list, take the absolute value (there are contexts where signed area is interesting - for example, some graphics application - but that was not a part of this task) and divide that by 2.

Julia

Works with: Julia version 0.6

Translation of: Python

<lang julia>""" Assumes x,y points go around the polygon in one direction. """ shoelacearea(x, y) =

   abs(sum(i * j for (i, j) in zip(x, append!(y[2:end], y[1]))) -
       sum(i * j for (i, j) in zip(append!(x[2:end], x[1]), y))) / 2

x, y = [3, 5, 12, 9, 5], [4, 11, 8, 5, 6] @show x y shoelacearea(x, y)</lang>

Output:

x = [3, 5, 12, 9, 5]
y = [4, 11, 8, 5, 6]
shoelacearea(x, y) = 30.0

Kotlin

<lang scala>// version 1.1.3

class Point(val x: Int, val y: Int) {

   override fun toString() = "($x, $y)"

}

fun shoelaceArea(v: List<Point>): Double {

   val n = v.size
   var a = 0.0
   for (i in 0 until n - 1) { 
       a += v[i].x * v[i + 1].y - v[i + 1].x * v[i].y
   }
   return Math.abs(a + v[n - 1].x * v[0].y - v[0].x * v[n -1].y) / 2.0

}

fun main(args: Array<String>) {

   val v = listOf(
       Point(3, 4), Point(5, 11), Point(12, 8), Point(9, 5), Point(5, 6)
   )
   val area = shoelaceArea(v) 
   println("Given a polygon with vertices at $v,")
   println("its area is $area")

}</lang>

Output:

Given a polygon with vertices at [(3, 4), (5, 11), (12, 8), (9, 5), (5, 6)],
its area is 30.0

Lua

<lang lua>function shoeArea(ps)

 local function det2(i,j)
   return ps[i][1]*ps[j][2]-ps[j][1]*ps[i][2]
 end
 local sum = #ps>2 and det2(#ps,1) or 0
 for i=1,#ps-1 do sum = sum + det2(i,i+1)end
 return math.abs(0.5 * sum)

end</lang> Using an accumulator helper inner function <lang lua>function shoeArea(ps)

 local function ssum(acc, p1, p2, ...)
   if not p2 or not p1 then
     return math.abs(0.5 * acc)
   else
     return ssum(acc + p1[1]*p2[2]-p1[2]*p2[1], p2, ...)
   end
 end
 return ssum(0, ps[#ps], table.unpack(ps))

end

local p = {{3,4}, {5,11}, {12,8}, {9,5}, {5,6}} print(shoeArea(p))-- 30 </lang> both version handle special cases of less than 3 point as 0 area result.

Mathematica

Geometry objects built-in in the Wolfram Language <lang Mathematica>Area[Polygon[{{3, 4}, {5, 11}, {12, 8}, {9, 5}, {5, 6}}]]</lang>

Output:

Perl 6

Index and mod offset

Works with: Rakudo version 2017.07

<lang perl6>sub area-by-shoelace(@p) {

   (^@p).map({@p[$_;0] * @p[($_+1)%@p;1] - @p[$_;1] * @p[($_+1)%@p;0]}).sum.abs / 2

}

say area-by-shoelace( [ (3,4), (5,11), (12,8), (9,5), (5,6) ] );</lang>

Output:

Slice and rotation

Works with: Rakudo version 2017.07

<lang perl6>sub area-by-shoelace ( @p ) {

   my @x := @p».[0];
   my @y := @p».[1];

   my $s := ( @x Z* @y.rotate( 1) ).sum
          - ( @x Z* @y.rotate(-1) ).sum;

   return $s.abs / 2;

}

say area-by-shoelace( [ (3,4), (5,11), (12,8), (9,5), (5,6) ] ); </lang>

Output:

Python

<lang python>>>> def area_by_shoelace(x, y):

   "Assumes x,y points go around the polygon in one direction"
   return abs( sum(i * j for i, j in zip(x,             y[1:] + y[:1]))
              -sum(i * j for i, j in zip(x[1:] + x[:1], y            ))) / 2

>>> points = [(3,4), (5,11), (12,8), (9,5), (5,6)] >>> x, y = zip(*points) >>> area_by_shoelace(x, y) 30.0 >>> </lang>

Racket

<lang racket>#lang racket/base

(struct P (x y))

(define (area . Ps)

 (define (A P-a P-b)
   (+ (for/sum ((p_i Ps)
                (p_i+1 (in-sequences (cdr Ps)
                                     (in-value (car Ps)))))
        (* (P-a p_i) (P-b p_i+1)))))
 (/ (abs (- (A P-x P-y) (A P-y P-x))) 2))

(module+ main

 (area (P 3 4) (P 5 11) (P 12 8) (P 9 5) (P 5 6)))</lang>

Output:

REXX

endpoints as exceptions

<lang rexx>/*REXX program uses a Shoelace formula to calculate the area of an N-sided polygon. */ parse arg pts /*obtain optional arguments from the CL*/ if pts= then pts= '(3,4),(5,11),(12,8),(9,5),(5,6)' /*Not specified? Use default. */ pts=space(pts, 0); @=pts /*elide extra blanks; save pts.*/

          do n=1  until @=                           /*perform destructive parse on @*/
          parse var @  '('  x.n  ","  y.n  ')'  ","  @ /*obtain  X  and  Y  coördinates*/
          end   /*n*/

A=0 /*initialize the area to zero.*/

          do j=1  for n;  jp=j+1;  if jp>n   then jp=1 /*adjust for  J  for overflow.  */
                          jm=j-1;  if jm==0  then jm=n /*   "    "   "   "  underflow. */
          A=A + x.j * (y.jp - y.jm)                    /*compute a part of the area.   */
          end   /*j*/

A=abs(A/2) /*obtain half of the │ A │ sum*/ say 'polygon area of ' n " points: " pts ' is ───► ' A /*stick a fork in it, we're done*/</lang>

output when using the default input:

polygon area of  5  points:  (3,4),(5,11),(12,8),(9,5),(5,6)   is ───►  30

endpoints as wrap-around

This REXX version uses a different method to define the X₀, Y₀, and X_n+1, Y_n+1 data points (and not treat them as exceptions).

When calculating the area for many polygons (or where the number of polygon sides is large), this method would be faster. <lang rexx>/*REXX program uses a Shoelace formula to calculate the area of an N-sided polygon. */ parse arg pts /*obtain optional arguments from the CL*/ if pts= then pts= '(3,4),(5,11),(12,8),(9,5),(5,6)' /*Not specified? Use default. */ pts=space(pts, 0); @=pts /*elide extra blanks; save pts.*/

          do n=1  until @=                           /*perform destructive parse on @*/
          parse var @  '('  x.n  ","  y.n  ')'  ","  @ /*obtain  X  and  Y  coördinates*/
          end   /*n*/

np=n+1 /*a variable to hold N+1 value*/ parse value x.1 y.1 x.n y.n with x.np y.np x.0 y.0 /*define X.0 & X.n+1 points.*/ A=0 /*initialize the area to zero.*/

          do j=1  for n;             jp=j+1;   jm=j-1  /*adjust for  J  for overflow.  */
          A=A + x.j * (y.jp - y.jm)                    /*compute a part of the area.   */
          end   /*j*/

A=abs(A/2) /*obtain half of the │ A │ sum*/ say 'polygon area of ' n " points: " pts ' is ───► ' A /*stick a fork in it, we're done*/</lang>

output is the same as the 1^st REXX version.

Scala

<lang scala>case class Point( x:Int,y:Int ) { override def toString = "(" + x + "," + y + ")" }

case class Polygon( pp:List[Point] ) {

 require( pp.size > 2, "A Polygon must consist of more than two points" )

 override def toString = "Polygon(" + pp.mkString(" ", ", ", " ") + ")"
 
 def area = {
 
   // Calculate using the Shoelace Formula
   val xx = pp.map( p => p.x )
   val yy = pp.map( p => p.y )
   val overlace = xx zip yy.drop(1)++yy.take(1)
   val underlace = yy zip xx.drop(1)++xx.take(1)
   
   (overlace.map( t => t._1 * t._2 ).sum - underlace.map( t => t._1 * t._2 ).sum).abs / 2.0
 }

}

// A little test... { val p = Polygon( List( Point(3,4), Point(5,11), Point(12,8), Point(9,5), Point(5,6) ) )

assert( p.area == 30.0 )

println( "Area of " + p + " = " + p.area ) } </lang>

Output:

Area of Polygon( (3,4), (5,11), (12,8), (9,5), (5,6) ) = 30.0

Sidef

Translation of: Perl 6

<lang ruby>func area_by_shoelace (*p) {

   var x = p.map{_[0]}
   var y = p.map{_[1]}

   var s = (
       (x ~Z* y.rotate(+1)).sum -
       (x ~Z* y.rotate(-1)).sum
   )

   s.abs / 2

}

say area_by_shoelace([3,4], [5,11], [12,8], [9,5], [5,6])</lang>

Output:

zkl

By the "book": <lang zkl>fcn areaByShoelace(points){ // ( (x,y),(x,y)...)

  xs,ys:=Utils.Helpers.listUnzip(points); // (x,x,...), (y,y,,,)
  ( xs.zipWith('*,ys[1,*]).sum(0) + xs[-1]*ys[0] - 
    xs[1,*].zipWith('*,ys).sum(0) - xs[0]*ys[-1] )
  .abs().toFloat()/2;

}</lang> or an iterative solution: <lang zkl>fcn areaByShoelace2(points){ // ( (x,y),(x,y)...)

  xs,ys:=Utils.Helpers.listUnzip(points); // (x,x,...), (y,y,,,)
  N:=points.len();
  N.reduce('wrap(s,n){ s + xs[n]*ys[(n+1)%N] - xs[(n+1)%N]*ys[n] },0)
  .abs().toFloat()/2;

}</lang> <lang zkl>points:=T(T(3,4), T(5,11), T(12,8), T(9,5), T(5,6)); areaByShoelace(points).println(); areaByShoelace2(points).println();</lang>

Output:

30
30