(49g 50g) Shoelace algorithm

08232018, 02:20 PM
(This post was last modified: 08262018 09:23 PM by John Keith.)
Post: #1




(49g 50g) Shoelace algorithm
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 twoelement vectors. Code:
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:
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. 

08232018, 03:33 PM
Post: #2




RE: ( HP49/50) Shoelace algorithm
(08232018 02:20 PM)John Keith Wrote: My program differs in that it takes a list of coordinates of the form Why is SWAP necessary ? Does ABS remove the sign of sum anyway ? Sorry if the question sound stupid, I do not know RPL. A minitutorial 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 ? 

08232018, 08:23 PM
(This post was last modified: 08242018 02:15 PM by John Keith.)
Post: #3




RE: ( HP49/50) Shoelace algorithm
(08232018 03:33 PM)Albert Chan Wrote: Why is SWAP necessary ? 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 speedup is due to the math being done on the stack rather than by higher level commands such as CROSS. That is one of the tradeoffs 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 

08232018, 09:49 PM
Post: #4




RE: ( HP49/50) Shoelace algorithm
(08232018 08:23 PM)John Keith Wrote: The SWAP is necessary for each group of coordinates because the answer depends on the order of subtraction. 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, ...])) 

08242018, 03:56 AM
Post: #5




RE: ( HP49/50) Shoelace algorithm
Alternatively we can use
\(A=\frac{1}{2}\sum_{i=1}^{n}x_{i}(y_{i+1}y_{i1})\) where \(y_{n+1}=y_{1}\) and \(y_{0}=y_{n}\). (xs ys  area) Code: « DUP TAIL OVER HEAD + SWAP 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 

08242018, 01:40 PM
Post: #6




RE: ( HP49/50) Shoelace algorithm
(08232018 09:49 PM)Albert Chan Wrote: I thought "UNROT * UNROT *" leave only 2 values on stack. Apparently you are correct! I tried several sets of coordinates and got the correct answer without the SWAP each time. 

08242018, 02:08 PM
(This post was last modified: 08242018 02:09 PM by John Keith.)
Post: #7




RE: ( HP49/50) Shoelace algorithm
(08242018 03:56 AM)Thomas Klemm Wrote: Alternatively we can use 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 areaandcentroid program I'm all ears. John 

08242018, 03:27 PM
(This post was last modified: 08302018 10:05 PM by Albert Chan.)
Post: #8




RE: ( HP49/50) Shoelace algorithm
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 

08242018, 05:02 PM
Post: #9




RE: ( HP49/50) Shoelace algorithm
(08242018 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 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 

08242018, 07:57 PM
Post: #10




RE: ( HP49/50) Shoelace algorithm
(08242018 05:02 PM)Thomas Klemm Wrote:(08242018 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. 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 builtin command is inexplicably slower than a User RPL equivalent. 

08252018, 01:51 AM
Post: #11




RE: ( HP49/50) Shoelace algorithm
(08242018 02:08 PM)John Keith Wrote: If you have any suggestions for improving the rest of my areaandcentroid program I'm all ears. I used these formulas found in Centroid: \(C_{\mathrm {x} }={\frac {1}{6A}}\sum _{i=0}^{n1}(x_{i}+x_{i+1})(x_{i}\ y_{i+1}x_{i+1}\ y_{i})\) \(C_{\mathrm {y} }={\frac {1}{6A}}\sum _{i=0}^{n1}(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 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 

08252018, 01:58 PM
(This post was last modified: 08252018 01:58 PM by John Keith.)
Post: #12




RE: ( HP49/50) Shoelace algorithm
(08242018 05:02 PM)Thomas Klemm Wrote: Does that mean that we could use this? Yes, it not only works, it is the smallest and fastest program. The complete list, timing in milliseconds for the Cat: Code:


08252018, 02:02 PM
(This post was last modified: 08252018 03:55 PM by Albert Chan.)
Post: #13




RE: ( HP49/50) Shoelace algorithm
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 = rotatedleft x, yr = rotatedleft 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 

08252018, 03:17 PM
Post: #14




RE: ( HP49/50) Shoelace algorithm  
08252018, 03:27 PM
Post: #15




RE: ( HP49/50) Shoelace algorithm
(08252018 01:51 AM)Thomas Klemm Wrote: I used these formulas found in Centroid: 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:
To convert in the opposite direction: Code:


08252018, 03:28 PM
Post: #16




RE: ( HP49/50) Shoelace algorithm
(08252018 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. 

08252018, 03:57 PM
Post: #17




RE: ( HP49/50) Shoelace algorithm  
08262018, 01:47 PM
(This post was last modified: 08272018 06:30 PM by John Keith.)
Post: #18




RE: (49g 50g) Shoelace algorithm
The only improvement I can see for Thomas's areaandcentroid 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 nonintersecting 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 

08272018, 05:20 PM
Post: #19




RE: (49g 50g) Shoelace algorithm
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 datastructure for the implementation of the algorithm. So it feels a bit like cheating. Best regards Thomas 

02282019, 04:26 AM
Post: #20




RE: (49g 50g) Shoelace algorithm
(08242018 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 I've read somewhere that it's faster to do the rotation on the stack. Cheers Thomas 

« Next Oldest  Next Newest »

User(s) browsing this thread: 1 Guest(s)