Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Avoid atan2 in minimum_rotated_rect #1094

Draft
wants to merge 3 commits into
base: main
Choose a base branch
from
Draft

Avoid atan2 in minimum_rotated_rect #1094

wants to merge 3 commits into from

Conversation

lnicola
Copy link
Member

@lnicola lnicola commented Oct 20, 2023

  • 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
Copy link
Member Author

lnicola commented Oct 20, 2023

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.

@lnicola lnicola changed the title Avoid atan2 in mrr Avoid atan2 in minimum_rotated_rect Oct 20, 2023
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;
Copy link
Member

@urschrei urschrei Oct 20, 2023

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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?

Copy link
Member Author

@lnicola lnicola Oct 20, 2023

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Ooh, good idea, forgot about FMA.

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Bah.

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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 }])

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Sorry to resurrect an old PR, but is herbie something we should consider adding to our CI lint?

I'd propose it should be non-failing, but advisory.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

If only! The Herbie lint in Rust hasn't worked since sometime in 2016 when the lint framework got rewritten, and nobody has had the energy to resurrect it. I actually asked about it a couple of years ago on URLO and everyone (including one of the Herbie authors) was enthusiastic, but it was clear that I was gonna be on my own. I generated the suggestion in this PR by manually writing it out and feeding it to the online tool.

@lnicola lnicola force-pushed the mbr-variant branch 2 times, most recently from bf8ea2a to 6f9708d Compare October 20, 2023 16:30
(14.649854163628408, 25.153412571095235),
(1.7000000000000006, 24.6),
(1.7000000000000004, 24.600000000000005),
(14.64985416362841, 25.153412571095235),
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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?

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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
Copy link
Member

frewsxcv commented Dec 4, 2023

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

@urschrei
Copy link
Member

urschrei commented Dec 4, 2023

(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
Copy link
Member

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?

@frewsxcv frewsxcv marked this pull request as draft January 28, 2024 18:43
@lnicola
Copy link
Member Author

lnicola commented Jan 30, 2024

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

image

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
Copy link
Member

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
Copy link
Member

urschrei commented Jan 31, 2024

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)
Screenshot 2024-01-31 at 12 35 23

geo vs GEOS (geo is blue)
Screenshot 2024-01-31 at 12 37 07

@urschrei
Copy link
Member

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
Copy link
Member Author

lnicola commented Jan 31, 2024

I have 3.12.1, can you send me some code?

@urschrei
Copy link
Member

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
Copy link
Member Author

lnicola commented Jan 31, 2024

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
Copy link
Member

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

geo vs GEOS (3.12.1):

Screenshot 2024-01-31 at 13 08 53

@lnicola
Copy link
Member Author

lnicola commented Jan 31, 2024

It's exactly the same output as with 3.12.0.

@urschrei
Copy link
Member

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

@lnicola
Copy link
Member Author

lnicola commented Jan 31, 2024

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

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

4 participants