(49g 50g) Shoelace algorithm

[John Keith]

EDIT: These programs have been largely obsoleted by Thomas Klemm's programs in post numbers 5 and 11
respectively. I am leaving them here for illustrative purposes only.

---

Inspired by Thomas Klemm's elegant RPL implementation of the shoelace method for calculating the area of a polygon, I am listing my version here. Though it is longer and not nearly as elegant, it is over 2.5 times faster.

My program differs in that it takes a list of coordinates of the form
{ X1 Y1 X2 Y2...Xn Yn }
instead of a list of two-element vectors.

Code:


\<< DUP 1. 2. SUB + 4.
  \<< NSUB 2. MOD { UNROT * UNROT * SWAP - } { DROP2 DROP2 } IFTE
  \>> DOSUBS \GSLIST ABS 2 /
\>>

The following extended version calculates the area and centroid (aka barycenter) of a polygon. The results are returned as a list of the form
{ Area X Y }.

Code:


\<< DUP 1. 2. SUB + DUP 4.
  \<< NSUB 2. MOD { UNROT * UNROT * SWAP - } { DROP2 DROP2 } IFTE
  \>> DOSUBS OVER 4.
  \<< NSUB 2. MOD { DROP NIP + } { DROP2 DROP2 } IFTE
  \>> DOSUBS OVER * \GSLIST SWAP ROT 4.
  \<< NSUB 2. MOD { NIP ROT DROP + } { DROP2 DROP2 } IFTE
  \>> DOSUBS OVER * \GSLIST SWAP \GSLIST 2 / DUP ABS EVAL 4. ROLLD 6 * INV DUP ROT * EVAL UNROT * EVAL SWAP 3. \->LIST
\>>

Both programs will return exact integer or rational results if the input coordinates are exact and the program is run in exact mode. If this is not desirable, the integer(s) in the last lines of the programs can be changed to reals, and the EVALs in the last line of the second program can be eliminated. This will make the programs run slightly faster.

[Albert Chan]

(08-23-2018 02:20 PM)John Keith Wrote:  My program differs in that it takes a list of coordinates of the form
{ X1 Y1 X2 Y2...Xn Yn }
instead of a list of two-element vectors.

Code:


\<< DUP 1. 2. SUB + 4.
  \<< NSUB 2. MOD { UNROT * UNROT * SWAP - } { DROP2 DROP2 } IFTE
  \>> DOSUBS \GSLIST ABS 2 /
\>>

Why is SWAP necessary ?
Does ABS remove the sign of sum anyway ?

Sorry if the question sound stupid, I do not know RPL.
A mini-tutorial of how above code work would be nice ...

Does "DUP 1. 2. SUB +" produce {X1 Y1 X2 Y2...Xn Yn X1 Y1} ?
Is this "flatten" list (pick 4 at a time) the reason for the 2.5X speedup ?

[John Keith]

(08-23-2018 03:33 PM)Albert Chan Wrote:  Why is SWAP necessary ?
Does ABS remove the sign of sum anyway ?

Sorry if the question sound stupid, I do not know RPL.
A mini-tutorial of how above code work would be nice ...

Does "DUP 1. 2. SUB +" produce {X1 Y1 X2 Y2...Xn Yn X1 Y1} ?
Is this "flatten" list (pick 4 at a time) the reason for the 2.5X speedup ?

To answer your questions more or less in order:

The SWAP is necessary for each group of coordinates because the answer depends on the order of subtraction. The ABS at end gives a positive answer even if the points are ordered clockwise. EDIT: apparently not, see post below.

The code works exactly the same way as Thomas's code in the linked thread. Most of the speed-up is due to the math being done on the stack rather than by higher level commands such as CROSS. That is one of the trade-offs in RPL even more than in other languages. Stack commands are less readable but much faster.

And yes, the DUP 1. 2. SUB + appends the first pair of coordinates onto the end of the list as the algorithm requires.

John

[Albert Chan]

(08-23-2018 08:23 PM)John Keith Wrote:  The SWAP is necessary for each group of coordinates because the answer depends on the order of subtraction.
The ABS at end gives a positive answer even if the points are ordered clockwise.

I thought "UNROT * UNROT *" leave only 2 values on stack.
If true, "SWAP -" is the same as "-" with opposite sign.

Opposite sign should not matter:

2*area = abs(sum([x1*y2 - x2*y1, x2*y3 - x3*y2, ...])) = abs(sum([x2*y1 - x1*y2, x3*y2 - x2*y3, ...]))

