[Scipy-svn] r4654 - in branches/interpolate: . interpSNd

scipy-svn@scip... scipy-svn@scip...
Tue Aug 19 14:12:22 CDT 2008


Author: fcady
Date: 2008-08-19 14:12:18 -0500 (Tue, 19 Aug 2008)
New Revision: 4654

Added:
   branches/interpolate/interpSNd/
   branches/interpolate/interpSNd/_triang.cpp
   branches/interpolate/interpSNd/dewall.py
   branches/interpolate/interpSNd/dewall.pyc
   branches/interpolate/interpSNd/example.py
   branches/interpolate/interpSNd/interp.py
   branches/interpolate/interpSNd/interpolateSNd.py
   branches/interpolate/interpSNd/interpolateSNd.pyc
   branches/interpolate/interpSNd/output.txt
   branches/interpolate/interpSNd/script.py
   branches/interpolate/interpSNd/setup.py
   branches/interpolate/interpSNd/test.py
Log:
added ND scattered interpolation

Added: branches/interpolate/interpSNd/_triang.cpp
===================================================================
--- branches/interpolate/interpSNd/_triang.cpp	2008-08-19 19:07:08 UTC (rev 4653)
+++ branches/interpolate/interpSNd/_triang.cpp	2008-08-19 19:12:18 UTC (rev 4654)
@@ -0,0 +1,85 @@
+/*This will become an extension module used to speed up the slow steps
+    in InterpolateSNd.
+    
+    The best step to speed up is calculation of the circumcircle.
+*/
+
+#include "Python.h"
+#include <stdlib.h>
+#include <map>
+#include <iostream>
+
+#include "numpy/arrayobject.h" // required to use arrays
+#include "numpy/ufuncobject.h"
+
+using namespace std;
+
+extern "C" {
+    
+/* test function */
+PyObject * sayHello(PyObject *self, PyObject *args)
+{
+    int x;
+    PyObject *other_arg;
+    if (!PyArg_ParseTuple(args, "iO", x, &other_arg)){
+        return NULL;
+    }
+    printf("\nhello world\n");
+    return other_arg;
+}
+
+
+
+/* merely return first element of a matrix */
+PyObject * upperleft(PyObject *self, PyObject *args)
+{
+    PyArrayObject *data_coordinates;
+    
+    // parse tuple
+    if (!PyArg_ParseTuple(args, "O!", &PyArray_Type, &data_coordinates)){
+        return NULL;
+    }
+    
+    // make sure floats
+    if (data_coordinates->descr->type_num != PyArray_DOUBLE){
+        PyErr_SetString(PyExc_ValueError, 
+            "array must be of type float");
+        return NULL;
+    }
+    
+    //first element of array
+    double ans;
+    ans = *(double *)(data_coordinates->data);
+    
+    return PyFloat_FromDouble(ans);
+    
+}
+
+
+/* tell Python which methods are available */
+static PyMethodDef triang_methods[]={
+    {"sayHello", (PyCFunction)sayHello, METH_VARARGS, "says hi"},
+    {"upperleft", (PyCFunction)upperleft, METH_VARARGS, ""},
+    //{"lookup", (PyCFunction)lookup, METH_VARARGS, ""},*/
+    {NULL, NULL, 0, NULL}
+};
+
+
+/* initialize the module */
+PyMODINIT_FUNC init_triang(void)
+{
+    PyObject* m;
+    m = Py_InitModule3("_triang", triang_methods,
+        "stupid, useless module that will get better\n"
+    );
+    if (m == NULL)
+        return;
+    import_array(); //required at initialization in order to use arrays
+}
+
+
+    
+    
+    
+    
+} //extern "C"
\ No newline at end of file

