+- HP Forums (https://www.hpmuseum.org/forum)
+-- Forum: HP Software Libraries (/forum-10.html)
+--- Forum: General Software Library (/forum-13.html)
+--- Thread: (49g 50g) Shoelace algorithm (/thread-11259.html)
Pages: 1 2
(49g 50g) Shoelace algorithm - John Keith - 08-23-2018 02:20 PM
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.
RE: ( HP49/50) Shoelace algorithm - Albert Chan - 08-23-2018 03:33 PM
(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 ?
RE: ( HP49/50) Shoelace algorithm - John Keith - 08-23-2018 08:23 PM
(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
RE: ( HP49/50) Shoelace algorithm - Albert Chan - 08-23-2018 09:49 PM
(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, ...]))
RE: ( HP49/50) Shoelace algorithm - Thomas Klemm - 08-24-2018 03:56 AM
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
RE: ( HP49/50) Shoelace algorithm - John Keith - 08-24-2018 01:40 PM
(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.
RE: ( HP49/50) Shoelace algorithm - John Keith - 08-24-2018 02:08 PM
(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
RE: ( HP49/50) Shoelace algorithm - Albert Chan - 08-24-2018 03:27 PM
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
»RE: ( HP49/50) Shoelace algorithm - Thomas Klemm - 08-24-2018 05:02 PM
(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
RE: ( HP49/50) Shoelace algorithm - John Keith - 08-24-2018 07:57 PM
(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.
RE: ( HP49/50) Shoelace algorithm - Thomas Klemm - 08-25-2018 01:51 AM
(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
RE: ( HP49/50) Shoelace algorithm - John Keith - 08-25-2018 01:58 PM
(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 "RE: ( HP49/50) Shoelace algorithm - Albert Chan - 08-25-2018 02:02 PM
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
RE: ( HP49/50) Shoelace algorithm - Thomas Klemm - 08-25-2018 03:17 PM
(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))
RE: ( HP49/50) Shoelace algorithm - John Keith - 08-25-2018 03:27 PM
(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
\>>RE: ( HP49/50) Shoelace algorithm - Thomas Klemm - 08-25-2018 03:28 PM
(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.
RE: ( HP49/50) Shoelace algorithm - Albert Chan - 08-25-2018 03:57 PM
(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.
RE: (49g 50g) Shoelace algorithm - John Keith - 08-26-2018 01:47 PM
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
RE: (49g 50g) Shoelace algorithm - Thomas Klemm - 08-27-2018 05:20 PM
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
RE: (49g 50g) Shoelace algorithm - Thomas Klemm - 02-28-2019 04:26 AM
(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