[Scipy-svn] r4654 - in branches/interpolate: . interpSNd
scipy-svn at scipy.org
scipy-svn at scipy.org
Tue Aug 19 15:12:22 EDT 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