Added: branches/interpolate/interpSNd/dewall.py
===================================================================
--- branches/interpolate/interpSNd/dewall.py	2008-08-19 19:07:08 UTC (rev 4653)
+++ branches/interpolate/interpSNd/dewall.py	2008-08-19 19:12:18 UTC (rev 4654)
@@ -0,0 +1,393 @@
+# implements algorithm found at 
+# http://vcg.isti.cnr.it/publications/papers/dewall.pdf
+
+# WARNING
+# This code is grotesquely inefficient and messy.
+# The goal is to have it be 1) correct, and 2) understandable
+# The next phase is to rewrite it using arrays and efficient algorithms.
+# After that, the slow parts should be translated into C.
+
+# In particular, calculation of the circumcircle is a purely
+# mathematical operation that really should be made into C.
+
+import numpy as np
+from numpy.linalg import norm
+
+eps =1e-5
+compare_first_elem = lambda t1, t2 : cmp(t1[0], t2[0])
+point_in_list = lambda p, L : np.any([ (p==elem).all() for elem in L ])
+face_in_list = lambda f, L :  np.any([ np.alltrue([point_in_list(p,face) for p in f]) for face in L])
+compare_pointlists = lambda L1, L2 : np.alltrue([ point_in_list(l1, L2) for l1 in L1]) and \
+                                                    np.alltrue([ point_in_list(l2, L1) for l2 in L2])
+
+def dewall (P, #set of points 
+                AFL = [], # list of faces: (d-1)face list
+                ):
+                
+    # checking input
+    assert isinstance(P, list)
+    if len(P)>0:
+        assert isinstance(P[0], np.ndarray)
+    assert isinstance(AFL, list)
+    if len(AFL)>0: 
+        #print "AFL[0]: ", AFL[0]
+        assert isinstance(AFL[0],tuple)
+        assert isinstance(AFL[0][0], list)
+        assert isinstance(AFL[0][0][0], np.ndarray)
+        assert isinstance(AFL[0][1], np.ndarray)
+    
+    # base case
+    if len(P) == 0: return []
+    if len(P) <= len(P[0]): return [] # better to return [] whenever points are co-hyper-planar
+    
+    # lists of active faces
+    # elem = ( [list of d points] , outvec )
+    AFL_alpha = []
+    AFL1 = []
+    AFL2 = []
+
+    # list of simplices
+    Sigma= []
+    
+    alpha = select_alpha(P, AFL)
+                            
+    # divide points into two sets separated by alpha
+    P1, P2 = pointset_partition(P, alpha) # both lists of points
+    
+    # Simplex Wall Construction
+    just_starting = False #source of problem?
+    if len(AFL) == 0:
+        just_starting = True
+        first_simplex = make_first_simplex(P, alpha)
+        AFL = [ (face, get_out_vec(face,first_simplex))\
+                for face in faces(first_simplex)] # d+1 of them
+        Sigma.append(first_simplex)
+    for face, outvec in AFL:
+        if is_intersected(face, alpha): # not counting as intersected
+            AFL_alpha.append((face, \
+                            get_out_vec(face,first_simplex) if just_starting \
+                            else outvec))
+        if is_subset(face, P1):
+            AFL1.append((face,outvec))
+        if is_subset(face, P2):
+            AFL2.append((face,outvec))
+    while len(AFL_alpha) != 0:
+        
+        face, outvec = AFL_alpha.pop()
+        
+        if outvec is not None:
+            outward_points = filter( lambda p: (np.dot(p,outvec)>np.dot(face[0],outvec)),\
+                                            P)
+        else:
+            outward_points = filter( lambda p: np.all([not point_in_list(p,vertex) for vertex in face]), P)
+                
+        t = make_simplex(face, outward_points) # make only over outer half space
+        if t is not None:
+            Sigma.append(t)
+            # for f0 != f in faces(t) , ie for new outward faces
+            for f0 in filter(lambda f: not face_in_list(f,[face]), faces(t)):
+                new_pair = (f0, get_out_vec(f0, t))
+                if is_intersected(f0, alpha):
+                    # continue building wall out in that direction
+                    AFL_alpha = update(new_pair, AFL_alpha)
+                if is_subset(f0, P1):
+                    AFL1 = update(new_pair, AFL1)
+                if is_subset(f0, P2):
+                    AFL2 = update(new_pair, AFL2)
+                    
+    # now Sigma contains all simplices that intersect alpha
+    
+    # Recursive Triangulation
+    if len(AFL2)!=0:
+        Sigma = Sigma + dewall(P2,AFL2)
+    if len(AFL1)!=0:
+        Sigma = Sigma + dewall(P1,AFL1) # not doing this branch of recursion
+    return Sigma
+
+def select_alpha(P, AFL):
+    # dividing plane.  This must divide points into 2 non-empty parts
+    # must intersect active face if there is one
+    # must not intersect any points
+    
+    #assert isinstance(P, list)
+    #if len(P)>0:
+    #    assert isinstance(P[0], np.ndarray)
+    #assert isinstance(AFL, list)
+    #if len(AFL)>0: 
+    #    assert isinstance(AFL[0],list)
+    #    assert isinstance(AFL[0][0], list)
+    #    assert isinstance(AFL[0][0][0], np.ndarray)
+    #    assert isinstance(AFL[0][1], np.ndarray)
+    
+    d = len(P[0])
+    
+    # plane through avg of cluster.  Guarantees separation
+    # must also make sure intersects face in AFL
+    if len(AFL) != 0:
+        mid = sum(AFL[0][0])/(d)
+    else:
+        mid = sum(P)/len(P)
+    if norm(mid)==0:
+        direction = np.random.random_sample((d))
+        alpha = ( direction/norm(direction), 0)
+    else: 
+        alpha =(mid/norm(mid), norm(mid))
+        
+    return alpha
+
+def get_out_vec(face, simplex): #SEEMS GOOD
+    # face is a face of simplex
+    # returns vector normal to face pointing out of simplex
+    assert len(face)==len(face[0]), "face needs d points"
+    assert len(simplex)-len(face)==1, "simplex needs one more point"
+    
+    d = len(face[0]) # dimension
+    
+    # get point in simplex that's not in face
+    # vector must face away from it
+    other_point = filter(lambda p: not point_in_list(p, face),
+                                simplex)[0]
+    
+    # determinants of submatrices give components of a vector
+    # that is orthogonal to each row
+    vecs_along_face = np.array([point-face[0] for point in face[1:]])
+    normal_vector = np.zeros(d)
+    
+    # normal vector must be orthogonal to each element of matrix_of_diffs
+    # components given by determinants of submatrices
+    for i in range(d):
+        normal_vector[i] = ((-1)**i)*\
+                np.linalg.det(vecs_along_face[:,np.arange(d)!=i])
+    unit_normal = normal_vector/norm(normal_vector)
+    
+    # return unit_normal point in correct direction
+    if np.dot(other_point, unit_normal) < np.dot(face[0], unit_normal):
+        return unit_normal
+    else:
+        return -1*unit_normal
+
+def is_subset(S1, S2):
+    # both are lists of arrays
+    return np.alltrue([ point_in_list(s1, S2) for s1 in S1])
+
+def update (face_pair, face_pair_list):
+    # returns face_list with face_pair added if it wasn't there, else deleted
+    face, outvec = face_pair
+    face_list = [face for face, outvec in face_pair_list]
+    if face_in_list(face, face_pair_list):
+        f_not_equal_face = lambda f :  not np.alltrue([ point_in_list(p, f[0]) for p in face ])
+        face_pair_list = filter(f_not_equal_face, face_pair_list)
+    else:
+        face_pair_list.append(face_pair)
+    return face_pair_list
+        
+def pointset_partition(P, alpha): #WORKS
+    P1 = [p for p in P if np.dot(p,alpha[0])<  alpha[1]]
+    P2 = [p for p in P if np.dot(p,alpha[0])>=alpha[1]]
+    return P1, P2
+    
+def is_intersected(f, alpha): #WORKS
+    assert isinstance(f, list), "the face must be a list: "+str(f)
+    assert isinstance(alpha, tuple), "alpha must be tuple: "+str(alpha)
+    assert isinstance(alpha[0],np.ndarray), "normal vector: "+str(alpha[0])
+    assert isinstance(f[0],np.ndarray), "point in f: "+str(f[0])
+    list_of_dot_prods = [np.dot(alpha[0], point) for point in f]
+    all_dots_nonneg = [ prod >= alpha[1] for prod in list_of_dot_prods]
+    all_dots_nonpos = [ prod <= alpha[1] for prod in list_of_dot_prods]
+    return not (np.all(all_dots_nonneg) or np.all(all_dots_nonpos))
+    
+def faces(simplex): #WORKS
+    # given simplex (as list of points) returns list of faces (each face a list of points)
+    faces = []
+    for point in simplex:
+        faces.append( [p for p in simplex if (p!=point).any()] )
+    return faces
+
+def make_first_simplex(P, alpha): #WORKS
+    # alpha = unit_normal_vec, distance
+    # assume no points on plane
+
+    d = len(P[0])
+    unit_normal, distance = alpha
+    points_and_distances = [(np.dot(unit_normal, point), point) for point in P]
+
+    points_and_dist_from_alpha = [ (norm(distance-dist), point) for (dist, point) in points_and_distances]
+
+    first_point = sorted(points_and_dist_from_alpha,
+                                cmp = compare_first_elem \
+                                )[0][1] # closest to alpha
+                                                    
+    possible_second_pts = [ (circumcircle([first_point, point])[1], point) \
+                                        for point in P if (point != first_point).any() and \
+                                        (np.dot(unit_normal, first_point)-distance)*(np.dot(unit_normal, point)-distance)<=0 \
+                                    ]
+    second_point = sorted(possible_second_pts, \
+                                    cmp = compare_first_elem \
+                                    )[0][1]
+    simplex = [first_point, second_point]
+    for i in range(d-1):
+        radii_of_circumcircles = [(circumcircle(  copy_list(simplex)+[point.copy()]  )[1], point) \
+                                        for point in P if not point_in_list(point, simplex) ]
+        new_point = sorted(radii_of_circumcircles, \
+                                    cmp = compare_first_elem \
+                                    )[0][1]
+        simplex.append(new_point)
+
+    return simplex
+    
+def copy_list(list_of_arrays): #WORKS
+    # returns a list of copies of the arrays
+    # use if want to modify an array list without 
+    result = [array.copy() for array in list_of_arrays]
+    return result
+    
+
+        
+def make_simplex(f,P): #WORKS
+    # must be only in the outer halfspace
+    
+    # returns the simlex
+    
+    # clean up by making sure only points not in hyperplane(f) are in P
+        
+    valid_points = [p for p in P if not point_in_list(p,f)]
+    if len(valid_points) == 0: 
+        return None
+    delaunay_distances = [(delaunay_dist(f,p) ,p) for p in valid_points]
+    new_point = sorted(delaunay_distances, cmp = compare_first_elem)[0][1]
+    simplex = f+[new_point]
+    return simplex
+
+def delaunay_dist(f, p): #WORKS
+    # |distance| = radius of circumcircles
+    # sign depends on whether center is on same side as p
+    
+    # FIXME : don't let pathological stuff come in
+    
+    normal, distance = spanning_hyperplane(f)
+    center, radius = circumcircle( [elem.copy() for elem in f]+[p] ) #need copy of f and add to it
+    if (np.dot(normal,center) >= distance and np.dot(normal,p) >= distance) \
+        or (np.dot(normal,center) <= distance and np.dot(normal,p) <= distance) \
+        or radius == np.Inf:
+        return radius
+    else:
+        return -1*radius
+        
+def spanning_hyperplane(f): #WORKS
+    # f=list of points
+    # return outward unit normal to hyperplane they define
+    # and distance to hyperplane from the origin
+    d = len(f[0])
+    
+    # determinants of submatrices give components of a vector
+    # that is orthogonal to each row
+    matrix_of_diffs = np.array([point-f[0] for point in f[1:]])
+    normal_vector = np.zeros(d)
+    
+    # normal vector must be orthogonal to each element of matrix_of_diffs
+    # components given by determinants of submatrices
+    for i in range(d):
+        normal_vector[i] = ((-1)**i)*np.linalg.det(\
+                                matrix_of_diffs[:,np.arange(d)!=i])
+    unit_normal = normal_vector/norm(normal_vector)
+    distance = np.dot(unit_normal, f[0])
+    
+    # want a positive distance
+    if distance < 0:
+        return -1*unit_normal, -1*distance
+    return unit_normal, distance
+    
+def circumcircle(Points): #WORKS
+    # P = list of distinct points
+    # n distinct points, with n<= d+1
+    # returns (center=array, radius)
+    
+    # FIXME : account for collinear/same and error check
+    
+    # calculates center by forming system of linear constraints
+    # on its components, and calculates radius after that
+    
+    d = len(Points[0]) # dimensionality of space
+    n = len(Points)
+    
+    # need "infinitely big" circle in this case
+    if not linearly_independent(Points):
+        return np.Inf*np.ones(d), np.Inf
+    
+    matrix = np.zeros((d,d), dtype=np.float64) # d constraints to find the center
+    vector = np.zeros((d,1), dtype=np.float64)
+    
+    difference_vecs = [(Points[i]-Points[0])/norm(Points[i]-Points[0]) \
+                                for i in range(1,n) if norm(Points[i]-Points[0])>eps] # n-1 difference vectors.  Point along hyperplane
+    n = len(difference_vecs) + 1 # correct for possible degeneracies in the input
+        
+    # for each P[0] -- P[i] line segment, center is along the middle of it
+    for i in range(n-1):
+        matrix[i,:] = difference_vecs[i]
+        vector[i,0] = np.dot( difference_vecs[i] , (Points[0]+Points[i+1])/2)
+    
+    # center is in the hyperplane defined by the points
+    # ie its componeny orthogonal to difference_vecs is same as that of all points
+    orthog_comp = orthogonal_complement(difference_vecs)
+    m = len(orthog_comp) # m+n = d
+    if m+(n-1) != d:
+        assert False, "messed up.  m=%i, and n=%i " % (m, n)
+    for i in range(m):
+        matrix[n-1+i,:] = orthog_comp[i]
+        vector[n-1+i,0] = np.dot(Points[0], orthog_comp[i])
+    
+    # calculating and returning
+    center = np.linalg.solve(matrix, vector)
+    radius = norm(center-np.reshape(Points[0],(d,1)))
+    return np.reshape(center,(d)), radius
+    
+def orthogonal_complement(P): #WORKS
+    # return orthogonal complement of P of unit length
+    # in = list of 1D arrays
+    # out = list of 1D arrays
+    
+    # nothing here can be pathological
+    
+    P = copy_list(P) # this way won't harm original P
+    d = len(P[0]) # dimensionality
+    n = len(P) # number of points.  <= d
+    
+    # Graham-Schmidt orthonormalization of input vectors
+    for j in range(n):
+        for k in range(j):
+            P[j] = P[j] - np.dot(P[k],P[j])*P[k]
+        if norm(P[j]) > eps:
+            P[j]=P[j]/norm(P[j])
+    
+    # remove linear dependencies and normalize
+    P = np.array([point/norm(point) for point in P if norm(point)>eps])
+    
+    orthog = np.eye(d,d)
+    for j in range(d):
+        # remove component of Oj in P
+        for k in range(len(P)):
+            orthog[j,:] = orthog[j,:] - np.dot(P[k,:],orthog[j,:])*P[k,:]
+            
+        # remove component of Oj in Oi with i<j
+        for m in range(j):
+            orthog[j,:] = orthog[j,:] - np.dot(orthog[m,:],orthog[j,:])*orthog[m,:]
+            
+        # normalize if vector not linearly dependent on previous
+        if norm(orthog[j,:]) > eps:
+            orthog[j,:] = orthog[j,:]/norm(orthog[j,:])
+    
+    # return each non-zero vector
+    return [np.reshape(orthog[j,:],(d)) for j in range(d) if norm(orthog[j,:])>eps]
+    
+def linearly_independent(P):
+    # return true if the vectors are linearly independent
+    # needs n points
+    d=len(P[0])
+    if len(P)>d: return False
+    matrix = np.array(P).reshape((d,len(P)))
+    return np.rank(matrix)==len(P)
+    
+    
+    
+    
+    
\ No newline at end of file

