[Scipy-svn] r6663 - in trunk/scipy: interpolate interpolate/tests spatial spatial/tests
scipy-svn at scipy.org
scipy-svn at scipy.org
Wed Sep 1 16:30:25 EDT 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