[Thomas Klemm]

Alternatively we can use

\(A=\frac{1}{2}\sum_{i=1}^{n}x_{i}(y_{i+1}-y_{i-1})\)

where \(y_{n+1}=y_{1}\) and \(y_{0}=y_{n}\).

(xs ys -- area)

Code:

« DUP TAIL OVER HEAD + SWAP
  REVLIST DUP TAIL SWAP HEAD + REVLIST
  - * ∑LIST 2 /
»

The coordinates of the polygon have to be entered as two separate lists xs and ys.
There might be better ways to rotate the elements of ys left and right.

Kind regards
Thomas

[John Keith]

(08-23-2018 09:49 PM)Albert Chan Wrote:  I thought "UNROT * UNROT *" leave only 2 values on stack.
If true, "SWAP -" is the same as "-" with opposite sign.

Opposite sign should not matter:

2*area = abs(sum([x1*y2 - x2*y1, x2*y3 - x3*y2, ...])) = abs(sum([x2*y1 - x1*y2, x3*y2 - x2*y3, ...]))

Apparently you are correct! I tried several sets of coordinates and got the correct answer without the SWAP each time.

[John Keith]

(08-24-2018 03:56 AM)Thomas Klemm Wrote:  Alternatively we can use

\(A=\frac{1}{2}\sum_{i=1}^{n}x_{i}(y_{i+1}-y_{i-1})\)

where \(y_{n+1}=y_{1}\) and \(y_{0}=y_{n}\).

(xs ys -- area)

Code:

« DUP TAIL OVER HEAD + SWAP
  REVLIST DUP TAIL SWAP HEAD + REVLIST
  - * ∑LIST 2 /
»

The coordinates of the polygon have to be entered as two separate lists xs and ys.
There might be better ways to rotate the elements of ys left and right.

Kind regards
Thomas

Excellent program, smaller and faster than mine. It also runs on the 48g.

The ListExt Library command LRLLD is the equivalent of DUP TAIL SWAP HEAD + but ListExt is HP49/50 only.

HEAD is very slow and I try to avoid it whenever possible. Though much longer, 1. DUP SUB EVAL is about three times as fast as HEAD!

If you have any suggestions for improving the rest of my area-and-centroid program I'm all ears.

John

[Albert Chan]

Is this a valid RPL program ?

\(A=\frac{1}{2}|\sum_{i=1}^{n}(x_{i+1}+x_{i})(y_{i+1}-y_{i})|\)

where \(x_{n+1}=x_{1}\) and \(y_{n+1}=y_{1}\).

(xs ys -- area)

Code:

« DUP HEAD + ΔLIST               @ do y diff
  SWAP DUP DUP HEAD + TAIL ADD   @ do x sums
  * ∑LIST ABS 2 /                @ do area
»

[Thomas Klemm]

(08-24-2018 02:08 PM)John Keith Wrote:  The ListExt Library command LRLLD is the equivalent of DUP TAIL SWAP HEAD + but ListExt is HP49/50 only.

Does that mean that we could use this?

Code:

« DUP LRLLD
  SWAP LROLL
  - * ∑LIST 2 /
»

Quote:HEAD is very slow and I try to avoid it whenever possible.

Do we know the reason? With a linked list HEAD and TAIL are usually fast.

Quote:Though much longer, 1. DUP SUB EVAL is about three times as fast as HEAD!

What about 1. GET?

Kind regards
Thomas

[John Keith]

(08-24-2018 05:02 PM)Thomas Klemm Wrote:  

(08-24-2018 02:08 PM)John Keith Wrote:  The ListExt Library command LRLLD is the equivalent of DUP TAIL SWAP HEAD + but ListExt is HP49/50 only.

Does that mean that we could use this?

Code:

« DUP LRLLD
  SWAP LROLL
  - * ∑LIST 2 /
»

Quote:HEAD is very slow and I try to avoid it whenever possible.

Do we know the reason? With a linked list HEAD and TAIL are usually fast.

Quote:Though much longer, 1. DUP SUB EVAL is about three times as fast as HEAD!

What about 1. GET?

Kind regards
Thomas

I won't have time to try the ListExt code until tomorrow but it looks about right.

I have no idea why HEAD is so slow, maybe some folks here with more System RPL knowledge can tell us.

1 GET is somewhere in between HEAD and 1. DUP SUB EVAL. I don't know the ratio off hand. I was using 1. DUP SUB EVAL mainly to illustrate how a built-in command is inexplicably slower than a User RPL equivalent.

