Avoid `atan2` in `minimum_rotated_rect` #1094

lnicola · 2023-10-20T08:48:13Z

I agree to follow the project's code of conduct.
I added an entry to CHANGES.md if knowledge of this change could be valuable to users.

lnicola · 2023-10-20T08:49:21Z

No benchmarks yet, and the output might be less precise, not sure if we want to merge this, but I wanted to put the code out.

urschrei · 2023-10-20T13:52:15Z

geo/src/algorithm/minimum_rotated_rect.rs

+            let (mut min_x, mut max_x) = (<T as Bounded>::max_value(), <T as Bounded>::min_value());
+            let (mut min_y, mut max_y) = (<T as Bounded>::max_value(), <T as Bounded>::min_value());
+            for point in hull.exterior().points() {
+                let x = point.y() * edge.1 + point.x() * edge.0;


Herbie thinks you can improve the accuracy of these two statements using fma:

let x = point.x().mul_add(edge.0, (point.y() * edge.1)); let y = point.y().mul_add(edge.0, (-point.x() * edge.1));

Does it make any difference?

Ooh, good idea, forgot about FMA.

Tried it, slightly better but still no luck :(.

Looks like I can still extract one ULP or so by dividing by norm after the multiplications. Reverting FMA still doesn't fix those two vertices.

left: LineString([Coord { x: 1.7000000000000002, y: 24.6 }, Coord { x: 14.649854163628412, y: 25.15341257109523 }, Coord { x: 14.400000000000002, y: 31.000000000000004 }, Coord { x: 1.4501458363715911, y: 30.446587428904774 }, Coord { x: 1.7000000000000002, y: 24.6 }]) right: LineString([Coord { x: 1.7000000000000004, y: 24.600000000000005 }, Coord { x: 14.64985416362841, y: 25.153412571095235 }, Coord { x: 14.400000000000002, y: 31.000000000000004 }, Coord { x: 1.4501458363715916, y: 30.446587428904774 }, Coord { x: 1.7000000000000004, y: 24.600000000000005 }])

michaelkirk · 2023-10-27T08:25:49Z

geo/src/algorithm/minimum_rotated_rect.rs

-                (14.649854163628408, 25.153412571095235),
-                (1.7000000000000006, 24.6),
+                (1.7000000000000004, 24.600000000000005),
+                (14.64985416362841, 25.153412571095235),


Maybe you noticed this but it looks like the y value for the 2nd and 4th point have been swapped. 😕

The points have been swapped. I think the vertices are in the reverse order.

Oh right, I stared at it too long and confused myself. 😄

Is the new ordering correct or is that something that needs to be fixed yet?

Is the new ordering correct or is that something that needs to be fixed yet?

It's the same rectangle in the opposite winding order (P4 P3 P2 P1 P4 instead of P1 P2 P3 P4 P1). I can probably match the order, but only up to the starting vertex, e.g. P2 P3 P4 P1 P2.

frewsxcv · 2023-12-04T19:23:35Z

@lnicola Can you add more context to the description of the PR (or in a comment) about why we're avoiding atan2?

urschrei · 2023-12-04T20:17:42Z

(This was all my idea). Atan2 is a comparatively very expensive operation (around 100 cpu cycles on x86, vs around 20 for something like sqrt). Using this approach (inspired by the approach from https://gis.stackexchange.com/a/22934) avoids that. Unfortunately, it's currently suffering from floating point precision issues.

frewsxcv · 2024-01-19T01:34:21Z

Do we see any path forward with the precision issues here? If not, @lnicola what do you think about converting this to a draft PR or closing it?

lnicola · 2024-01-30T19:26:50Z

I finally got to run some benchmarks (on the plain version, not the one using FMA):

minimum rotated rect f64/old/10
                        time:   [1.2758 µs 1.2777 µs 1.2798 µs]
minimum rotated rect f64/new/10
                        time:   [648.70 ns 649.96 ns 651.92 ns]
Found 3 outliers among 100 measurements (3.00%)
  2 (2.00%) low mild
  1 (1.00%) high severe
minimum rotated rect f64/old/100
                        time:   [32.838 µs 32.873 µs 32.908 µs]
Found 5 outliers among 100 measurements (5.00%)
  2 (2.00%) high mild
  3 (3.00%) high severe
minimum rotated rect f64/new/100
                        time:   [18.851 µs 18.854 µs 18.857 µs]
Found 1 outliers among 100 measurements (1.00%)
  1 (1.00%) high severe
Benchmarking minimum rotated rect f64/old/1000: Warming up for 3.0000 s
Warning: Unable to complete 100 samples in 5.0s. You may wish to increase target time to 9.7s, enable flat sampling, or reduce sample count to 50.
minimum rotated rect f64/old/1000
                        time:   [1.9078 ms 1.9094 ms 1.9111 ms]
Found 15 outliers among 100 measurements (15.00%)
  11 (11.00%) low mild
  2 (2.00%) high mild
  2 (2.00%) high severe
Benchmarking minimum rotated rect f64/new/1000: Warming up for 3.0000 s
Warning: Unable to complete 100 samples in 5.0s. You may wish to increase target time to 6.5s, enable flat sampling, or reduce sample count to 60.
minimum rotated rect f64/new/1000
                        time:   [1.2859 ms 1.2866 ms 1.2873 ms]
Found 12 outliers among 100 measurements (12.00%)
  8 (8.00%) low severe
  2 (2.00%) low mild
  2 (2.00%) high severe

So it's about 33% faster, which is noticeable, but not earthshaking.

As for the precision issue, I can't prove it, but I suspect the new method is more precise in some cases and less in others. Looking at the numbers, they're something like:

1.7000000000000006 24.600000000000000
1.7000000000000004 24.600000000000005

14.649854163628408 25.153412571095235
14.649854163628410 25.153412571095235

14.400000000000000 31.000000000000000
14.400000000000002 31.000000000000004 # <- here

1.4501458363715918 30.446587428904767
1.4501458363715916 30.446587428904774

The only reason I suspect the new method is less accurate is that one vertex which is part of the input.

Anyway, I don't know how to improve on this, but I don't insist on merging it either, so we can close this if nobody wants it.

urschrei · 2024-01-31T00:42:40Z

From a practical perspective MRR calculation tends to be something done at scale (e.g. building footprints in a city etc), so a 33 % improvement on something that might take 10 seconds is pretty respectable. I'm gonna try some more comparisons with JTS and GEOS tomorrow and see whether it's a consistent accuracy issue, although it's hard to see what else we could do without re-introducing complexity.

urschrei · 2024-01-31T12:38:03Z

input:

MULTIPOINT ((3.3 30.4), (1.7 24.6), (13.4 25.1), (14.4 31), (3.3 30.4))

geo:

POLYGON ((1.7000000000000004 24.600000000000005, 14.64985416362841 25.153412571095235, 14.400000000000002 31.000000000000004, 1.4501458363715916 30.446587428904774, 1.7000000000000004 24.600000000000005))

JTS (HEAD):

POLYGON ((1.3920611798980265 30.296868171886377, 1.7071376547705703 24.467953386744355, 14.715076474872582 25.171085214857957, 14.40000000000005 30.999999999999982, 1.3920611798980265 30.296868171886377))

GEOS (3.12.0)

POLYGON((14.649854163628415 25.153412571095224,1.7 24.6,1.450145836371595 30.44658742890477,14.4 31,14.649854163628415 25.153412571095224))

geo vs JTS (geo is blue)

geo vs GEOS (geo is blue)

urschrei · 2024-01-31T12:59:50Z

Oh but wait, GEOS had a bug until 3.12.1:
libgeos/geos@2efd5ed

I can't easily access 3.12.1, anyone else got it?

lnicola · 2024-01-31T13:00:54Z

I have 3.12.1, can you send me some code?

urschrei · 2024-01-31T13:01:55Z

I have 3.12.1, can you send me some code?

SELECT ST_AsText(ST_OrientedEnvelope('MULTIPOINT ((3.3 30.4), (1.7 24.6), (13.4 25.1), (14.4 31.0), (3.3 30.4))'));

lnicola · 2024-01-31T13:03:52Z

Same as before:

POLYGON((14.649854163628415 25.153412571095224,1.7 24.6,1.450145836371595 30.44658742890477,14.4 31,14.649854163628415 25.153412571095224))

POSTGIS="3.4.1 ca035b9" [EXTENSION] PGSQL="160" GEOS="3.12.1-CAPI-1.18.1" PROJ="9.3.1 NETWORK_ENABLED=OFF URL_ENDPOINT=https://cdn.proj.org USER_WRITABLE_DIRECTORY=/var/lib/postgresql/.local/share/proj DATABASE_PATH=/usr/local/share/proj/proj.db" LIBXML="2.9.14" LIBJSON="0.16" LIBPROTOBUF="1.4.1" WAGYU="0.5.0 (Internal)" TOPOLOGY

urschrei · 2024-01-31T13:10:50Z

Even though there's no visual difference with 3.12.0 at this resolution:

geo vs GEOS (3.12.1):

lnicola · 2024-01-31T13:17:31Z

It's exactly the same output as with 3.12.0.

urschrei · 2024-01-31T13:19:12Z

Yeah I just saw 🤦‍♂️. Did you run it via PostGIS or using the cpp function?

lnicola · 2024-01-31T13:19:51Z

Using docker.io/imresamu/postgis:16-recent

lnicola changed the title ~~Avoid atan2 in mrr~~ Avoid atan2 in minimum_rotated_rect Oct 20, 2023

lnicola force-pushed the mbr-variant branch from 20225da to 7c6a641 Compare October 20, 2023 08:51

urschrei reviewed Oct 20, 2023

View reviewed changes

lnicola force-pushed the mbr-variant branch 2 times, most recently from bf8ea2a to 6f9708d Compare October 20, 2023 16:30

michaelkirk reviewed Oct 27, 2023

View reviewed changes

frewsxcv marked this pull request as draft January 28, 2024 18:43

lnicola added 2 commits February 3, 2024 18:42

Avoid atan2 in mrr

5429bd7

Use FMA

626dce6

lnicola force-pushed the mbr-variant branch from 6f9708d to 626dce6 Compare February 3, 2024 16:42

Update changelog

3ed0cfa

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Avoid `atan2` in `minimum_rotated_rect` #1094

Avoid `atan2` in `minimum_rotated_rect` #1094

lnicola commented Oct 20, 2023

lnicola commented Oct 20, 2023

urschrei Oct 20, 2023 •

edited

lnicola Oct 20, 2023 •

edited

lnicola Oct 20, 2023

urschrei Oct 20, 2023

lnicola Oct 20, 2023

michaelkirk Oct 27, 2023

lnicola Oct 27, 2023

michaelkirk Oct 27, 2023

lnicola Jan 30, 2024

frewsxcv commented Dec 4, 2023

urschrei commented Dec 4, 2023

frewsxcv commented Jan 19, 2024

lnicola commented Jan 30, 2024

urschrei commented Jan 31, 2024

urschrei commented Jan 31, 2024 •

edited

urschrei commented Jan 31, 2024

lnicola commented Jan 31, 2024

urschrei commented Jan 31, 2024

lnicola commented Jan 31, 2024 •

edited

urschrei commented Jan 31, 2024

lnicola commented Jan 31, 2024

urschrei commented Jan 31, 2024

lnicola commented Jan 31, 2024

Avoid atan2 in minimum_rotated_rect #1094

Are you sure you want to change the base?

Avoid atan2 in minimum_rotated_rect #1094

Conversation

lnicola commented Oct 20, 2023

lnicola commented Oct 20, 2023

urschrei Oct 20, 2023 • edited

Choose a reason for hiding this comment

lnicola Oct 20, 2023 • edited

Choose a reason for hiding this comment

lnicola Oct 20, 2023

Choose a reason for hiding this comment

urschrei Oct 20, 2023

Choose a reason for hiding this comment

lnicola Oct 20, 2023

Choose a reason for hiding this comment

michaelkirk Oct 27, 2023

Choose a reason for hiding this comment

lnicola Oct 27, 2023

Choose a reason for hiding this comment

michaelkirk Oct 27, 2023

Choose a reason for hiding this comment

lnicola Jan 30, 2024

Choose a reason for hiding this comment

frewsxcv commented Dec 4, 2023

urschrei commented Dec 4, 2023

frewsxcv commented Jan 19, 2024

lnicola commented Jan 30, 2024

urschrei commented Jan 31, 2024

urschrei commented Jan 31, 2024 • edited

urschrei commented Jan 31, 2024

lnicola commented Jan 31, 2024

urschrei commented Jan 31, 2024

lnicola commented Jan 31, 2024 • edited

urschrei commented Jan 31, 2024

lnicola commented Jan 31, 2024

urschrei commented Jan 31, 2024

lnicola commented Jan 31, 2024

Avoid `atan2` in `minimum_rotated_rect` #1094

Avoid `atan2` in `minimum_rotated_rect` #1094

urschrei Oct 20, 2023 •

edited

lnicola Oct 20, 2023 •

edited

urschrei commented Jan 31, 2024 •

edited

lnicola commented Jan 31, 2024 •

edited