Shoelace formula for polygonal area: Difference between revisions

Content added Content deleted

Inline

Revision as of 10:56, 25 September 2017

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.

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.

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.

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

version 1

<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

version 2

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

@@ Line 330: / Line 330: @@
 ===version 2===
-This REXX version uses a different method to define the &nbsp; <big>'''X<sub>0</sub>, &nbsp; Y<sub>0</sub>'''</big>, &nbsp; and &nbsp; <big>'''X<sub>n+1</sub>''', &nbsp; <big>'''Y<sub>n+1</sub>'''</big> &nbsp; data points &nbsp; (and not treat them as exceptions).
+This REXX version uses a different method to define the &nbsp; <big>'''X<sub>0</sub>, &nbsp; Y<sub>0</sub>'''</big>, &nbsp; and &nbsp; <big>'''X<sub>n+1</sub>''', &nbsp; '''Y<sub>n+1</sub>'''</big> &nbsp; data points &nbsp; (and not treat them as exceptions).
 When calculating the area for many polygons &nbsp; (or where the number of polygon sides is large), &nbsp; this method would be faster.