[Thomas Klemm]

(08-24-2018 02:08 PM)John Keith Wrote:  If you have any suggestions for improving the rest of my area-and-centroid program I'm all ears.

I used these formulas found in Centroid:

\(C_{\mathrm {x} }={\frac {1}{6A}}\sum _{i=0}^{n-1}(x_{i}+x_{i+1})(x_{i}\ y_{i+1}-x_{i+1}\ y_{i})\)

\(C_{\mathrm {y} }={\frac {1}{6A}}\sum _{i=0}^{n-1}(y_{i}+y_{i+1})(x_{i}\ y_{i+1}-x_{i+1}\ y_{i})\)

where \(x_{n}=x_{0}\) and \(y_{n}=y_{0}\).

(xs ys -- x y area)

Code:

« DUP TAIL OVER HEAD + ROT
  DUP TAIL OVER HEAD +
  → y v x u
  « x v * u y * -
    x u ADD OVER * ∑LIST
    OVER y v ADD * ∑LIST
    ROT ∑LIST
  »
  ROT OVER 3 * / SWAP
  ROT OVER 3 * / SWAP
  2 /
»

Example

For the cat from this video:

xs: { 4 0 -2 -6 -1 5 }
ys: { 4 1 5 0 -4 -2 }

The result is:

3:       -.169696969697
2:        .030303030303
1:                   55

Cheers
Thomas

[John Keith]

(08-24-2018 05:02 PM)Thomas Klemm Wrote:  Does that mean that we could use this?

Code:

« DUP LRLLD
  SWAP LROLL
  - * ∑LIST 2 /
»

Yes, it not only works, it is the smallest and fastest program.

The complete list, timing in milliseconds for the Cat:

Code:


Your vector version: 446ms  57 bytes
My program above:    164ms 107    "
Your program above:  157ms  76    "
ListExt version:     140ms  41.5  "

[Albert Chan]

Hi, Thomas Klemm. Nice post.

Here is another way, calculations can be done in parallel.
If only x's or y's changed, only 1 side need the update.

Let xr = rotated-left x, yr = rotated-left y. Borrowing your example:
x = [4, 0, -2, -6, -1, 5]
xr= [0, -2, -6, -1, 5, 4]
y = [4, 1, 5, 0, -4, -2]
yr= [1, 5, 0, -4, -2, 4]

2*A = sum((x*yr - xr*y))
= sum(xr*yr - xr*y + x*yr - x*y)
= sum((xr + x)(yr - y))

6*A*Cx = sum( (xr + x) (x*yr - xr*y) )
= sum((xr + x) * ((xr + x)(yr - y) - (xr*yr - x*y)))
= sum((xr + x)(xr + x)(yr - y) - (xr^2*yr - x*xr*y + x*xr*yr - x^2*y))
= sum(((xr + x)^2 - x*xr) (yr - y))

Confirm the numbers ...

xr + x = [4, -2, -8, -7, 4, 9]
yr - y = [-3, 4, -5, -4, 2, 6]
A = 1/2 sum([-12, -8, 40, 28, 8, 54]) = 110/2 = 55

(xr + x)^2 = [16, 4, 64, 49, 16, 81]
x * xr = [0, 0, 12, 6, -5, 20]
xdiff = [16, 4, 52, 43, 21, 61]
6A*Cx = sum(xdiff * (yr - y)) = sum([-48, 16, -260, -172, 42, 366]) = -56
Cx = -56/6/55 = -28/165 ~ -0.1696969697

(yr + y)^2 = [25, 36, 25, 16, 36, 4]
y * yr = [ 4, 5, 0, 0, 8, -8]
ydiff = [21, 31, 25, 16, 28, 12]
xr - x = [-4, -2, -4, 5, 6, -1]
6A*Cy = sum(ydiff * (xr - x)) = sum([-84, -62, -100, 80, 168, -12]) = -10
Cy = -10/6/55 = -1/33 ~ -0.0303030303

[Thomas Klemm]

(08-25-2018 02:02 PM)Albert Chan Wrote:  2*A = sum((x*yr - xr*y))
= sum(xr^2 - xr*y + x*yr + x^2)
= sum((xr + x)(yr - y))

I assume that this should rather be:

2*A = sum(x*yr - xr*y)
= sum(xr*yr - xr*y + x*yr - x*y)
= sum((xr + x)(yr - y))

[John Keith]

(08-25-2018 01:51 AM)Thomas Klemm Wrote:  I used these formulas found in Centroid:

