mbc breaks on windows and solaris

[edzer]

See https://cran.r-project.org/web/checks/check_results_lwgeom.html ; this was already outcommented because it broke on windows (appveyor).

@rundel since this was your PR, could you please take a look what is going on here?

[rundel]

I'll take a look, do you recall what the windows error looked like (similar to solaris?)

The build failure on macos box looks like a failure to find libgeos_c, any idea about what is causing that? configure.ac looks fine as far as I can tell.

[edzer]

Yes, identical. I don't have access to windows, did the development/testing through appveyor.

[rundel]

Weird, I just rebuilt from github on a windows vm and I'm not seeing any issue. Will try with the cran version now and see what happens.

[rundel]

So on 32-bit R on Windows I get the following (but not on 64 bit),

> lwgeom:::CPL_minimum_bounding_circle(state)
$center
$center[[1]]
  x   y 
NaN NaN 


$radius
[1] NaN

seems like it may be a lwgeom bug?

[edzer]

Weird indeed!

[edzer]

But why on 32 bits? Solaris is also 32 bits.

[rundel]

Seems like the short term fix is to put a check into the R code to provide an error on NaN results from the C++ code. I'll put up a pull request shortly.

I also opened r-spatial/sf#561 since it seems like sf should also handle the empty geometrys with either an empty plot or at least a more explicit error message.

[rundel]

Definitely a 32 bit issue, I was able to get a 32 bit ubuntu docker container up and running and I'm seeing the exact same behavior. Time to poke at it with lldb and see what I can see.

[rundel]

So gdb/lldb didn't like docker for some reason, but many many printf statements later it turns out the issue was a trivial floating point comparison issue.

When checking if a new vertex was within the bounding circle it occasionally happened that the new vertex was identical to one of the points used as the circle's support and hence fell exactly on the edge of the circle. Due to the 32 bit double precision it was possible then that (distance2d_pt_pt(p, c->center) > c->radius) evaluated as true even though (distance2d_pt_pt(p, c->center) - c->radius) == 0 . This erroneous vertex would then be used as a support point which mucked up the circumcenter calculation due to a divide by zero which results in NaNs everywhere.

Changing distance2d_pt_pt(p, c->center) > c->radius to distance2d_pt_pt(p, c->center) - c->radius > DBL_EPSILON along with including float.h fixes the issue as far as I can tell. Change is included in #9.

Clearly this is an upstream postgis issue, and I don't see any mention of it in their trac system. Unfortunately I don't have an osgeo account and their request system is borked at the moment so I haven't been able to post the bug report.

R based repex:

library(lwgeom)
library(sf)

nc = st_read(system.file("shape/nc.shp", package="sf"))
state = st_union(st_geometry(nc))

st_minimum_bounding_circle(state)

Sample WKT that produces the issue on a 32 bit systems can be found here.

[edzer]

Good catch! I'd suggest doing a PR on the github repo of postgis/liblwgeom, referring here.

[rundel]

Opened a ticket on postgis' trac, https://trac.osgeo.org/postgis/ticket/3930

[rundel]

Officially added to the postgis 2.4.3 milestone - https://trac.osgeo.org/postgis/query?id=3827%2C3700%2C3906%2C3787%2C3899%2C3914%2C3930%2C2960