[Scipy-svn] r6663 - in trunk/scipy: interpolate interpolate/tests spatial spatial/tests

scipy-svn@scip... scipy-svn@scip...
Wed Sep 1 15:30:25 CDT 2010


Author: ptvirtan
Date: 2010-09-01 15:30:25 -0500 (Wed, 01 Sep 2010)
New Revision: 6663

Modified:
   trunk/scipy/interpolate/interpnd.pyx
   trunk/scipy/interpolate/tests/test_interpnd.py
   trunk/scipy/spatial/qhull.pxd
   trunk/scipy/spatial/qhull.pyx
   trunk/scipy/spatial/tests/test_qhull.py
Log:
BUG: spatial/qhull: fix bugs in RidgeIter2D near boundaries

Also, expose RidgeIter2D to the Python side, so that it can be tested.
Add tests for it.

Modified: trunk/scipy/interpolate/interpnd.pyx
===================================================================
--- trunk/scipy/interpolate/interpnd.pyx	2010-08-31 22:00:43 UTC (rev 6662)
+++ trunk/scipy/interpolate/interpnd.pyx	2010-09-01 20:30:25 UTC (rev 6663)
@@ -394,7 +394,7 @@
             # walk over neighbours of given point
             qhull._RidgeIter2D_init(&it, d, ipoint)
 
-            while it.edge != -1:
+            while it.index != -1:
                 # edge
                 ex = d.points[2*it.vertex2 + 0] - d.points[2*it.vertex + 0]
                 ey = d.points[2*it.vertex2 + 1] - d.points[2*it.vertex + 1]