\(C_{\mathrm {x} }={\frac {1}{6A}}\sum _{i=0}^{n-1}(x_{i}+x_{i+1})(x_{i}\ y_{i+1}-x_{i+1}\ y_{i})\)

\(C_{\mathrm {y} }={\frac {1}{6A}}\sum _{i=0}^{n-1}(y_{i}+y_{i+1})(x_{i}\ y_{i+1}-x_{i+1}\ y_{i})\)

where \(x_{n}=x_{0}\) and \(y_{n}=y_{0}\).

(xs ys -- x y area)

Code:

« DUP TAIL OVER HEAD + ROT
  DUP TAIL OVER HEAD +
  → y v x u
  « x v * u y * -
    x u ADD OVER * ∑LIST
    OVER y v ADD * ∑LIST
    ROT ∑LIST
  »
  ROT OVER 3 * / SWAP
  ROT OVER 3 * / SWAP
  2 /
»

Example

For the cat from this video:

xs: { 4 0 -2 -6 -1 5 }
ys: { 4 1 5 0 -4 -2 }

The result is:

3:       -.169696969697
2:        .030303030303
1:                   55

Cheers
Thomas

Fantastic! Once again smaller and faster.

My program from Post 1: 534ms, 290 bytes
Your program: 441ms, 194 bytes.

I may try a ListExt version when I have some free time.

By the way, ListExt commands make it easy to convert between the flat list format used by my programs and the two separate lists used by yours.

To covert from { X1 Y1 X2 Y2...Xn Yn } to
{ X1 X2...Xn }
{ Y1 Y2...Yn } :

Code:


\<< DUP SIZE 2. / LDIST EVAL
\>>

To convert in the opposite direction:

Code:


\<< 2. \->LIST LCLLT
\>>

[Thomas Klemm]

(08-25-2018 01:58 PM)John Keith Wrote:  Yes, it not only works, it is the smallest and fastest program.

Thanks for taking the time to measure it!
It would be interesting to see the timings for larger examples than the cat.
Just to get a feeling for what happens if we increase the amount of points of the polygon.

Best regards
Thomas

---

A bit off topic but you might still enjoy Ian's Shoelace Site.

[Albert Chan]

(08-25-2018 03:17 PM)Thomas Klemm Wrote:  I assume that this should rather be:

2*A = sum(x*yr - xr*y)
= sum(xr*yr - xr*y + x*yr - x*y)
= sum((xr + x)(yr - y))

thanks. corrected.

[John Keith]

The only improvement I can see for Thomas's area-and-centroid program using ListExt commands is to simply replace TAIL OVER HEAD + with LRLLD so that the first two lines become
DUP LRLLD ROT DUP LRLLD.

This saves about 10ms for the Cat and reduces the size from 194 to 172.5 bytes, at the cost of HP48g compatibility.

Congratulations again to Thomas for a fine example of RPL programming!

I can't think of a way to make a non-intersecting polygon with a large number of points other than points spread evenly around a circle, or tediously typing hundreds of points in. Based on experience I would guess that the speed ratios we have seen would hold pretty well for larger polygons.

EDIT: I created a random (improper, multiply intersecting) polygon of 500 points, and a polygon made from 500 points evenly spread around a circle. For both polygons, my program from the first post takes about 42 seconds, and Thomas's program from post #11 takes about 18 seconds. The version using LRLLD is only about 100 milliseconds faster.

John

[Thomas Klemm]

John, thanks again for measuring the timings of these different implementations.

One of the reason for the difference could be that my solution uses two lists xs and ys while your solution uses a single list with the xs and ys interleaved. Thus you're forced to loop over twice the length of the list and throw half of them out.

That's when a 2nd parameter step for DOSUBS would be handy that controls how many elements to advance with each pass. In this case we would use step = 2.

And then I have no idea how DOSUBS behaves compared to ∑LIST and the other list operations ADD, – and *.

I've optimised the input data-structure for the implementation of the algorithm. So it feels a bit like cheating.

Best regards
Thomas

[Thomas Klemm]

(08-24-2018 03:56 AM)Thomas Klemm Wrote:  There might be better ways to rotate the elements of ys left and right.

(xs ys -- area)

Code:

« DUP SIZE → xs ys s
  « ys OBJ→ ROLLD s → LIST
    ys OBJ→ ROLL s →LIST
    - xs * ΣLIST 2 /
  »
»

I've read somewhere that it's faster to do the rotation on the stack.

Cheers
Thomas