Added: branches/interpolate/interpSNd/dewall.pyc
===================================================================
(Binary files differ)


Property changes on: branches/interpolate/interpSNd/dewall.pyc
___________________________________________________________________
Name: svn:mime-type
   + application/octet-stream

Added: branches/interpolate/interpSNd/example.py
===================================================================
--- branches/interpolate/interpSNd/example.py	2008-08-19 19:07:08 UTC (rev 4653)
+++ branches/interpolate/interpSNd/example.py	2008-08-19 19:12:18 UTC (rev 4654)
@@ -0,0 +1,12 @@
+# example sript
+import _triang as _t
+import numpy as np
+
+print "dir(_t): ", dir(_t)
+a = _t.sayHello(5 ,'foo')
+print "it worked?  I won't see this"
+print "a equals: ", a
+
+arr = np.array([1.0, 2.3, 4.8])
+print "other answer: ", _t.triangulate(arr)
+print "now at end"
\ No newline at end of file

Added: branches/interpolate/interpSNd/interp.py
===================================================================
--- branches/interpolate/interpSNd/interp.py	2008-08-19 19:07:08 UTC (rev 4653)
+++ branches/interpolate/interpSNd/interp.py	2008-08-19 19:12:18 UTC (rev 4654)
@@ -0,0 +1,18 @@
+import interpolateSNd as SN
+import numpy as np
+
+reload(SN)
+
+points = np.array([[ 1.,0.,1.,0.],[0.,0.,1.,1.]])
+z = np.array([1.,0.,2.,1.])
+interp = SN.InterpolateSNd(points,z)
+
+print "triangulation:\n",interp._triangulation
+
+X=np.array([[1.4,.1,.55],[.4,.1,.3]])
+print "X values:\n",X
+# last component is .1 too much
+
+output = interp(X)
+
+print "output:\n", output
\ No newline at end of file

