fitodic / centerline

Calculate the polygon's centerline
https://centerline.readthedocs.io
MIT License
266 stars 55 forks source link

Ran Into this issue while running - Not enough points to construct initial simplex #9

Closed 3vivekb closed 5 years ago

3vivekb commented 7 years ago

(/Users/vivek/anaconda3/envs/centerline27) bash-3.2$ shp2centerline 84/sidewalks_w_id/sidewalks_w_id.shp 84/output/sidewalk_output.shp
Importing polygons from: 84/sidewalks_w_id/sidewalks_w_id.shp
Calculating centerlines.
Traceback (most recent call last):
  File "/Users/vivek/anaconda3/envs/centerline27/bin/shp2centerline", line 94, in <module>
    Shp2centerline(args.SRC, args.DEST, args.BORDENS)
  File "/Users/vivek/anaconda3/envs/centerline27/bin/shp2centerline", line 20, in __init__
    self.run()
  File "/Users/vivek/anaconda3/envs/centerline27/bin/shp2centerline", line 38, in run
    self.dct_centerlines[key] = centerlineObj.createCenterline()
  File "/Users/vivek/anaconda3/envs/centerline27/lib/python2.7/site-packages/centerline/centerline.py", line 32, in createCenterline
    vor = Voronoi(border)
  File "scipy/spatial/qhull.pyx", line 2544, in scipy.spatial.qhull.Voronoi.__init__ (scipy/spatial/qhull.c:25511)
  File "scipy/spatial/qhull.pyx", line 434, in scipy.spatial.qhull._Qhull.__init__ (scipy/spatial/qhull.c:5799)
scipy.spatial.qhull.QhullError: QH6214 qhull input error: not enough points(2) to construct initial simplex (need 4)

