root/trunk/delaunay/tests/test_triangulate.py

import sys

from numpy import random
import numpy as np

from scikits import delaunay as dlny

def onright(x0, y0, x1, y1, x, y):
    """Return True if (x,y) is to the right of the vector from (x0,y0) to
    (x1,y1).
    """
    return (y0-y)*(x1-x) > (x0-x)*(y1-y)

def incircle(cx, cy, r, x, y):
    """Return True if (x,y) is strictly inside the circle centered at (cx,cy)
    with radius r.
    """
    r2 = np.hypot(x-cx, y-cy)
    assert r2 < r


class TestSanity(object):
    def setup(self):
        self.rs = random.RandomState(1234567890)

    def test_counts(self):
        for n in (10, 30, 100, 300, 1000, 3000):
            x, y = self.rs.uniform(0, 1, size=(2, n))
            tri = dlny.Triangulation(x, y)
            k = len(tri.hull)
            ntriangles = 2*n - 2 - k
            nedges = 3*n - 3 - k
            assert tri.triangle_nodes.shape == (ntriangles, 3)
            assert tri.triangle_neighbors.shape == (ntriangles, 3)
            assert tri.edge_db.shape == (nedges, 2)
            assert tri.circumcenters.shape == (ntriangles, 2)
            assert np.sum((tri.triangle_neighbors == -1).astype(np.int32).flat) == k

    def test_ccw_triangles(self):
        for n in (10, 30, 100, 300, 1000, 3000):
            x, y = self.rs.uniform(0, 1, size=(2, n))
            tri = dlny.Triangulation(x, y)

            for i,j,k in tri.triangle_nodes:
                assert not onright(x[i], y[i], x[j], y[j], x[k], y[k])

    def test_ccw_hull(self):
        for n in (10, 30, 100, 300, 1000, 3000):
            x, y = self.rs.uniform(0, 1, size=(2, n))
            tri = dlny.Triangulation(x, y)

            hull = list(tri.hull)
            hull.append(hull[0])
            hull.append(hull[1])

            for i,j,k in zip(hull[:-2], hull[1:-1], hull[2:]):
                assert not onright(x[i], y[i], x[j], y[j], x[k], y[k])

    def test_circle_condition(self):
        assert False, "doesn't work perfectly given floating point imprecision"
        for n in (10, 30, 100, 300, 1000, 3000):
            x, y = self.rs.uniform(0, 1, size=(2, n))
            tri = dlny.Triangulation(x, y)

            i = tri.triangle_nodes[:,0]
            r2 = ((x[i] - tri.circumcenters[:,0])**2
                + (y[i] - tri.circumcenters[:,1])**2)
            alldist2 = (np.subtract.outer(x, tri.circumcenters[:,0])**2
                      + np.subtract.outer(y, tri.circumcenters[:,1])**2)
            good = (r2 <= (alldist2 + 1e-12))
            assert np.alltrue(good)


data = 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],
])

class TestRegression(object):

    def test_ticket_376(self):
        x  = np.array([0, 1, 0, 1], dtype=np.float64)
        y  = np.array([0, 0, 1, 1], dtype=np.float64)
        z  = np.array([1, 2, 3, 4], dtype=np.float64)

        tri1 = dlny.Triangulation(x, y)
        nn1 = tri1.nn_interpolator(z)

        xp = np.r_[x, x]
        yp = np.r_[y, y]
        zp = np.r_[z, z]

        # shouldn't fail on duplicate points
        tri2 = dlny.Triangulation(xp, yp)
        nn2 = tri2.nn_interpolator(zp)

        assert nn1(0.3, 0.5) == nn2(0.3, 0.5)

    def test_ticket_376_2(self):
        tri = dlny.Triangulation(data[:,0], data[:,1])
        # should pass

    def test_ticket_382(self):
        x, y = np.meshgrid(np.linspace(0,1,4), np.linspace(0,1,3))

        tri = dlny.Triangulation(x.flatten(), y.flatten())
        refcount1 = sys.getrefcount(tri.x)
        interp = tri.nn_interpolator(x.flatten())
        refcount2 = sys.getrefcount(tri.x)
        interp(x.flatten(), y.flatten())
        refcount3 = sys.getrefcount(tri.x)

        # shouldn't leak references
        assert refcount1 == refcount2
        assert refcount1 == refcount3

if __name__ == '__main__':
    import nose
    nose.main()