Added: branches/interpolate/interpSNd/interpolateSNd.py
===================================================================
--- branches/interpolate/interpSNd/interpolateSNd.py	2008-08-19 19:07:08 UTC (rev 4653)
+++ branches/interpolate/interpSNd/interpolateSNd.py	2008-08-19 19:12:18 UTC (rev 4654)
@@ -0,0 +1,200 @@
+""" ND Scattered interpolation
+"""
+
+import dewall as dw
+import numpy as np
+
+class InterpolateSNd:
+    """ Interpolation of scatter data in arbitrarily many dimensions
+        
+        *******
+        WARNING : this code is currently extremely slow.  It is very much still
+        *******
+                        in development, but because it functions correctly it is being
+                        presented.
+                        
+        Parameters:
+        ------------
+        
+            points : 2D numpy array
+                    dxn array, where d is dimensionality of the data and
+                    there are n data points.  Each column of the array is
+                    a data point.
+                    
+            fvals : 1D numpy array
+                    Must have length n.  fvals[j] is the value of the function
+                    at point point[:,j]
+                    
+            newX (when calling) : 2D numpy array
+                    shape dxm.  Each column is a point at which to interpolate
+                    the function value.
+                    
+        The valid interpolation range is the convex hull of the points.  All in-range
+        points will be linearly interolated, and all out-of-bounds points will return
+        NaN.
+    """
+    
+    def __init__(self, P, fvals):
+        # P = array of points, where each column is a point
+        #       points must be distinct
+        # fvals = array of known function values
+        
+        assert P.ndim == 2, "P must be 2-dimensional"
+        d, n = P.shape
+        assert len(fvals)==n, "fvals must have length n,\
+                    where n is number of points"
+        
+        # remember dimensionality of space
+        self.ndim = d
+                    
+        # store list of points
+        self.points = P.copy()
+        
+        # store function values
+        self.fvals = fvals.copy()
+        
+        # calculate and store delaunay triangulation
+        list_of_points = [P[:,i].reshape(d) for i in range(n)]
+        self._triangulation = dw.dewall(list_of_points)
+        
+        # for each simplex, indices of its points in self.points
+        # lists of tuples(lists) of integers
+        indices_list = []
+        for simplex in self._triangulation:
+            indices_for_current_point = []
+            for point in simplex:
+                for i in range(n):
+                    # if ith col = point, append i to list
+                    if (point.reshape((d)) == np.reshape(P[:,i],(d))).all():
+                        indices_for_current_point.append(i)
+            indices_list.append(indices_for_current_point)
+        self._indices_list = indices_list
+        
+        # FIXME : store list of (simplex, vals) tuples ????
+            
+    def __call__(self, X):
+        # X = dxm array of points at which to interpolate
+        d, m = X.shape
+        assert d==self.ndim, "input must have correct dimensionality"
+        
+        # break input into a list of points for list comprehensions.
+        # yeah, it's very slow, but it's elegant for now
+        list_of_points = [X[:,i].reshape((d,1)) for i in range(m)]
+                                        
+        #  Tuples of form (simplex , vals_at_vertices, point)
+        simplices_vals_points = [(self._get_simplex(pt), self._get_fvals(pt), pt) \
+                                        for pt in list_of_points]
+        
+        result = [self.linear_interp( simplex, vals, point ) \
+                        for simplex, vals, point in simplices_vals_points]
+                            
+        return result
+        
+    def _get_fvals(self, pt):
+        assert pt.shape == (self.ndim, 1), "point shape: "+str(pt.shape)
+        indices = self._get_simplex_indices(pt)
+        if indices==None:
+            return None
+        return self.fvals[indices]
+                
+    def _get_simplex(self, pt):
+        assert pt.shape == (self.ndim, 1), "point shape: "+str(pt.shape)
+        indices = self._get_simplex_indices(pt)
+        if indices==None:
+            return None
+        return self.points[:,indices]
+        
+    def _get_simplex_indices(self, pt):
+        # returns tuple indicating vertices of simplex containing pt
+        assert pt.shape == (self.ndim, 1), "point shape: "+str(pt.shape)
+        indices = None
+        for i in range(len(self._indices_list)):
+            if self.point_is_in_simplex(pt.reshape((self.ndim,1)), self.points[:,self._indices_list[i]]):
+                indices = self._indices_list[i]
+                break
+        return indices
+        
+    def point_is_in_simplex(self, point, simplex):
+        # point = array
+        # simplex = matrix
+        assert point.shape == (self.ndim, 1), "wrong shape of point: "+str(point.shape)
+        assert simplex.shape == (self.ndim, self.ndim+1), "wrong shape of simplex: "+str(simplex.shape)
+        weights_vec = self.calculate_weights(simplex, point)
+        weight_in_range = (0<= weights_vec) & (weights_vec <= 1)
+        return np.alltrue(weight_in_range, axis=0)
+    
+    def linear_interp(self, simplex, vals, point):
+        if simplex == None or vals == None: # point is out of range
+            return np.NaN
+        
+        # testing
+        assert point.shape == (self.ndim, 1), \
+                "wrong shape of point: "+str(point.shape)
+        assert simplex.shape == (self.ndim, self.ndim+1), \
+                "wrong shape of simplex: "+str(simplex.shape)
+        assert vals.shape == (self.ndim+1,), \
+                "vals wrong shape: "+str(vals.shape)+"need: "+str((self.ndim,))
+        
+        weights = self.calculate_weights(simplex, point)
+        return np.dot(np.ravel(weights), np.ravel(vals))
+    
+    def calculate_weights(self, simplex, points):
+        """ Each column in points is a weighted average
+            of columns in simplex.  Returns matrix where
+            jth column gives these weights for the jth column
+            of points.
+        """
+        assert simplex.shape == (self.ndim, self.ndim+1), "simplex shape: "+str(simplex.shape)
+        d, V = simplex.shape
+        d, P = points.shape
+        
+        matrix_to_solve = np.concatenate((simplex, \
+                                                        np.ones((1,V))
+                                                        ),
+                                                        axis = 0
+                                                        )
+        vec_to_solve = np.concatenate((points, \
+                                                        np.ones((1,P))
+                                                        ),
+                                                        axis = 0
+                                                        )
+        
+        weights_vecs = np.linalg.solve(matrix_to_solve, vec_to_solve)
+        
+        return weights_vecs
+        
+def point_is_in_simplex(point, simplex):
+    # point = array
+    # simplex = matrix
+    #assert point.shape == (self.ndim, 1), "wrong shape of point: "+str(point.shape)
+    #assert simplex.shape == (self.ndim, self.ndim+1), "wrong shape of simplex: "+str(simplex.shape)
+    weights_vec = calculate_weights(simplex, point)
+    print "weights:\n", weights_vec
+    weight_in_range = (0.<= weights_vec) & \
+                                    (weights_vec <= 1.)
+    print "in_range:\n", weight_in_range
+    return np.alltrue(weight_in_range, axis=0)
+    
+def calculate_weights(simplex, points):
+        """ Each column in points is a weighted average
+            of columns in simplex.  Returns matrix where
+            jth column gives these weights for the jth column
+            of points.
+        """
+        N, V = simplex.shape
+        N, P = points.shape
+        
+        matrix_to_solve = np.concatenate((simplex, \
+                                                        np.ones((1,V))
+                                                        ),
+                                                        axis = 0
+                                                        )
+        vec_to_solve = np.concatenate((points, \
+                                                        np.ones((1,P))
+                                                        ),
+                                                        axis = 0
+                                                        )
+        
+        weights_vecs = np.linalg.solve(matrix_to_solve, vec_to_solve)
+        
+        return weights_vecs
\ No newline at end of file

