Error in convex_hull_image: QH6228 Qhull internal error (qh_findbestlower)

[jmetz]

I came across the following bug when calling convex_hull_objects, which I found was caused by a call to convex_hull_image on one of the labelled objects.

Below is a minimal(ish) snippet to reproduce the error.

I get

Success on cropped
Success on Expanded and shifted
QH6228 
Qhull internal error (qh_findbestlower): all neighbors of facet 3 are flipped or upper Delaunay.

import numpy as np
import skimage.morphology as skmo

nonzeros = (([1367, 1368, 1368, 1368, 1369, 1369, 1369, 1369, 1369, 1370, 1370,
       1370, 1370, 1370, 1370, 1370, 1371, 1371, 1371, 1371, 1371, 1371,
       1371, 1371, 1371, 1372, 1372, 1372, 1372, 1372, 1372, 1372, 1372,
       1372, 1373, 1373, 1373, 1373, 1373, 1373, 1373, 1373, 1373, 1374,
       1374, 1374, 1374, 1374, 1374, 1374, 1375, 1375, 1375, 1375, 1375,
       1376, 1376, 1376, 1377]),
    ([151, 150, 151, 152, 149, 150, 151, 152, 153, 148, 149, 150, 151,
       152, 153, 154, 147, 148, 149, 150, 151, 152, 153, 154, 155, 146,
       147, 148, 149, 150, 151, 152, 153, 154, 146, 147, 148, 149, 150,
       151, 152, 153, 154, 147, 148, 149, 150, 151, 152, 153, 148, 149,
       150, 151, 152, 149, 150, 151, 150]))

image = np.zeros((1392, 1040))
image[ nonzeros ] = 1

imageCropped = image[1365:1380, 144:160].copy()
imageCropped = imageCropped.astype(bool)
chull_im = skmo.convex_hull_image(imageCropped)
print "Success on cropped"

imageBigShift = np.zeros((5000,5000), dtype=image.dtype)
imageBigShift[1365:1380, 1144:1160] = imageCropped.copy()
chull_im = skmo.convex_hull_image(imageBigShift)
print "Success on Expanded and shifted"

imageBig = np.zeros((5000,5000), dtype=image.dtype)
imageBig[1365:1380, 144:160] = imageCropped.copy()

chull_im = skmo.convex_hull_image(imageBig)
print "Success on Expanded"

chull_im = skmo.convex_hull_image(image)
print "Success on original"

[stefanv]

This works for me with SciPy 0.12.0. Which version are you using?

[jmetz]

Also 0.12.0. The issue is being caused by the call to qhull, perhaps my version has a bug that has since been fixed... I have qhull version (2009.1-3) (from ubuntu 13.04 repos).

[jmetz]

I just upgraded to scipy 0.13.2, and also found that scipy 0.12.0 was also using qhull 2012.1 2012/02/18 (as is 0.13.2) - so it's not an old qhull bug... But I get the same error regardless.

The full qhull error message is:

Success on cropped
Success on Expanded and shifted
QH6228 
Qhull internal error (qh_findbestlower): all neighbors of facet 3 are flipped or upper Delaunay.
Please report this error to qhull_bug@qhull.org with the input and all of the output.
ERRONEOUS FACET:
- f3
    - flags: top simplicial upperDelaunay
    - normal:   -0.6912   -0.654   0.3075
    - offset:    1043.28
    - vertices: p156(v3) p5(v2) p78(v0)
    - neighboring facets: f26 f17 f7

While executing:  | qhull d Qbb Qz Qc Qt
Options selected for Qhull 2012.1 2012/02/18:
  run-id 535251457  delaunay  Qbbound-last  Qz-infinity-point  Qcoplanar-keep
  Qtriangulate  _pre-merge  _zero-centrum  Qinterior-keep  Pgood  _max-width 11
  Error-roundoff 1.3e-12  _one-merge 9.4e-12  Visible-distance 2.7e-12
  U-coplanar-distance 2.7e-12  Width-outside 5.4e-12  _wide-facet 1.6e-11
Last point added to hull was p97.

At error exit:

Delaunay triangulation by the convex hull of 157 points in 3-d:

  Number of input sites and at-infinity: 12
  Number of nearly incident points: 12
  Number of Delaunay regions: 0

Statistics for:  | qhull d Qbb Qz Qc Qt

  Number of points processed: 12
  Number of hyperplanes created: 39
  Number of facets in hull: 20
  Number of distance tests for qhull: 1308

precision problems (corrected unless 'Q0' or an error)
     11 coplanar points during partitioning
Traceback (most recent call last):
  File "convex_hull_error.py", line 37, in <module>
    chull_im = skmo.convex_hull_image(imageBig)
  File "/usr/local/lib/python2.7/dist-packages/scikit_image-0.9dev-py2.7-linux-x86_64.egg/skimage/morphology/convex_hull.py", line 55, in convex_hull_image
    chull = Delaunay(coords).convex_hull
  File "qhull.pyx", line 1701, in scipy.spatial.qhull.Delaunay.__init__ (scipy/spatial/qhull.c:13548)
  File "qhull.pyx", line 323, in scipy.spatial.qhull._Qhull.__init__ (scipy/spatial/qhull.c:3592)
scipy.spatial.qhull.QhullError: Qhull error

[jmetz]

I've found a workaround; as the full error implies this is a precision issue, so I tried subtracting the mean as follows and it works (NB as well as subtracting the mean a full fix should probably scale the set of points into [-1,1]).

My temporary fix is essentially:

# skimage/morphology/convex_hull.py line 53

    # Find the convex hull
    offset = coords.mean(axis=0)
    coords -= offset
    chull = Delaunay(coords).convex_hull
    v = coords[np.unique(chull)]

    # Sort vertices clock-wise
    v_centred = v - v.mean(axis=0)
    angles = np.arctan2(v_centred[:, 0], v_centred[:, 1])
    v = v[np.argsort(angles)]

    # Add back here 
    v += offset

As I hinted, as well as subtracting the offset, a full fix would possibly also scale ( ~ /= coords.max(axis=0) ) though this could introduce other issues.

I'm going to wait and see if anyone else confirms this bug before submitting the patch...

[stefanv]

Is this a bug in QHull, or is that expected behavior?

[jmetz]

I don't know enough to be sure; but I also submitted a query to qhull_bug@qhull.org . I'd consider it a bug given the results from the test code, at least for the purposes of scikit-image.

[jmetz]

This has also been reported here: https://gist.github.com/automata/5447635, so I'm going to put in a pull request with my workaround.