Modified: trunk/scipy/interpolate/tests/test_interpnd.py
===================================================================
--- trunk/scipy/interpolate/tests/test_interpnd.py	2010-08-31 22:00:43 UTC (rev 6662)
+++ trunk/scipy/interpolate/tests/test_interpnd.py	2010-09-01 20:30:25 UTC (rev 6663)
@@ -98,7 +98,8 @@
         np.random.seed(1234)
         if x is None:
             x = np.array([(0, 0), (0, 1),
-                          (1, 0), (1, 1), (0.25, 0.75), (0.6, 0.8)],
+                          (1, 0), (1, 1), (0.25, 0.75), (0.6, 0.8),
+                          (0.5, 0.2)],
                          dtype=float)
 
         ip = interpnd.CloughTocher2DInterpolator(x, func(x[:,0], x[:,1]),
@@ -128,7 +129,7 @@
             self._check_accuracy(func, tol=1e-13, atol=1e-7, rtol=1e-7,
                                  err_msg="Function %d" % j)
 
-    def test_quadratic_dense_smoketest(self):
+    def test_quadratic_smoketest(self):
         # Should be reasonably accurate for quadratic functions
         funcs = [
             lambda x, y: x**2,
@@ -156,7 +157,7 @@
                      np.random.rand(30*30, 2)]
 
         for j, func in enumerate(funcs):
-            self._check_accuracy(func, x=grid, tol=1e-9, atol=1e-2, rtol=1e-2,
+            self._check_accuracy(func, x=grid, tol=1e-9, atol=5e-3, rtol=1e-2,
                                  err_msg="Function %d" % j)
 
 if __name__ == "__main__":

Modified: trunk/scipy/spatial/qhull.pxd
===================================================================
--- trunk/scipy/spatial/qhull.pxd	2010-08-31 22:00:43 UTC (rev 6662)
+++ trunk/scipy/spatial/qhull.pxd	2010-09-01 20:30:25 UTC (rev 6663)
@@ -77,12 +77,13 @@
 
 ctypedef struct RidgeIter2D_t:
     DelaunayInfo_t *info
+    int index
     int vertex
-    int edge
     int vertex2
     int triangle
     int start_triangle
-    int start_edge
+    int start_index
+    int restart
 
 cdef void _RidgeIter2D_init(RidgeIter2D_t *it, DelaunayInfo_t *d,
                             int vertex) nogil

Modified: trunk/scipy/spatial/qhull.pyx
===================================================================
--- trunk/scipy/spatial/qhull.pyx	2010-08-31 22:00:43 UTC (rev 6662)
+++ trunk/scipy/spatial/qhull.pyx	2010-09-01 20:30:25 UTC (rev 6663)
@@ -485,19 +485,20 @@
     it.vertex = vertex
     it.triangle = d.vertex_to_simplex[vertex]
     it.start_triangle = it.triangle
+    it.restart = 0
 
     if it.triangle != -1:
-        # find some edge connected to this vertex; start from -1 edges
+        # find some edge connected to this vertex
         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
+                it.index = k
+                it.start_index = k
                 break
     else:
-        it.start_edge = -1
-        it.edge = -1
+        it.start_index = -1
+        it.index = -1
 
 cdef void _RidgeIter2D_next(RidgeIter2D_t *it) nogil:
     cdef int itri, k, ivertex
@@ -515,27 +516,48 @@
     #         `+------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:
+    if it.restart:
+        if it.start_index == -1:
             # we already did that -> we have iterated over everything
-            it.edge = -1
+            it.index = -1
             return
 
+        # restart to opposite direction
+        it.triangle = it.start_triangle
         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
+            if ivertex != it.vertex and k != it.start_index:
+                it.index = k
                 it.vertex2 = ivertex
                 break
+        it.start_index = -1
+        it.restart = 0
 
-        it.start_edge = -1
+        if it.info.neighbors[it.triangle*3 + it.index] == -1:
+            it.index = -1
+            return
+        else:
+            _RidgeIter2D_next(it)
+            if it.index == -1:
+                return
+
+    # jump to the next triangle
+    itri = it.info.neighbors[it.triangle*3 + it.index]
+
+    # if it's outside triangulation, take the last edge, and signal
+    # restart to the opposite direction
+    if itri == -1:
+        for k in xrange(3):
+            ivertex = it.info.vertices[it.triangle*3 + k]
+            if ivertex != it.vertex and k != it.index:
+                it.index = k
+                it.vertex2 = ivertex
+                break
+
+        it.restart = 1
         return
 
-    # Find at which edge we are now:
+    # Find at which index we are now:
     #
     # it.vertex
     #      O-------k------.
@@ -547,14 +569,14 @@
     #
     # A = it.triangle
     # B = itri
-    # E = it.edge
+    # E = it.index
     # 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.index = k
             it.vertex2 = ivertex
             break
 
@@ -562,9 +584,39 @@
 
     # check termination
     if it.triangle == it.start_triangle:
-        it.edge = -1
+        it.index = -1
         return
 
+cdef class RidgeIter2D(object):
+    cdef RidgeIter2D_t it
+    cdef object delaunay
+    cdef DelaunayInfo_t *info
+
+    def __init__(self, delaunay, ivertex):
+        self.info = NULL
+        if delaunay.ndim != 2:
+            raise ValueError("RidgeIter2D supports only 2-D")
+        self.delaunay = delaunay
+        self.info = _get_delaunay_info(delaunay, 0, 1)
+        _RidgeIter2D_init(&self.it, self.info, ivertex)
+
+    def __del__(self):
+        if self.info != NULL:
+            free(self.info)
+            self.info = NULL
+        self.delaunay = None
+
+    def __iter__(self):
+        return self
+
+    def __next__(self):
+        if self.it.index == -1:
+            raise StopIteration()
+        ret = (self.it.vertex, self.it.vertex2, self.it.index, self.it.triangle)
+        _RidgeIter2D_next(&self.it)
+        return ret
+
+
 #------------------------------------------------------------------------------
 # Finding simplices
 #------------------------------------------------------------------------------

Modified: trunk/scipy/spatial/tests/test_qhull.py
===================================================================
--- trunk/scipy/spatial/tests/test_qhull.py	2010-08-31 22:00:43 UTC (rev 6662)
+++ trunk/scipy/spatial/tests/test_qhull.py	2010-09-01 20:30:25 UTC (rev 6663)
@@ -70,7 +70,77 @@
 
         assert_equal(tri.convex_hull, [[1, 2], [3, 2], [1, 0], [3, 0]])
 
+class TestRidgeIter2D(object):
 
+    def _check_ridges(self, tri, vertex, expected):
+        got = [(v1, v2) for v1, v2, i, t in qhull.RidgeIter2D(tri, vertex)]
+        got.sort()
+        expected.sort()
+        assert_equal(got, expected, err_msg="%d: %r != %r" % (
+            vertex, got, expected))
+
+    def test_triangle(self):
+        points = np.array([(0,0), (0,1), (1,0)], dtype=np.double)
+        tri = qhull.Delaunay(points)
+
+        # 1
+        # +
+        # |\
+        # | \
+        # |0 \ 
+        # +---+
+        # 0   2
+
+        self._check_ridges(tri, 0, [(0, 1), (0, 2)])
+        self._check_ridges(tri, 1, [(1, 0), (1, 2)])
+        self._check_ridges(tri, 2, [(2, 0), (2, 1)])
+
+    def test_rectangle(self):
+        points = np.array([(0,0), (0,1), (1,1), (1,0)], dtype=np.double)
+        tri = qhull.Delaunay(points)
+
+        # 1   2
+        # +---+
+        # |\ 0|
+        # | \ |
+        # |1 \|
+        # +---+
+        # 0   3
+
+        self._check_ridges(tri, 0, [(0, 1), (0, 3)])
+        self._check_ridges(tri, 1, [(1, 0), (1, 3), (1, 2)])
+        self._check_ridges(tri, 2, [(2, 1), (2, 3)])
+        self._check_ridges(tri, 3, [(3, 0), (3, 1), (3, 2)])
+
+    def test_complicated(self):
+        points = np.array([(0,0), (0,1), (1,1), (1,0),
+                           (0.5, 0.5), (0.9, 0.5)], dtype=np.double)
+        tri = qhull.Delaunay(points)
+
+        #  1                       2
+        #  +-----------------------+
+        #  | \-                 /-||
+        #  |   \-      0      /-  /|
+        #  |     \-         /-   / |
+        #  |       \-     /-    |  |
+        #  |         \-4/-  4  5/  |
+        #  |   1       +-------+  3|
+        #  |         -/  \- 5   \  |
+        #  |      --/      \--   \ |
+        #  |   --/     2      \- | |
+        #  | -/                 \-\|
+        #  +-----------------------+
+        #  0                       3
+        #
+
+        self._check_ridges(tri, 0, [(0, 1), (0, 3), (0, 4)])
+        self._check_ridges(tri, 1, [(1, 0), (1, 2), (1, 4)])
+        self._check_ridges(tri, 2, [(2, 1), (2, 4), (2, 5), (2, 3)])
+        self._check_ridges(tri, 3, [(3, 0), (3, 4), (3, 5), (3, 2)])
+        self._check_ridges(tri, 4, [(4, 0), (4, 1), (4, 2), (4, 3), (4, 5)])
+        self._check_ridges(tri, 5, [(5, 2), (5, 3), (5, 4)])
+
+
 class TestTriangulation(object):
     """
     Check that triangulation works.



More information about the Scipy-svn mailing list