Added: branches/interpolate/interpSNd/interpolateSNd.pyc
===================================================================
(Binary files differ)


Property changes on: branches/interpolate/interpSNd/interpolateSNd.pyc
___________________________________________________________________
Name: svn:mime-type
   + application/octet-stream

Added: branches/interpolate/interpSNd/output.txt
===================================================================
--- branches/interpolate/interpSNd/output.txt	2008-08-19 19:07:08 UTC (rev 4653)
+++ branches/interpolate/interpSNd/output.txt	2008-08-19 19:12:18 UTC (rev 4654)
@@ -0,0 +1,40 @@
+
+face:  [array([ 0.78256165,  1.25794206]), array([ 5.5867473 ,  0.30561524])]
+points:  []
+
+face:  [array([ 5.5867473 ,  0.30561524]), array([ 0.84738355,  0.4417885 ])]
+points:  [array([ 0.49630621,  0.38945595])]
+delaunay distances:
+[(14.481574504701088, array([ 0.49630621,  0.38945595]))]
+
+face:  [array([ 5.5867473 ,  0.30561524]), array([ 0.49630621,  0.38945595])]
+points:  []
+
+face:  [array([ 0.78256165,  1.25794206]), array([ 0.84738355,  0.4417885 ])]
+points:  [array([ 0.45697308,  1.26267539]), array([ 0.49630621,  0.38945595]), array([ 0.15363154,  0.794239  ])]
+delaunay distances:
+[(0.45650505613020009, array([ 0.45697308,  1.26267539])), (0.45830428089834485, array([ 0.49630621,  0.38945595])), (0.45808830610514073, array([ 0.15363154,  0.794239  ]))]
+
+face:  [array([ 0.84738355,  0.4417885 ]), array([ 0.45697308,  1.26267539])]
+points:  [array([ 0.49630621,  0.38945595]), array([ 0.15363154,  0.794239  ])]
+delaunay distances:
+[(0.45691821155399881, array([ 0.49630621,  0.38945595])), (0.45699647524242137, array([ 0.15363154,  0.794239  ]))]
+
+face:  [array([ 0.45697308,  1.26267539]), array([ 0.49630621,  0.38945595])]
+points:  [array([ 0.15363154,  0.794239  ])]
+delaunay distances:
+[(-0.45659643493482976, array([ 0.15363154,  0.794239  ]))]
+
+face:  [array([ 0.45697308,  1.26267539]), array([ 0.15363154,  0.794239  ])]
+points:  []
+
+face:  [array([ 0.84738355,  0.4417885 ]), array([ 0.49630621,  0.38945595])]
+points:  []
+
+face:  [array([ 0.84738355,  0.4417885 ]), array([ 0.49630621,  0.38945595])]
+points:  [array([ 0.15363154,  0.794239  ])]
+delaunay distances:
+[(0.45765190740816952, array([ 0.15363154,  0.794239  ]))]
+
+face:  [array([ 0.84738355,  0.4417885 ]), array([ 0.15363154,  0.794239  ])]
+points:  []

