libgeos
minimumrotatedrectangle not working properly
I found out that my own implementation (via spatialite) finds smaller area's for minimumrotatedrectangles when applying to all buildings in the Netherlands. Sometimes even 29% smaller. So significantly wrong. Only a few buildings are smaller via GEOS compared to my spatialite method, but the difference is 0.1% (to my opinion a precision difference).
Here a link to the building with 29% difference. http://bag.basisregistraties.overheid.nl/bag/id/pand/0502100000021153
I am willing to look at the implementation to see if I can find the bug, but miss the knowhow on where this function is implemented. Issue could also be part of the conversion of postgis to geom , since I am using postgis ST_orientedenvelope to perform this operation.
Any thoughts or insight is appreciated and with tips I am eager to help fixing this.
Can you provide WKT (and an image if possible) of the output from your implementation?
Here's the result from JTS (which should be the same as GEOS). It's obviously sub-optimal - a rectangle aligned to the corners of the geometry would have smaller area. In fact, even the geometry envelope has smaller area!
So now the question is: where did that algorithm come from, and is it basically wrong or is this a bug?
Here's a simple reproducer:
POLYGON ((150 300, 200 300, 300 300, 300 250, 280 120, 210 100, 100 100, 100 240, 150 300))
The problem is that the algorithm used is fundamentally wrong! It should be using the Rotating Calipers algorithm, but is using a simpler, incorrect approach.
We'll put this in the development queue.
Hi Martin,
Thanks for the response. Seeing your responses, I think an image of my results is no longer needed? My algorithm indeed finds the bounding box that is aligned to the corners of the ‘2 corner rounded square’.
My algorithm is similar to chapter 2’s algorithm of the rotating calipier document (The Exhaustive Search Algorithms), and therefor accurate but slow. It takes about 1 hour and 15 minutes for little over 6.1 million geometries. Although I think it can be made a little bit faster when I spend some time in optimizing my queries.
I think an image of my results is no longer needed?
Yes, no longer needed.
My algorithm is similar to chapter 2’s algorithm of the rotating calipier document (The Exhaustive Search Algorithms), and therefor accurate but slow. It takes about 1 hour and 15 minutes for little over 6.1 million geometries. Although I think it can be made a little bit faster when I spend some time in optimizing my queries.
Hard to say if that is "slow" - 6.1 million geometries is a lot! We'll see how a correct implementation in GEOS can do.
Hi Martin and others in need of a solution,
The query below does the correct implementation in finding the minimal bounding box. It is not as fast as using the st_orientedenvelope (with the error), but at least it returns the correct results. On my machine it took me about 29 minutes for 6.1 million objects/geometries, compared to little over 2 minutes for st_orientedenvelope (both including calculating the bounding box dimensions).
To apply the query on your data, replace “my_table” with your source table and “identification” and “geom” with the primary key and geometry column of your source table (all in bold).
Since I needed the measures of the boundingbox in a separate table containing measures as well, I am using the yellow marked parts to set the response as a new table and getting the dimensions.
Kind regards, Pim
/* the query / create table myBoundingBoxResultTable as with segments as ( SELECT distinct identificatie, st_azimuth(sp, ep)::numeric % (0.5 * pi())::numeric as angle, geom FROM ( SELECT identificatie, ST_PointN(geom, generate_series(1, ST_NPoints(geom)-1)) as sp, ST_PointN(geom, generate_series(2, ST_NPoints(geom) )) as ep, geom FROM ( SELECT identificatie as identificatie, st_exteriorring(st_convexhull(geom)) as geom FROM my_table ) AS ls ) AS seg ), envelopes as ( select identificatie, angle, st_envelope(st_rotate(geom, angle)) as envelope from segments ), selectedboundingbox as ( select distinct on (identificatie) identificatie, angle, envelope, st_area(envelope) as min_area from envelopes order by identificatie, st_area(envelope) asc ) select identificatie, st_rotate(envelope, -angle) as boundingbox, min_area, st_y(ST_PointN(ST_ExteriorRing(envelope), 2)) - st_y(ST_PointN(ST_ExteriorRing(envelope), 1)) AS dimension1, st_x(ST_PointN(ST_ExteriorRing(envelope), 3)) - st_x(ST_PointN(ST_ExteriorRing(envelope), 2)) AS dimension2 from selectedboundingbox /*********/
From: Martin Davis @.> Sent: woensdag 14 september 2022 00:56 To: libgeos/geos @.> Cc: Helm, P.W. (Pim) van den @.>; Author @.> Subject: Re: [libgeos/geos] minimumrotatedrectangle not working properly (Issue #679)
I think an image of my results is no longer needed?
Yes, no longer needed.
My algorithm is similar to chapter 2’s algorithm of the rotating calipier document (The Exhaustive Search Algorithms), and therefor accurate but slow. It takes about 1 hour and 15 minutes for little over 6.1 million geometries. Although I think it can be made a little bit faster when I spend some time in optimizing my queries.
Hard to say if that is "slow" - 6.1 million geometries is a lot! We'll see how a correct implementation in GEOS can do.
— Reply to this email directly, view it on GitHubhttps://github.com/libgeos/geos/issues/679#issuecomment-1246028891, or unsubscribehttps://github.com/notifications/unsubscribe-auth/AA4QLTBJVINWS2ULNVBPPSTV6EA7RANCNFSM6AAAAAAQLRH7WA. You are receiving this because you authored the thread.Message ID: @.***> This message may contain information that is not intended for you. If you are not the addressee or if this message was sent to you by mistake, you are requested to inform the sender and delete the message. TNO accepts no liability for the content of this e-mail, for the manner in which you use it and for damage of any kind resulting from the risks inherent to the electronic transmission of messages.
