[Scipy-svn] r6653 - in trunk/scipy/spatial: . tests
scipy-svn at scipy.org
scipy-svn at scipy.org
Tue Aug 31 17:55:50 EDT 2010
Author: ptvirtan
Date: 2010-08-31 16:55:50 -0500 (Tue, 31 Aug 2010)
New Revision: 6653
Added:
trunk/scipy/spatial/generate_qhull.py
trunk/scipy/spatial/qhull.pxd
trunk/scipy/spatial/qhull.pyx
trunk/scipy/spatial/qhull_blas.h
trunk/scipy/spatial/tests/test_qhull.py
Modified:
trunk/scipy/spatial/SConscript
trunk/scipy/spatial/__init__.py
trunk/scipy/spatial/setup.py
Log:
spatial: add wrappers for qhull's Delaunay tesselation
Add a wrapper for the Delaunay tesselation routine from Qhull.
Modified: trunk/scipy/spatial/SConscript
===================================================================
--- trunk/scipy/spatial/SConscript 2010-08-31 21:54:55 UTC (rev 6652)
+++ trunk/scipy/spatial/SConscript 2010-08-31 21:55:50 UTC (rev 6653)
@@ -10,3 +10,12 @@
env.NumpyPythonExtension('_distance_wrap',
source = [join('src', 'distance_wrap.c'),
join('src', 'distance.c')])
+
+# Build qhull
+src = [pjoin('qhull', 'src', s) for s in [
+ 'geom.c', 'geom2.c', 'global.c', 'io.c', 'mem.c',
+ 'merge.c', 'poly.c', 'poly2.c', 'qset.c', 'user.c',
+ 'stat.c', 'qhull.c']]
+qhullsrc = env.DistutilsStaticExtLibrary('qhullsrc', source=src)
+
+env.NumpyPythonExtension('qhull', source = 'qhull.c', LIBS="qhullsrc")
Modified: trunk/scipy/spatial/__init__.py
===================================================================
--- trunk/scipy/spatial/__init__.py 2010-08-31 21:54:55 UTC (rev 6652)
+++ trunk/scipy/spatial/__init__.py 2010-08-31 21:55:50 UTC (rev 6653)
@@ -5,6 +5,7 @@
from info import __doc__
from kdtree import *
from ckdtree import *
+from qhull import *
__all__ = filter(lambda s:not s.startswith('_'),dir())
__all__ += ['distance']
Added: trunk/scipy/spatial/generate_qhull.py
===================================================================
--- trunk/scipy/spatial/generate_qhull.py (rev 0)
+++ trunk/scipy/spatial/generate_qhull.py 2010-08-31 21:55:50 UTC (rev 6653)
@@ -0,0 +1,29 @@
+#!/usr/bin/env python
+import tempfile
+import subprocess
+import os
+import sys
+import re
+import shutil
+
+tmp_dir = tempfile.mkdtemp()
+try:
+ # Run Cython
+ dst_fn = os.path.join(tmp_dir, 'qhull.c')
+ ret = subprocess.call(['cython', '-o', dst_fn, 'qhull.pyx'])
+ if ret != 0:
+ sys.exit(ret)
+
+ # Strip comments
+ f = open(dst_fn, 'r')
+ text = f.read()
+ f.close()
+
+ r = re.compile(r'/\*(.*?)\*/', re.S)
+
+ text = r.sub('', text)
+ f = open('qhull.c', 'w')
+ f.write(text)
+ f.close()
+finally:
+ shutil.rmtree(tmp_dir)
Property changes on: trunk/scipy/spatial/generate_qhull.py
___________________________________________________________________
Name: svn:executable
+ *
Added: trunk/scipy/spatial/qhull.pxd
===================================================================
--- trunk/scipy/spatial/qhull.pxd (rev 0)
+++ trunk/scipy/spatial/qhull.pxd 2010-08-31 21:55:50 UTC (rev 6653)
@@ -0,0 +1,89 @@
+# -*-cython-*-
+"""
+Qhull shared definitions, for use by other Cython modules
+
+"""
+#
+# Copyright (C) Pauli Virtanen, 2010.
+#
+# Distributed under the same BSD license as Scipy.
+#
+
+cdef extern from "stdlib.h":
+ void *malloc(int size)
+ void free(void *ptr)
+
+cdef extern from "numpy/ndarraytypes.h":
+ cdef enum:
+ NPY_MAXDIMS
+
+ctypedef struct DelaunayInfo_t:
+ int ndim
+ int npoints
+ int nsimplex
+ double *points
+ int *vertices
+ int *neighbors
+ double *equations
+ double *transform
+ int *vertex_to_simplex
+ double paraboloid_scale
+ double paraboloid_shift
+ double *max_bound
+ double *min_bound
+
+cdef DelaunayInfo_t *_get_delaunay_info(obj, int compute_transform,
+ int compute_vertex_to_simplex)
+
+#
+# N-D geometry
+#
+
+cdef int _barycentric_inside(int ndim, double *transform,
+ double *x, double *c, double eps) nogil
+
+cdef void _barycentric_coordinate_single(int ndim, double *transform,
+ double *x, double *c, int i) nogil
+
+cdef void _barycentric_coordinates(int ndim, double *transform,
+ double *x, double *c) nogil
+
+#
+# N+1-D geometry
+#
+
+cdef void _lift_point(DelaunayInfo_t *d, double *x, double *z) nogil
+
+cdef double _distplane(DelaunayInfo_t *d, int isimplex, double *point) nogil
+
+#
+# Finding simplices
+#
+
+cdef int _is_point_fully_outside(DelaunayInfo_t *d, double *x, double eps) nogil
+
+cdef int _find_simplex_bruteforce(DelaunayInfo_t *d, double *c, double *x,
+ double eps) nogil
+
+cdef int _find_simplex_directed(DelaunayInfo_t *d, double *c, double *x,
+ int *start, double eps) nogil
+
+cdef int _find_simplex(DelaunayInfo_t *d, double *c, double *x, int *start,
+ double eps) nogil
+
+#
+# Walking ridges connected to a vertex
+#
+
+ctypedef struct RidgeIter2D_t:
+ DelaunayInfo_t *info
+ int vertex
+ int edge
+ int vertex2
+ int triangle
+ int start_triangle
+ int start_edge
+
+cdef void _RidgeIter2D_init(RidgeIter2D_t *it, DelaunayInfo_t *d,
+ int vertex) nogil
+cdef void _RidgeIter2D_next(RidgeIter2D_t *it) nogil
Added: trunk/scipy/spatial/qhull.pyx
===================================================================
--- trunk/scipy/spatial/qhull.pyx (rev 0)
+++ trunk/scipy/spatial/qhull.pyx 2010-08-31 21:55:50 UTC (rev 6653)
@@ -0,0 +1,1160 @@
+"""
+Wrappers for Qhull triangulation, plus some additional N-D geometry utilities
+
+.. versionadded:: 0.9
+
+"""
+#
+# Copyright (C) Pauli Virtanen, 2010.
+#
+# Distributed under the same BSD license as Scipy.
+#
+
+import threading
+import numpy as np
+cimport numpy as np
+cimport cython
+cimport qhull
+
+__all__ = ['Delaunay', 'tsearch']
+
+#------------------------------------------------------------------------------
+# Qhull interface
+#------------------------------------------------------------------------------
+
+cdef extern from "stdio.h":
+ extern void *stdin
+ extern void *stderr
+ extern void *stdout
+
+cdef extern from "qhull/src/qset.h":
+ ctypedef union setelemT:
+ void *p
+ int i
+
+ ctypedef struct setT:
+ int maxsize
+ setelemT e[1]
+
+cdef extern from "qhull/src/qhull.h":
+ ctypedef double realT
+ ctypedef double coordT
+ ctypedef double pointT
+ ctypedef int boolT
+ ctypedef unsigned int flagT
+
+ ctypedef struct facetT:
+ coordT offset
+ coordT *center
+ coordT *normal
+ facetT *next
+ facetT *previous
+ unsigned id
+ setT *vertices
+ setT *neighbors
+ flagT simplicial
+ flagT flipped
+ flagT upperdelaunay
+
+ ctypedef struct vertexT:
+ vertexT *next
+ vertexT *previous
+ unsigned int id, visitid
+ pointT *point
+ setT *neighbours
+
+ ctypedef struct qhT:
+ boolT DELAUNAY
+ boolT SCALElast
+ boolT KEEPcoplanar
+ boolT MERGEexact
+ boolT NOerrexit
+ boolT PROJECTdelaunay
+ boolT ATinfinity
+ int normal_size
+ char *qhull_command
+ facetT *facet_list
+ facetT *facet_tail
+ int num_facets
+ unsigned int facet_id
+ pointT *first_point
+ pointT *input_points
+ realT last_low
+ realT last_high
+ realT last_newhigh
+ realT max_outside
+ realT MINoutside
+ realT DISTround
+
+
+ extern qhT qh_qh
+ extern int qh_PRINToff
+ extern int qh_ALL
+
+ void qh_init_A(void *inp, void *out, void *err, int argc, char **argv)
+ void qh_init_B(realT *points, int numpoints, int dim, boolT ismalloc)
+ void qh_checkflags(char *, char *)
+ void qh_initflags(char *)
+ void qh_option(char *, char*, char* )
+ void qh_freeqhull(boolT)
+ void qh_memfreeshort(int *curlong, int *totlong)
+ void qh_qhull()
+ void qh_check_output()
+ void qh_produce_output()
+ void qh_triangulate()
+ void qh_checkpolygon()
+ void qh_findgood_all()
+ void qh_appendprint(int format)
+ realT *qh_readpoints(int* num, int *dim, boolT* ismalloc)
+ int qh_new_qhull(int dim, int numpoints, realT *points,
+ boolT ismalloc, char* qhull_cmd, void *outfile,
+ void *errfile)
+ int qh_pointid(pointT *point)
+
+# Qhull is not threadsafe: needs locking
+_qhull_lock = threading.Lock()
+
+
+#------------------------------------------------------------------------------
+# LAPACK interface
+#------------------------------------------------------------------------------
+
+cdef extern from "qhull_blas.h":
+ void qh_dgesv(int *n, int *nrhs, double *a, int *lda, int *ipiv,
+ double *b, int *ldb, int *info)
+
+
+#------------------------------------------------------------------------------
+# Delaunay triangulation using Qhull
+#------------------------------------------------------------------------------
+
+def _construct_delaunay(np.ndarray[np.double_t, ndim=2] points):
+ """
+ Perform Delaunay triangulation of the given set of points.
+
+ """
+
+ # Run qhull with the options
+ #
+ # - d: perform delaunay triangulation
+ # - Qbb: scale last coordinate for Delaunay
+ # - Qz: reduces Delaunay precision errors for cospherical sites
+ # - Qt: output only simplical facets (can produce degenerate 0-area ones)
+ #
+ cdef char *options = "qhull d Qz Qbb Qt"
+ cdef int curlong, totlong
+ cdef int dim
+ cdef int numpoints
+ cdef int exitcode
+
+ points = np.ascontiguousarray(points)
+ numpoints = points.shape[0]
+ dim = points.shape[1]
+
+ if numpoints <= 0:
+ raise ValueError("No points to triangulate")
+
+ if dim < 2:
+ raise ValueError("Need at least 2-D data to triangulate")
+
+ _qhull_lock.acquire()
+ try:
+ qh_qh.NOerrexit = 1
+ exitcode = qh_new_qhull(dim, numpoints, <realT*>points.data, 0,
+ options, NULL, stderr)
+ try:
+ if exitcode != 0:
+ raise RuntimeError("Qhull error")
+
+ qh_triangulate() # get rid of non-simplical facets
+
+ if qh_qh.SCALElast:
+ paraboloid_scale = qh_qh.last_newhigh / (
+ qh_qh.last_high - qh_qh.last_low)
+ paraboloid_shift = - qh_qh.last_low * paraboloid_scale
+ else:
+ paraboloid_scale = 1.0
+ paraboloid_shift = 0.0
+
+ vertices, neighbors, equations = \
+ _qhull_get_facet_array(dim, numpoints)
+
+ return (vertices, neighbors, equations,
+ paraboloid_scale, paraboloid_shift)
+ finally:
+ qh_freeqhull(0)
+ qh_memfreeshort(&curlong, &totlong)
+ if curlong != 0 or totlong != 0:
+ raise RuntimeError("qhull: did not free %d bytes (%d pieces)" %
+ (totlong, curlong))
+ finally:
+ _qhull_lock.release()
+
+
+def _qhull_get_facet_array(int ndim, int numpoints):
+ """
+ Return array of simplical facets currently in Qhull.
+
+ Returns
+ -------
+ vertices : array of int, shape (nfacets, ndim+1)
+ Indices of coordinates of vertices forming the simplical facets
+ neighbors : array of int, shape (nfacets, ndim)
+ Indices of neighboring facets. The kth neighbor is opposite
+ the kth vertex, and the first neighbor is the horizon facet
+ for the first vertex.
+
+ Facets extending to infinity are denoted with index -1.
+
+ """
+
+ cdef facetT* facet
+ cdef facetT* neighbor
+ cdef vertexT *vertex
+ cdef int i, j, point
+ cdef np.ndarray[np.int_t, ndim=2] vertices
+ cdef np.ndarray[np.int_t, ndim=2] neighbors
+ cdef np.ndarray[np.double_t, ndim=2] equations
+ cdef np.ndarray[np.int_t, ndim=1] id_map
+
+ id_map = np.empty((qh_qh.facet_id,), dtype=np.int)
+ id_map.fill(-1)
+
+ # Compute facet indices
+ facet = qh_qh.facet_list
+ j = 0
+ while facet and facet.next:
+ if facet.simplicial and not facet.upperdelaunay:
+ id_map[facet.id] = j
+ j += 1
+ facet = facet.next
+
+ # Allocate output
+ vertices = np.zeros((j, ndim+1), dtype=np.int)
+ neighbors = np.zeros((j, ndim+1), dtype=np.int)
+ equations = np.zeros((j, ndim+2), dtype=np.double)
+
+ # Retrieve facet information
+ facet = qh_qh.facet_list
+ j = 0
+ while facet and facet.next:
+ if not facet.simplicial:
+ raise ValueError("non-simplical facet encountered")
+
+ if facet.upperdelaunay:
+ facet = facet.next
+ continue
+
+ # Save vertex info
+ for i in xrange(ndim+1):
+ vertex = <vertexT*>facet.vertices.e[i].p
+ point = qh_pointid(vertex.point)
+ vertices[j, i] = point
+
+ # Save neighbor info
+ for i in xrange(ndim+1):
+ neighbor = <facetT*>facet.neighbors.e[i].p
+ neighbors[j,i] = id_map[neighbor.id]
+
+ # Save simplex equation info
+ for i in xrange(ndim+1):
+ equations[j,i] = facet.normal[i]
+ equations[j,ndim+1] = facet.offset
+
+ j += 1
+ facet = facet.next
+
+ return vertices, neighbors, equations
+
+
+#------------------------------------------------------------------------------
+# Barycentric coordinates
+#------------------------------------------------------------------------------
+
+ at cython.boundscheck(False)
+def _get_barycentric_transforms(np.ndarray[np.double_t, ndim=2] points,
+ np.ndarray[np.int_t, ndim=2] vertices):
+ """
+ Compute barycentric affine coordinate transformations for given
+ simplices.
+
+ Returns
+ -------
+ Tinvs : array, shape (nsimplex, ndim+1, ndim)
+ Barycentric transforms for each simplex.
+
+ Tinvs[i,:ndim,:ndim] contains inverse of the matrix ``T``,
+ and Tinvs[i,ndim,:] contains the vector ``r_n`` (see below).
+
+ Notes
+ -----
+ Barycentric transform from ``x`` to ``c`` is defined by::
+
+ T c = x - r_n
+
+ where the ``r_1, ..., r_n`` are the vertices of the simplex.
+ The matrix ``T`` is defined by the condition::
+
+ T e_j = r_j - r_n
+
+ where ``e_j`` is the unit axis vector, e.g, ``e_2 = [0,1,0,0,...]``
+ This implies that ``T_ij = (r_j - r_n)_i``.
+
+ For the barycentric transforms, we need to compute the inverse
+ matrix ``T^-1`` and store the vectors ``r_n`` for each vertex.
+ These are stacked into the `Tinvs` returned.
+
+ """
+
+ cdef np.ndarray[np.double_t, ndim=2] T
+ cdef np.ndarray[np.double_t, ndim=3] Tinvs
+ cdef int ivertex
+ cdef int i, j, n, nrhs, lda, ldb, info
+ cdef int ipiv[NPY_MAXDIMS+1]
+ cdef int ndim, nvertex
+ cdef double nan
+
+ cdef double x1, x2, x3
+ cdef double y1, y2, y3
+ cdef double det
+
+ nan = np.nan
+ ndim = points.shape[1]
+ nvertex = vertices.shape[0]
+
+ T = np.zeros((ndim, ndim), dtype=np.double)
+ Tinvs = np.zeros((nvertex, ndim+1, ndim), dtype=np.double)
+
+ for ivertex in xrange(nvertex):
+ if ndim == 2:
+ # Manual unrolling of the generic barycentric transform
+ # code below. This is roughly 3.5x faster than the generic
+ # implementation; however, the time taken here is in any
+ # case a small fraction of `qh_new_qhull`, so optimization
+ # here is probably not very important.
+
+ x1 = points[vertices[ivertex,0],0]
+ x2 = points[vertices[ivertex,1],0]
+ x3 = points[vertices[ivertex,2],0]
+
+ y1 = points[vertices[ivertex,0],1]
+ y2 = points[vertices[ivertex,1],1]
+ y3 = points[vertices[ivertex,2],1]
+
+ x1 -= x3
+ x2 -= x3
+
+ y1 -= y3
+ y2 -= y3
+
+ det = x1*y2 - x2*y1
+
+ if det == 0:
+ info = 1
+ else:
+ info = 0
+ Tinvs[ivertex,0,0] = y2/det
+ Tinvs[ivertex,0,1] = -x2/det
+ Tinvs[ivertex,1,0] = -y1/det
+ Tinvs[ivertex,1,1] = x1/det
+ Tinvs[ivertex,2,0] = x3
+ Tinvs[ivertex,2,1] = y3
+ else:
+ # General dimensions
+
+ for i in xrange(ndim):
+ Tinvs[ivertex,ndim,i] = points[vertices[ivertex,ndim],i]
+ for j in xrange(ndim):
+ T[i,j] = (points[vertices[ivertex,j],i]
+ - Tinvs[ivertex,ndim,i])
+ Tinvs[ivertex,i,i] = 1
+
+ n = ndim
+ nrhs = ndim
+ lda = ndim
+ ldb = ndim
+ qh_dgesv(&n, &nrhs, <double*>T.data, &lda, ipiv,
+ (<double*>Tinvs.data) + ndim*(ndim+1)*ivertex, &ldb, &info)
+
+ if info != 0:
+ # degenerate simplex
+ for i in xrange(ndim+1):
+ for j in xrange(ndim):
+ Tinvs[ivertex,i,j] = nan
+
+ return Tinvs
+
+cdef int _barycentric_inside(int ndim, double *transform,
+ double *x, double *c, double eps) nogil:
+ """
+ Check whether point is inside a simplex, using barycentric
+ coordinates. `c` will be filled with barycentric coordinates, if
+ the point happens to be inside.
+
+ """
+ cdef int i, j
+ c[ndim] = 1.0
+ for i in xrange(ndim):
+ c[i] = 0
+ for j in xrange(ndim):
+ c[i] += transform[ndim*i + j] * (x[j] - transform[ndim*ndim + j])
+ c[ndim] -= c[i]
+
+ if not (-eps <= c[i] <= 1 + eps):
+ return 0
+ if not (-eps <= c[ndim] <= 1 + eps):
+ return 0
+ return 1
+
+cdef void _barycentric_coordinate_single(int ndim, double *transform,
+ double *x, double *c, int i) nogil:
+ """
+ Compute a single barycentric coordinate.
+
+ Before the ndim+1'th coordinate can be computed, the other must have
+ been computed earlier.
+
+ """
+ cdef int j
+
+ if i == ndim:
+ c[ndim] = 1.0
+ for j in xrange(ndim):
+ c[ndim] -= c[j]
+ else:
+ c[i] = 0
+ for j in xrange(ndim):
+ c[i] += transform[ndim*i + j] * (x[j] - transform[ndim*ndim + j])
+
+cdef void _barycentric_coordinates(int ndim, double *transform,
+ double *x, double *c) nogil:
+ """
+ Compute barycentric coordinates.
+
+ """
+ cdef int i, j
+ c[ndim] = 1.0
+ for i in xrange(ndim):
+ c[i] = 0
+ for j in xrange(ndim):
+ c[i] += transform[ndim*i + j] * (x[j] - transform[ndim*ndim + j])
+ c[ndim] -= c[i]
+
+
+#------------------------------------------------------------------------------
+# N-D geometry
+#------------------------------------------------------------------------------
+
+cdef void _lift_point(DelaunayInfo_t *d, double *x, double *z) nogil:
+ cdef int i
+ z[d.ndim] = 0
+ for i in xrange(d.ndim):
+ z[i] = x[i]
+ z[d.ndim] += x[i]**2
+ z[d.ndim] *= d.paraboloid_scale
+ z[d.ndim] += d.paraboloid_shift
+
+cdef double _distplane(DelaunayInfo_t *d, int isimplex, double *point) nogil:
+ """
+ qh_distplane
+ """
+ cdef double dist
+ cdef int k
+ dist = d.equations[isimplex*(d.ndim+2) + d.ndim+1]
+ for k in xrange(d.ndim+1):
+ dist += d.equations[isimplex*(d.ndim+2) + k] * point[k]
+ return dist
+
+
+#------------------------------------------------------------------------------
+# Iterating over ridges connected to a vertex in 2D
+#------------------------------------------------------------------------------
+
+cdef void _RidgeIter2D_init(RidgeIter2D_t *it, DelaunayInfo_t *d,
+ int vertex) nogil:
+ """
+ Start iteration over all triangles connected to the given vertex.
+
+ """
+
+ cdef double c[3]
+ cdef int k, ivertex, start
+
+ start = 0
+ it.info = d
+ it.vertex = vertex
+ it.triangle = d.vertex_to_simplex[vertex]
+ it.start_triangle = it.triangle
+
+ if it.triangle != -1:
+ # find some edge connected to this vertex; start from -1 edges
+ for k in xrange(3):
+ ivertex = it.info.vertices[it.triangle*3 + k]
+ if ivertex != vertex:
+ it.vertex2 = ivertex
+ it.edge = k
+ it.start_edge = k
+ break
+ else:
+ it.start_edge = -1
+ it.edge = -1
+
+cdef void _RidgeIter2D_next(RidgeIter2D_t *it) nogil:
+ cdef int itri, k, ivertex
+
+ #
+ # Remember: k-th edge and k-th neigbour are opposite vertex k;
+ # imagine now we are iterating around vertex `O`
+ #
+ # .O------,
+ # ./ |\. |
+ # ./ | \. |
+ # \ | \. |
+ # \ |k \. |
+ # \ | \.|
+ # `+------k
+ #
+
+ # jump to the next triangle
+ itri = it.info.neighbors[it.triangle*3 + it.edge]
+
+ # if it's outside triangulation, restart it to the opposite direction
+ if itri == -1:
+ if it.start_edge == -1:
+ # we already did that -> we have iterated over everything
+ it.edge = -1
+ return
+
+ for k in xrange(3):
+ ivertex = it.info.vertices[it.triangle*3 + k]
+ if ivertex != it.vertex and k != it.start_edge:
+ it.edge = k
+ it.vertex2 = ivertex
+ break
+
+ it.start_edge = -1
+ return
+
+ # Find at which edge we are now:
+ #
+ # it.vertex
+ # O-------k------.
+ # | \- /
+ # | \- E B /
+ # | \- /
+ # | A \- /
+ # +---------´
+ #
+ # A = it.triangle
+ # B = itri
+ # E = it.edge
+ # O = it.vertex
+ #
+ for k in xrange(3):
+ ivertex = it.info.vertices[itri*3 + k]
+ if it.info.neighbors[itri*3 + k] != it.triangle and \
+ ivertex != it.vertex:
+ it.edge = k
+ it.vertex2 = ivertex
+ break
+
+ it.triangle = itri
+
+ # check termination
+ if it.triangle == it.start_triangle:
+ it.edge = -1
+ return
+
+#------------------------------------------------------------------------------
+# Finding simplices
+#------------------------------------------------------------------------------
+
+cdef int _is_point_fully_outside(DelaunayInfo_t *d, double *x,
+ double eps) nogil:
+ """
+ Is the point outside the bounding box of the triangulation?
+
+ """
+
+ cdef int i
+ for i in xrange(d.ndim):
+ if x[i] < d.min_bound[i] - eps or x[i] > d.max_bound[i] + eps:
+ return 1
+ return 0
+
+cdef int _find_simplex_bruteforce(DelaunayInfo_t *d, double *c,
+ double *x, double eps) nogil:
+ """
+ Find simplex containing point `x` by going through all simplices.
+
+ """
+ cdef int inside, isimplex
+
+ if _is_point_fully_outside(d, x, eps):
+ return -1
+
+ for isimplex in xrange(d.nsimplex):
+ inside = _barycentric_inside(
+ d.ndim,
+ d.transform + isimplex*d.ndim*(d.ndim+1),
+ x, c, eps)
+
+ if inside:
+ return isimplex
+ return -1
+
+cdef int _find_simplex_directed(DelaunayInfo_t *d, double *c,
+ double *x, int *start, double eps) nogil:
+ """
+ Find simplex containing point `x` via a directed walk in the tesselation.
+
+ If the simplex is found, the array `c` is filled with the corresponding
+ barycentric coordinates.
+
+ Notes
+ -----
+
+ The idea here is the following:
+
+ 1) In a simplex, the k-th neighbour is opposite the k-th vertex.
+ Call the ridge between them the k-th ridge.
+
+ 2) If the k-th barycentric coordinate of the target point is negative,
+ then the k-th vertex and the target point lie on the opposite sides
+ of the k-th ridge.
+
+ 3) Consequently, the k-th neighbour simplex is *closer* to the target point
+ than the present simplex, if projected on the normal of the k-th ridge.
+
+ 4) In a regular tesselation, hopping to any such direction is OK.
+
+ Also, if one of the negative-coordinate neighbors happens to be -1,
+ then the target point is outside the tesselation (because the
+ tesselation is convex!).
+
+ 5) If all barycentric coordinates are in [-eps, 1+eps], we have found the
+ simplex containing the target point.
+
+ 6) If all barycentric coordinates are non-negative but 5) is not true,
+ we are in an inconsistent situation -- this should never happen.
+
+ """
+ cdef int k, m, ndim, inside, isimplex
+ cdef double *transform
+ cdef double v
+
+ ndim = d.ndim
+ isimplex = start[0]
+
+ if isimplex < 0 or isimplex >= d.nsimplex:
+ isimplex = 0
+
+ while isimplex != -1:
+ transform = d.transform + isimplex*ndim*(ndim+1)
+
+ inside = 1
+ for k in xrange(ndim+1):
+ _barycentric_coordinate_single(ndim, transform, x, c, k)
+
+ if c[k] < -eps:
+ # The target point is in the direction of neighbor `k`!
+
+ m = d.neighbors[(ndim+1)*isimplex + k]
+ if m == -1:
+ # The point is outside the triangulation: bail out
+ start[0] = isimplex
+ return -1
+
+ # Check that the target simplex is not degenerate.
+ v = d.transform[m*ndim*(ndim+1)]
+ if v != v:
+ # nan
+ continue
+ else:
+ isimplex = m
+ inside = -1
+ break
+ elif c[k] > 1 + eps:
+ # we're outside this simplex
+ inside = 0
+
+ if inside == -1:
+ # hopped to another simplex
+ continue
+ elif inside == 1:
+ # we've found the right one!
+ break
+ else:
+ # we've failed utterly (degenerate simplices in the way).
+ # fall back to brute force
+ isimplex = _find_simplex_bruteforce(d, c, x, eps)
+ break
+
+ start[0] = isimplex
+ return isimplex
+
+cdef int _find_simplex(DelaunayInfo_t *d, double *c,
+ double *x, int *start, double eps) nogil:
+ """
+ Find simplex containing point `x` by walking the triangulation.
+
+ Notes
+ -----
+ This algorithm is similar as used by ``qh_findbest``. The idea
+ is the following:
+
+ 1. Delaunay triangulation is a projection of the lower half of a convex
+ hull, of points lifted on a paraboloid.
+
+ Simplices in the triangulation == facets on the convex hull.
+
+ 2. If a point belongs to a given simplex in the triangulation,
+ its image on the paraboloid is on the positive side of
+ the corresponding facet.
+
+ 3. However, it is not necessarily the *only* such facet.
+
+ 4. Also, it is not necessarily the facet whose hyperplane distance
+ to the point on the paraboloid is the largest.
+
+ ..note::
+
+ If I'm not mistaken, `qh_findbestfacet` finds a facet for
+ which the plane distance is maximized -- so it doesn't always
+ return the simplex containing the point given. For example:
+
+ >>> p = np.array([(1 - 1e-4, 0.1)])
+ >>> points = np.array([(0,0), (1, 1), (1, 0), (0.99189033, 0.37674127),
+ ... (0.99440079, 0.45182168)], dtype=np.double)
+ >>> tri = qhull.delaunay(points)
+ >>> tri.vertices
+ array([[4, 1, 0],
+ [4, 2, 1],
+ [3, 2, 0],
+ [3, 4, 0],
+ [3, 4, 2]])
+ >>> dist = qhull.plane_distance(tri, p)
+ >>> dist
+ array([[-0.12231439, 0.00184863, 0.01049659, -0.04714842,
+ 0.00425905]])
+ >>> tri.vertices[dist.argmax()]
+ array([3, 2, 0]
+
+ Now, the maximally positive-distant simplex is [3, 2, 0], although
+ the simplex containing the point is [4, 2, 1].
+
+ In this algorithm, we walk around the tesselation trying to locate
+ a positive-distant facet. After finding one, we fall back to a
+ directed search.
+
+ """
+ cdef int isimplex, i, j, k, inside, ineigh, neighbor_found
+ cdef int ndim
+ cdef double z[NPY_MAXDIMS+1]
+ cdef double best_dist, dist
+ cdef int changed
+
+ if _is_point_fully_outside(d, x, eps):
+ return -1
+ if d.nsimplex <= 0:
+ return -1
+
+ ndim = d.ndim
+ isimplex = start[0]
+
+ if isimplex < 0 or isimplex >= d.nsimplex:
+ isimplex = 0
+
+ # Lift point to paraboloid
+ _lift_point(d, x, z)
+
+ # Walk the tesselation searching for a facet with a positive planar distance
+ best_dist = _distplane(d, isimplex, z)
+ changed = 1
+ while changed:
+ if best_dist > 0:
+ break
+ changed = 0
+ for k in xrange(ndim+1):
+ ineigh = d.neighbors[(ndim+1)*isimplex + k]
+ if ineigh == -1:
+ continue
+ dist = _distplane(d, ineigh, z)
+ if dist > best_dist:
+ # Note: this is intentional: we jump in the middle of the cycle,
+ # and continue the cycle from the next k.
+ #
+ # This apparently sweeps the different directions more
+ # efficiently. We don't need full accuracy, since we do
+ # a directed search afterwards in any case.
+ isimplex = ineigh
+ best_dist = dist
+ changed = 1
+
+ # We should now be somewhere near the simplex containing the point,
+ # locate it with a directed search
+ start[0] = isimplex
+ return _find_simplex_directed(d, c, x, start, eps)
+
+
+#------------------------------------------------------------------------------
+# Delaunay triangulation interface, for Python
+#------------------------------------------------------------------------------
+
+class Delaunay(object):
+ """
+ Delaunay(points)
+
+ Delaunay tesselation in N dimensions
+
+ .. versionadded:: 0.9
+
+ Parameters
+ ----------
+ points : ndarray of floats, shape (npoints, ndim)
+ Coordinates of points to triangulate
+
+ Attributes
+ ----------
+ points : ndarray of double, shape (npoints, ndim)
+ Points in the triangulation
+ vertices : ndarray of ints, shape (nsimplex, ndim+1)
+ Indices of vertices forming simplices in the triangulation
+ neighbors : ndarray of ints, shape (nsimplex, ndim+1)
+ Indices of neighbor simplices for each simplex.
+ The kth neighbor is opposite to the kth vertex.
+ For simplices at the boundary, -1 denotes no neighbor.
+ equations : ndarray of double, shape (nsimplex, ndim+2)
+ [normal, offset] forming the hyperplane equation of the facet
+ on the paraboloid. (See [Qhull]_ documentation for more.)
+ paraboloid_scale, paraboloid_shift : float
+ Scale and shift for the extra paraboloid dimension.
+ (See [Qhull]_ documentation for more.)
+ transform : ndarray of double, shape (nsimplex, ndim+1, ndim)
+ Affine transform from ``x`` to the barycentric coordinates ``c``.
+ This is defined by::
+
+ T c = x - r
+
+ At vertex ``j``, ``c_j = 1`` and the other coordinates zero.
+
+ For simplex ``i``, ``transform[i,:ndim,:ndim]`` contains
+ inverse of the matrix ``T``, and ``transform[i,ndim,:]``
+ contains the vector ``r``.
+ vertex_to_simplex : ndarray of int, shape (npoints,)
+ Lookup array, from a vertex, to some simplex which it is a part of.
+ convex_hull : ndarray of int, shape (nfaces, ndim)
+ Vertices of facets forming the convex hull of the point set.
+ The array contains the indices of the points belonging to
+ the (N-1)-dimensional facets that form the convex hull
+ of the triangulation.
+
+ Notes
+ -----
+ The tesselation is computed using the Qhull libary [Qhull]_.
+
+ References
+ ----------
+
+ .. [Qhull] http://www.qhull.org/
+
+ """
+
+ def __init__(self, points):
+ vertices, neighbors, equations, paraboloid_scale, paraboloid_shift = \
+ _construct_delaunay(points)
+
+ self.ndim = points.shape[1]
+ self.npoints = points.shape[0]
+ self.nsimplex = vertices.shape[0]
+ self.points = points
+ self.vertices = vertices
+ self.neighbors = neighbors
+ self.equations = equations
+ self.paraboloid_scale = paraboloid_scale
+ self.paraboloid_shift = paraboloid_shift
+ self.min_bound = self.points.min(axis=0)
+ self.max_bound = self.points.max(axis=0)
+ self._transform = None
+ self._vertex_to_simplex = None
+
+ @property
+ def transform(self):
+ """
+ Affine transform from ``x`` to the barycentric coordinates ``c``.
+
+ :type: ndarray of double, shape (nsimplex, ndim+1, ndim)
+
+ This is defined by::
+
+ T c = x - r
+
+ At vertex ``j``, ``c_j = 1`` and the other coordinates zero.
+
+ For simplex ``i``, ``transform[i,:ndim,:ndim]`` contains
+ inverse of the matrix ``T``, and ``transform[i,ndim,:]``
+ contains the vector ``r``.
+
+ """
+ if self._transform is None:
+ self._transform = _get_barycentric_transforms(self.points,
+ self.vertices)
+ return self._transform
+
+ @property
+ def vertex_to_simplex(self):
+ """
+ Lookup array, from a vertex, to some simplex which it is a part of.
+
+ :type: ndarray of int, shape (npoints,)
+ """
+ cdef int isimplex, k, ivertex, nsimplex, ndim
+ cdef np.ndarray[np.int_t, ndim=2] vertices
+ cdef np.ndarray[np.int_t, ndim=1] arr
+
+ if self._vertex_to_simplex is None:
+ self._vertex_to_simplex = np.empty((self.npoints,), dtype=int)
+ self._vertex_to_simplex.fill(-1)
+
+ arr = self._vertex_to_simplex
+ vertices = self.vertices
+
+ nsimplex = self.nsimplex
+ ndim = self.ndim
+
+ for isimplex in xrange(nsimplex):
+ for k in xrange(ndim+1):
+ ivertex = vertices[isimplex, k]
+ if arr[ivertex] == -1:
+ arr[ivertex] = isimplex
+
+ return self._vertex_to_simplex
+
+ @property
+ @cython.boundscheck(False)
+ def convex_hull(self):
+ """
+ Vertices of facets forming the convex hull of the point set.
+
+ :type: ndarray of int, shape (nfaces, ndim)
+
+ The array contains the indices of the points
+ belonging to the (N-1)-dimensional facets that form the convex
+ hull of the triangulation.
+
+ """
+ cdef int isimplex, k, j, ndim, nsimplex, m, msize
+ cdef np.ndarray[np.int_t, ndim=2] arr
+ cdef np.ndarray[np.int_t, ndim=2] neighbors
+ cdef np.ndarray[np.int_t, ndim=2] vertices
+
+ neighbors = self.neighbors
+ vertices = self.vertices
+ ndim = self.ndim
+ nsimplex = self.nsimplex
+
+ msize = 10
+ out = np.empty((msize, ndim), dtype=int)
+ arr = out
+
+ m = 0
+ for isimplex in xrange(nsimplex):
+ for k in xrange(ndim+1):
+ if neighbors[isimplex,k] == -1:
+ for j in xrange(ndim+1):
+ if j < k:
+ arr[m,j] = vertices[isimplex,j]
+ elif j > k:
+ arr[m,j-1] = vertices[isimplex,j]
+ m += 1
+
+ if m >= msize:
+ arr = None
+ msize = 2*msize + 1
+ out.resize(msize, ndim)
+ arr = out
+
+ arr = None
+ out.resize(m, ndim)
+ return out
+
+ def find_simplex(self, xi, bruteforce=False):
+ """
+ find_simplex(xi, bruteforce=False)
+
+ Find the simplices containing the given points.
+
+ Parameters
+ ----------
+ tri : DelaunayInfo
+ Delaunay triangulation
+ xi : ndarray of double, shape (..., ndim)
+ Points to locate
+ bruteforce : bool, optional
+ Whether to only perform a brute-force search
+
+ Returns
+ -------
+ i : ndarray of int, same shape as `xi`
+ Indices of simplices containing each point.
+ Points outside the triangulation get the value -1.
+
+ Notes
+ -----
+ This uses an algorithm adapted from Qhull's qh_findbestfacet,
+ which makes use of the connection between a convex hull and a
+ Delaunay triangulation. After finding the simplex closest to
+ the point in N+1 dimensions, the algorithm falls back to
+ directed search in N dimensions.
+
+ """
+ cdef DelaunayInfo_t *info
+ cdef int isimplex
+ cdef double c[NPY_MAXDIMS]
+ cdef double eps
+ cdef int start
+ cdef int k
+ cdef np.ndarray[np.double_t, ndim=2] x
+ cdef np.ndarray[np.int_t, ndim=1] out_
+
+ xi = np.asanyarray(xi)
+
+ if xi.shape[-1] != self.ndim:
+ raise ValueError("wrong dimensionality in xi")
+
+ xi_shape = xi.shape
+ xi = xi.reshape(np.prod(xi.shape[:-1]), xi.shape[-1])
+ x = np.ascontiguousarray(xi.astype(np.double))
+
+ start = 0
+
+ eps = np.finfo(np.double).eps * 10
+ out = np.zeros((xi.shape[0],), dtype=np.int)
+ out_ = out
+ info = _get_delaunay_info(self, 1, 0)
+
+ if bruteforce:
+ for k in xrange(x.shape[0]):
+ isimplex = _find_simplex_bruteforce(
+ info, c,
+ <double*>x.data + info.ndim*k,
+ eps)
+ out_[k] = isimplex
+ else:
+ for k in xrange(x.shape[0]):
+ isimplex = _find_simplex(info, c, <double*>x.data + info.ndim*k,
+ &start, eps)
+ out_[k] = isimplex
+
+ free(info)
+
+ return out.reshape(xi_shape[:-1])
+
+ def plane_distance(self, xi):
+ """
+ plane_distance(xi)
+
+ Compute hyperplane distances to the point `xi` from all simplices.
+
+ """
+ cdef np.ndarray[np.double_t, ndim=2] x
+ cdef np.ndarray[np.double_t, ndim=2] out_
+ cdef DelaunayInfo_t *info
+ cdef double z[NPY_MAXDIMS+1]
+ cdef int i, j, k
+
+ if xi.shape[-1] != self.ndim:
+ raise ValueError("xi has different dimensionality than "
+ "triangulation")
+
+ xi_shape = xi.shape
+ xi = xi.reshape(np.prod(xi.shape[:-1]), xi.shape[-1])
+ x = np.ascontiguousarray(xi.astype(np.double))
+
+ info = _get_delaunay_info(self, 0, 0)
+
+ out = np.zeros((x.shape[0], info.nsimplex), dtype=np.double)
+ out_ = out
+
+ for i in xrange(x.shape[0]):
+ for j in xrange(info.nsimplex):
+ _lift_point(info, (<double*>x.data) + info.ndim*i, z)
+ out_[i,j] = _distplane(info, j, z)
+
+ free(info)
+
+ return out.reshape(xi_shape[:-1] + (self.nsimplex,))
+
+ def lift_points(tri, x):
+ """
+ lift_points(tri, x)
+
+ Lift points to the Qhull paraboloid.
+
+ """
+ z = np.zeros(x.shape[:-1] + (x.shape[-1]+1,), dtype=np.double)
+ z[...,:-1] = x
+ z[...,-1] = (x**2).sum(axis=-1)
+ z[...,-1] *= tri.paraboloid_scale
+ z[...,-1] += tri.paraboloid_shift
+ return z
+
+# Alias familiar from other environments
+def tsearch(tri, xi):
+ """
+ tsearch(tri, xi)
+
+ Find simplices containing the given points. This function does the
+ same thing as Delaunay.find_simplex.
+
+ .. versionadded:: 0.9
+
+ See Also
+ --------
+ Delaunay.find_simplex
+
+ """
+ return tri.find_simplex(xi)
+
+
+#------------------------------------------------------------------------------
+# Delaunay triangulation interface, for low-level C
+#------------------------------------------------------------------------------
+
+cdef DelaunayInfo_t *_get_delaunay_info(obj,
+ int compute_transform,
+ int compute_vertex_to_simplex):
+ cdef DelaunayInfo_t *info
+ cdef np.ndarray[np.double_t, ndim=3] transform
+ cdef np.ndarray[np.int_t, ndim=1] vertex_to_simplex
+ cdef np.ndarray[np.double_t, ndim=2] points = obj.points
+ cdef np.ndarray[np.int_t, ndim=2] vertices = obj.vertices
+ cdef np.ndarray[np.int_t, ndim=2] neighbors = obj.neighbors
+ cdef np.ndarray[np.double_t, ndim=2] equations = obj.equations
+ cdef np.ndarray[np.double_t, ndim=1] min_bound = obj.min_bound
+ cdef np.ndarray[np.double_t, ndim=1] max_bound = obj.max_bound
+
+ info = <DelaunayInfo_t*>malloc(sizeof(DelaunayInfo_t))
+ info.ndim = points.shape[1]
+ info.npoints = points.shape[0]
+ info.nsimplex = vertices.shape[0]
+ info.points = <double*>points.data
+ info.vertices = <int*>vertices.data
+ info.neighbors = <int*>neighbors.data
+ info.equations = <double*>equations.data
+ info.paraboloid_scale = obj.paraboloid_scale
+ info.paraboloid_shift = obj.paraboloid_shift
+ if compute_transform:
+ transform = obj.transform
+ info.transform = <double*>transform.data
+ else:
+ info.transform = NULL
+ if compute_vertex_to_simplex:
+ vertex_to_simplex = obj.vertex_to_simplex
+ info.vertex_to_simplex = <int*>vertex_to_simplex.data
+ else:
+ info.vertex_to_simplex = NULL
+ info.min_bound = <double*>min_bound.data
+ info.max_bound = <double*>max_bound.data
+
+ return info
Added: trunk/scipy/spatial/qhull_blas.h
===================================================================
--- trunk/scipy/spatial/qhull_blas.h (rev 0)
+++ trunk/scipy/spatial/qhull_blas.h 2010-08-31 21:55:50 UTC (rev 6653)
@@ -0,0 +1,19 @@
+/*
+ * Handle different Fortran conventions.
+ */
+
+#if defined(NO_APPEND_FORTRAN)
+#if defined(UPPERCASE_FORTRAN)
+#define F_FUNC(f,F) F
+#else
+#define F_FUNC(f,F) f
+#endif
+#else
+#if defined(UPPERCASE_FORTRAN)
+#define F_FUNC(f,F) F##_
+#else
+#define F_FUNC(f,F) f##_
+#endif
+#endif
+
+#define qh_dgesv F_FUNC(dgesv,DGESV)
Modified: trunk/scipy/spatial/setup.py
===================================================================
--- trunk/scipy/spatial/setup.py 2010-08-31 21:54:55 UTC (rev 6652)
+++ trunk/scipy/spatial/setup.py 2010-08-31 21:55:50 UTC (rev 6653)
@@ -4,10 +4,28 @@
def configuration(parent_package = '', top_path = None):
from numpy.distutils.misc_util import Configuration, get_numpy_include_dirs
+ from numpy.distutils.system_info import get_info
+
config = Configuration('spatial', parent_package, top_path)
config.add_data_dir('tests')
+ qhull_src = ['geom.c', 'geom2.c', 'global.c', 'io.c', 'mem.c',
+ 'merge.c', 'poly.c', 'poly2.c', 'qset.c', 'user.c',
+ 'stat.c', 'qhull.c']
+
+ config.add_library('qhull',
+ sources=[join('qhull', 'src', x) for x in qhull_src],
+ # XXX: GCC dependency!
+ #extra_compiler_args=['-fno-strict-aliasing'],
+ )
+
+ lapack = get_info('lapack_opt')
+ config.add_extension('qhull',
+ sources=['qhull.c'],
+ libraries=['qhull'] + lapack['libraries'],
+ )
+
config.add_extension('ckdtree', sources=['ckdtree.c']) # FIXME: cython
config.add_extension('_distance_wrap',
Added: trunk/scipy/spatial/tests/test_qhull.py
===================================================================
--- trunk/scipy/spatial/tests/test_qhull.py (rev 0)
+++ trunk/scipy/spatial/tests/test_qhull.py 2010-08-31 21:55:50 UTC (rev 6653)
@@ -0,0 +1,158 @@
+import numpy as np
+from numpy.testing import *
+
+import scipy.spatial.qhull as qhull
+
+class TestUtilities(object):
+ """
+ Check that utility functions work.
+
+ """
+
+ def test_find_simplex(self):
+ # Simple check that simplex finding works
+ points = np.array([(0,0), (0,1), (1,1), (1,0)], dtype=np.double)
+ tri = qhull.Delaunay(points)
+
+ # +---+
+ # |\ 0|
+ # | \ |
+ # |1 \|
+ # +---+
+
+ assert_equal(tri.vertices, [[3, 1, 2], [3, 1, 0]])
+
+ for p in [(0.25, 0.25, 1),
+ (0.75, 0.75, 0),
+ (0.3, 0.2, 1)]:
+ i = tri.find_simplex(p[:2])
+ assert_equal(i, p[2], err_msg='%r' % (p,))
+ j = qhull.tsearch(tri, p[:2])
+ assert_equal(i, j)
+
+ def test_plane_distance(self):
+ # Compare plane distance from hyperplane equations obtained from Qhull
+ # to manually computed plane equations
+ x = np.array([(0,0), (1, 1), (1, 0), (0.99189033, 0.37674127),
+ (0.99440079, 0.45182168)], dtype=np.double)
+ p = np.array([0.99966555, 0.15685619], dtype=np.double)
+
+ tri = qhull.Delaunay(x)
+
+ z = tri.lift_points(x)
+ pz = tri.lift_points(p)
+
+ dist = tri.plane_distance(p)
+
+ for j, v in enumerate(tri.vertices):
+ x1 = z[v[0]]
+ x2 = z[v[1]]
+ x3 = z[v[2]]
+
+ n = np.cross(x1 - x3, x2 - x3)
+ n /= np.sqrt(np.dot(n, n))
+ n *= -np.sign(n[2])
+
+ d = np.dot(n, pz - x3)
+
+ assert_almost_equal(dist[j], d)
+
+ def test_convex_hull(self):
+ # Simple check that the convex hull seems to works
+ points = np.array([(0,0), (0,1), (1,1), (1,0)], dtype=np.double)
+ tri = qhull.Delaunay(points)
+
+ # +---+
+ # |\ 0|
+ # | \ |
+ # |1 \|
+ # +---+
+
+ assert_equal(tri.convex_hull, [[1, 2], [3, 2], [1, 0], [3, 0]])
+
+
+class TestTriangulation(object):
+ """
+ Check that triangulation works.
+
+ """
+
+ def test_nd_simplex(self):
+ # simple smoke test: triangulate a n-dimensional simplex
+ for nd in xrange(2, 8):
+ points = np.zeros((nd+1, nd))
+ for j in xrange(nd):
+ points[j,j] = 1.0
+ points[-1,:] = 1.0
+
+ tri = qhull.Delaunay(points)
+
+ tri.vertices.sort()
+
+ assert_equal(tri.vertices, np.arange(nd+1, dtype=np.int)[None,:])
+ assert_equal(tri.neighbors, -1 + np.zeros((nd+1), dtype=np.int)[None,:])
+
+ def test_2d_square(self):
+ # simple smoke test: 2d square
+ points = np.array([(0,0), (0,1), (1,1), (1,0)], dtype=np.double)
+ tri = qhull.Delaunay(points)
+
+ assert_equal(tri.vertices, [[3, 1, 2], [3, 1, 0]])
+ assert_equal(tri.neighbors, [[-1, -1, 1], [-1, -1, 0]])
+
+ def test_duplicate_points(self):
+ x = np.array([0, 1, 0, 1], dtype=np.float64)
+ y = np.array([0, 0, 1, 1], dtype=np.float64)
+
+ xp = np.r_[x, x]
+ yp = np.r_[y, y]
+
+ # shouldn't fail on duplicate points
+ tri = qhull.Delaunay(np.c_[x, y])
+ tri2 = qhull.Delaunay(np.c_[xp, yp])
+
+ pathological_data_1 = np.array([
+ [-3.14,-3.14], [-3.14,-2.36], [-3.14,-1.57], [-3.14,-0.79],
+ [-3.14,0.0], [-3.14,0.79], [-3.14,1.57], [-3.14,2.36],
+ [-3.14,3.14], [-2.36,-3.14], [-2.36,-2.36], [-2.36,-1.57],
+ [-2.36,-0.79], [-2.36,0.0], [-2.36,0.79], [-2.36,1.57],
+ [-2.36,2.36], [-2.36,3.14], [-1.57,-0.79], [-1.57,0.79],
+ [-1.57,-1.57], [-1.57,0.0], [-1.57,1.57], [-1.57,-3.14],
+ [-1.57,-2.36], [-1.57,2.36], [-1.57,3.14], [-0.79,-1.57],
+ [-0.79,1.57], [-0.79,-3.14], [-0.79,-2.36], [-0.79,-0.79],
+ [-0.79,0.0], [-0.79,0.79], [-0.79,2.36], [-0.79,3.14],
+ [0.0,-3.14], [0.0,-2.36], [0.0,-1.57], [0.0,-0.79], [0.0,0.0],
+ [0.0,0.79], [0.0,1.57], [0.0,2.36], [0.0,3.14], [0.79,-3.14],
+ [0.79,-2.36], [0.79,-0.79], [0.79,0.0], [0.79,0.79],
+ [0.79,2.36], [0.79,3.14], [0.79,-1.57], [0.79,1.57],
+ [1.57,-3.14], [1.57,-2.36], [1.57,2.36], [1.57,3.14],
+ [1.57,-1.57], [1.57,0.0], [1.57,1.57], [1.57,-0.79],
+ [1.57,0.79], [2.36,-3.14], [2.36,-2.36], [2.36,-1.57],
+ [2.36,-0.79], [2.36,0.0], [2.36,0.79], [2.36,1.57],
+ [2.36,2.36], [2.36,3.14], [3.14,-3.14], [3.14,-2.36],
+ [3.14,-1.57], [3.14,-0.79], [3.14,0.0], [3.14,0.79],
+ [3.14,1.57], [3.14,2.36], [3.14,3.14],
+ ])
+
+ pathological_data_2 = np.array([
+ [-1, -1 ], [-1, 0], [-1, 1],
+ [ 0, -1 ], [ 0, 0], [ 0, 1],
+ [ 1, -1 - np.finfo(np.float_).eps], [ 1, 0], [ 1, 1],
+ ])
+
+ def test_pathological(self):
+ # both should succeed
+ tri = qhull.Delaunay(self.pathological_data_1)
+ assert_equal(tri.points[tri.vertices].max(),
+ self.pathological_data_1.max())
+ assert_equal(tri.points[tri.vertices].min(),
+ self.pathological_data_1.min())
+
+ tri = qhull.Delaunay(self.pathological_data_2)
+ assert_equal(tri.points[tri.vertices].max(),
+ self.pathological_data_2.max())
+ assert_equal(tri.points[tri.vertices].min(),
+ self.pathological_data_2.min())
+
+if __name__ == "__main__":
+ run_module_suite()
More information about the Scipy-svn
mailing list