While executing:  | qhull v Qbb Qz Qc
Options selected for Qhull 2015.2.r 2016/01/18:
  run-id 1309006505  voronoi  Qbbound-last  Qz-infinity-point  Qcoplanar-keep
  _pre-merge  _zero-centrum  Qinterior-keep```
fitodic commented 7 years ago

Have you identified the geometric feature that raises this error? How many points does its border consist of?

On the other hand, have you tried modifying the border density parameter?

nachomat commented 6 years ago

Hi fitodic, mind you emailing me? I have a few questions to do but I'm affraid here is not the best forum to do.

fitodic commented 6 years ago

@nachomat Of course not, but other users might be interested in what you have to say. As far as I can see, there are at least two other users facing similar issues.

Plus, there is no email address on your profile :)

HenryW95 commented 6 years ago

I have been getting a similar Qhull error ("QH6228 - Qhull internal error") for one polygon in a shapefile. The shp2centerline command works perfectly for thousands of other polygons in the shapefile, but throws an error on this one. Attached is the command line error message. Additionally, the error only occurs with a specified border density of 6. It works with 5, 7, 6.1, etc.

Here is a link to the polygon (extracted from my larger shapefile) that is causing the error: https://www.dropbox.com/s/rd340o5t43u2wcj/error_feature.zip?dl=0

qhullerror
fitodic commented 6 years ago

@HenryW95 Apologies for the late response and thank you for the detailed explanation as well as providing the troublesome geometry. 👍

This issue might be connected to the griddata issue, but I'll know more once I debug it. I'll take a look as soon as I find the time.

I'm currently in the process of reorganizing the codebase and improving test coverage. If I were able to reproduce the bug, would you mind if I add the data you've provided to the repository to be used as input test data?

HenryW95 commented 6 years ago

@fitodic Thanks for the info! Definitely feel free to use the data.

fitodic commented 6 years ago

@HenryW95 I've tried reproducing this issue, but unfortunately without much success. I am using the latest centerline package (v0.3) and have run the following command:

$ create_centerlines tests/data/issue_9/error_feature.shp tests/data/issue_9/centerline.shp 6

which successfully produced the centerline Shapefile: issue_9.zip.

What versions of scipy and numpy are you using? Have you tried updating the centerline package to v0.3+ and rerunning the command? I would also recommend updating the rest of the requirements, especially scipy and numpy.

MattFerraro commented 5 years ago

Just weighing in months later that I'm also seeing this problem when I'm inputing an extremely simple shape that I traced by hand on geojson.io.

Shape is here: https://gist.github.com/MattFerraro/227c13924f60288caaf970cd8c95b143

I'm on a mac with gdal 2.1 and running in a virtualenv with centerline==0.5

fitodic commented 5 years ago

@MattFerraro What border density were you using?

MattFerraro commented 5 years ago

just the default

fitodic commented 5 years ago

@MattFerraro Have you tried decreasing it below the default (0.5)?

MattFerraro commented 5 years ago

Sure, I get the same error at 0.4, 0.3, 0.2, 0.1, 0.05, 0.04, 0.03, 0.02, then finally at 0.01 I get a different error:

WARNING:root:ignoring record that could not be processed: QH6154 Qhull precision error: Initial simplex is flat (facet 2 is coplanar with the interior point)

While executing:  | qhull v Qbb Qz Qc
Options selected for Qhull 2015.2.r 2016/01/18:
  run-id 1335308354  voronoi  Qbbound-last  Qz-infinity-point  Qcoplanar-keep
  _pre-merge  _zero-centrum  Qinterior-keep  Pgood  _max-width 0.0006
  Error-roundoff 6.2e-16  _one-merge 4.4e-15  Visible-distance 1.2e-15
  U-coplanar-distance 1.2e-15  Width-outside 2.5e-15  _wide-facet 7.5e-15

precision problems (corrected unless 'Q0' or an error)
      2 flipped facets
      2 degenerate hyperplanes recomputed with gaussian elimination
      4 nearly singular or axis-parallel hyperplanes
      2 zero divisors during back substitute
      2 zero divisors during gaussian elimination

The input to qhull appears to be less than 3 dimensional, or a
computation has overflowed.

Qhull could not construct a clearly convex simplex from points:
- p2(v3): -0.45  0.39     0
- p3(v2): -0.45  0.39 0.0006
- p0(v1): -0.45  0.39     0
- p1(v0): -0.45  0.39 6.4e-06

The center point is coplanar with a facet, or a vertex is coplanar
with a neighboring facet.  The maximum round off error for
computing distances is 6.2e-16.  The center point, facets and distances
to the center point are as follows:

center point  -0.4485   0.3879 0.0001512

facet p3 p0 p1 distance= -5.4e-17
facet p2 p0 p1 distance=    0
facet p2 p3 p1 distance= -1.1e-16
facet p2 p3 p0 distance=    0

These points either have a maximum or minimum x-coordinate, or
they maximize the determinant for k coordinates.  Trial points
are first selected from points that maximize a coordinate.

The min and max coordinates for each dimension are:
  0:   -0.4489  -2.225e-308  difference= 0.4489
  1:    0.3878     0.388  difference= 0.0002036
  2:         0  0.0005984  difference= 0.0005984

If the input should be full dimensional, you have several options that
may determine an initial simplex:
  - use 'QJ'  to joggle the input and make it full dimensional
  - use 'QbB' to scale the points to the unit cube
  - use 'QR0' to randomly rotate the input for different maximum points
  - use 'Qs'  to search all points for the initial simplex
  - use 'En'  to specify a maximum roundoff error less than 6.2e-16.
  - trace execution with 'T3' to see the determinant for each point.

If the input is lower dimensional:
  - use 'QJ' to joggle the input and make it full dimensional
  - use 'Qbk:0Bk:0' to delete coordinate k from the input.  You should
    pick the coordinate with the least range.  The hull will have the
    correct topology.
  - determine the flat containing the points, rotate the points
    into a coordinate plane, and delete the other coordinates.
  - add one or more points to make the input full dimensional.
fitodic commented 5 years ago

@MattFerraro I've successfully created the centerline based on the geometry you've provided: centerline.zip

The result was achieved using the following command:

$ create_centerlines input.geojson centerline.geojson 0.00001

This seems to be an issue of inconsistency in units. The input data is in EPSG:4326 which uses degrees, whereas the create_centerline script's default border density is in meters, as stated in the documentation.

That leaves us with two options:

  1. Specify the border density in degrees
  2. Transform the input data into a metric system (e.g. EPSG:3857)

I'll update the documentation to reflect this issue.