Added: branches/interpolate/interpSNd/script.py
===================================================================
--- branches/interpolate/interpSNd/script.py	2008-08-19 19:07:08 UTC (rev 4653)
+++ branches/interpolate/interpSNd/script.py	2008-08-19 19:12:18 UTC (rev 4654)
@@ -0,0 +1,49 @@
+# example script
+
+import dewall as dw
+import numpy as np
+from numpy import array
+import matplotlib.pyplot as pyplot
+
+Pa=[array([1.1,1.1]), array([-1.,1.]), array([-1.,-1.]), \
+        array([1.,-1.]), array([1.5,1.5])]
+
+def segments(T):
+    seg0 = [ [T[0][0],T[1][0]] , [T[0][1],T[1][1]] ]
+    seg1 = [ [T[0][0],T[2][0]] , [T[0][1],T[2][1]] ]
+    seg2 = [ [T[1][0],T[2][0]] , [T[1][1],T[2][1]] ]
+    return [seg0, seg1, seg2]
+    
+def plot_circle(circle):
+    center, radius = circle
+    t = np.linspace(0,2*np.pi,500)
+    x_offset = radius*np.sin(t)
+    y_offset = radius*np.cos(t)
+    x = center[0] + x_offset
+    y = center[1] + y_offset
+    pyplot.plot(x,y)
+    
+Pb=[array([.25, -.25]), array([0,.75])]
+
+P = Pa+Pb
+
+P = [ np.array([np.random.gamma(1), np.random.gamma(1)]) \
+        for j in range(6) ]
+
+triangul = dw.dewall(P)
+
+# plotting the known data points
+for p in P: pyplot.scatter([p[0]],[p[1]])
+
+# plotting the triangulation
+for tri in triangul:
+    for seg in segments(tri):
+        pyplot.plot(seg[0],seg[1])
+        
+# plotting the circumcircles
+for tri in triangul:
+    plot_circle(dw.circumcircle(tri))
+        
+pyplot.show()
+
+print triangul
\ No newline at end of file

