[Scipy-svn] r6653 - in trunk/scipy/spatial: . tests

scipy-svn@scip... scipy-svn@scip...
Tue Aug 31 16:55:50 CDT 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
+#------------------------------------------------------------------------------
+
+@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