Skip to content
Browse files

BUG: spatial: ensure that qhull._find_simplex_directed always terminates

If the triangulation contains very narrow simplices, this may introduce
rounding error in the barycentric transformations, making the algorithm
think that two triangles are mutually closer to the target point. This
causes an infinite loop.

The solution is to limit the iteration count, and fall back to brute
force if the algorithm fails. More sophisticated cycle detection might
also be made --- however, the brute force O(N) cost probably overwhelms
any gains from such approaches, and adding cycle detection would
slow down the fast path, so it's probably best to keep it simple.
  • Loading branch information...
1 parent 1e7bea8 commit bb9ea9e37a4e569e9e2ec6d43ebd9b5787e536be @pv pv committed
Showing with 17 additions and 2 deletions.
  1. +17 −2 scipy/spatial/qhull.pyx
View
19 scipy/spatial/qhull.pyx
@@ -686,8 +686,12 @@ cdef int _find_simplex_directed(DelaunayInfo_t *d, double *c,
6) If all barycentric coordinates are non-negative but 5) is not true,
we are in an inconsistent situation -- this should never happen.
+ This may however enter an infinite loop due to rounding errors in
+ the computation of the barycentric coordinates, so the iteration
+ count needs to be limited, and a fallback to brute force provided.
+
"""
- cdef int k, m, ndim, inside, isimplex
+ cdef int k, m, ndim, inside, isimplex, cycle_k
cdef double *transform
cdef double v
@@ -697,7 +701,15 @@ cdef int _find_simplex_directed(DelaunayInfo_t *d, double *c,
if isimplex < 0 or isimplex >= d.nsimplex:
isimplex = 0
- while isimplex != -1:
+ # The maximum iteration count: it should be large enough so that
+ # the algorithm usually succeeds, but smaller than nsimplex so
+ # that for the cases where the algorithm fails, the main cost
+ # still comes from the brute force search.
+
+ for cycle_k in range(1 + d.nsimplex//4):
+ if isimplex == -1:
+ break
+
transform = d.transform + isimplex*ndim*(ndim+1)
inside = 1
@@ -733,6 +745,9 @@ cdef int _find_simplex_directed(DelaunayInfo_t *d, double *c,
# fall back to brute force
isimplex = _find_simplex_bruteforce(d, c, x, eps)
break
+ else:
+ # the algorithm failed to converge -- fall back to brute force
+ isimplex = _find_simplex_bruteforce(d, c, x, eps)
start[0] = isimplex
return isimplex

0 comments on commit bb9ea9e

Please sign in to comment.
Something went wrong with that request. Please try again.