Added: branches/interpolate/interpSNd/setup.py
===================================================================
--- branches/interpolate/interpSNd/setup.py	2008-08-19 19:07:08 UTC (rev 4653)
+++ branches/interpolate/interpSNd/setup.py	2008-08-19 19:12:18 UTC (rev 4654)
@@ -0,0 +1,19 @@
+"""setup"""
+
+import os, sys, setuptools
+
+def configuration(parent_package='', top_path=None):
+    from numpy.distutils.misc_util import Configuration
+    config = Configuration('', parent_package, top_path)
+    config.add_extension('_triang',
+                                    sources = ['_triang.cpp'],
+                                    depends = [],
+                                    )
+    
+    return config
+    
+if __name__ == '__main__':
+    from numpy.distutils.core import setup
+    setup(
+        **configuration().todict()
+    )

Added: branches/interpolate/interpSNd/test.py
===================================================================
--- branches/interpolate/interpSNd/test.py	2008-08-19 19:07:08 UTC (rev 4653)
+++ branches/interpolate/interpSNd/test.py	2008-08-19 19:12:18 UTC (rev 4654)
@@ -0,0 +1,53 @@
+# P for pathological
+import unittest
+import numpy as np
+from numpy import array
+import dewall as dw
+import interpolateSNd as SNd
+reload(dw)
+
+class Test(unittest.TestCase):
+    
+    def compare_array(self, a, b):
+        return np.allclose(a,b)
+    
+    ## test Delaunay triangulation itself
+    
+    # testing pathological cases with
+    
+    def test_square(self):
+        P = [array([0.,1.]), array([0.,0.]), array([1.,0.]), array([1.,1.])]
+        tri = dw.dewall(P)
+        self.assert_( len(tri)==2 )
+        self.assert_( len(tri[0])==3 )
+        self.assert_( len(tri[1])==3 )
+        print "square triangulation:\n", tri
+        
+    def test_linear(self):
+        P = [array([0.,1.]), array([0.,0.]), array([0.,-1.])]
+        tri = dw.dewall(P)
+        print "line triang:\n", tri
+        
+    # testing general case using random data
+        
+    ## test interpolation, and thus also triangulation by extension
+    def test_linear_on_cube(self):
+        x = array([0., 1, 0, 1, 0, 1, 0, 1])
+        y = array([0., 0, 1, 1, 0, 0, 1, 1])
+        z = array([0., 0, 0, 0, 1, 1, 1, 1])
+        points = array([x,y,z]).reshape((3,8))
+        fvals = x+y-z
+        
+        interp = SNd.InterpolateSNd(points, fvals)
+        
+        newdata = np.random.random_sample((3,20))
+        interpvals = interp(newdata)
+        
+        realvals = newdata[0,:]+newdata[1,:]-newdata[2,:]
+        
+        self.assert_(compare_array(np.ravel(interpvals), np.ravel(realvals)))
+        
+        
+        
+if __name__ == "__main__":
+    unittest.main()
\ No newline at end of file



More information about the Scipy-svn mailing list