source: inundation/fit_interpolate/interpolate.py @ 2187

Last change on this file since 2187 was 2187, checked in by duncan, 19 years ago

refactoring interpolate
File size: 11.1 KB
Line 
1"""Least squares interpolation.
2
3   Implements a least-squares interpolation.
4
5   Ole Nielsen, Stephen Roberts, Duncan Gray, Christopher Zoppou
6   Geoscience Australia, 2004.
7
8DESIGN ISSUES
9* what variables should be global?
10- if there are no global vars functions can be moved around alot easier
11
12* What will be the public interface to this class?
13"""
14
15import time
16
17from Numeric import zeros, array, Float, Int, dot, transpose, concatenate, \
18     ArrayType, allclose, take
19
20from pyvolution.mesh import Mesh
21from pyvolution.sparse import Sparse, Sparse_CSR
22from pyvolution.cg_solve import conjugate_gradient, VectorShapeError
23from coordinate_transforms.geo_reference import Geo_reference
24
25from pyvolution.quad import build_quadtree
26from utilities.numerical_tools import ensure_numeric
27from utilities.polygon import inside_polygon
28
29
30class Interpolate:
31   
32    def __init__(self,
33                 vertex_coordinates,
34                 triangles,
35                 point_coordinates = None,
36                 mesh_origin = None,
37                 verbose=False,
38                 max_points_per_cell=30):
39
40
41        """ Build interpolation matrix mapping from
42        function values at vertices to function values at data points
43
44        Inputs:
45
46          vertex_coordinates: List of coordinate pairs [xi, eta] of
47          points constituting mesh (or a an m x 2 Numeric array)
48          Points may appear multiple times
49          (e.g. if vertices have discontinuities)
50
51          triangles: List of 3-tuples (or a Numeric array) of
52          integers representing indices of all vertices in the mesh.
53
54          point_coordinates: List of coordinate pairs [x, y] of
55          data points (or an nx2 Numeric array)
56          If point_coordinates is absent, only smoothing matrix will
57          be built
58
59          alpha: Smoothing parameter
60
61          data_origin and mesh_origin are 3-tuples consisting of
62          UTM zone, easting and northing. If specified
63          point coordinates and vertex coordinates are assumed to be
64          relative to their respective origins.
65
66        """
67       
68        from pyvolution.util import ensure_numeric
69
70        #Convert input to Numeric arrays
71        triangles = ensure_numeric(triangles, Int)
72        vertex_coordinates = ensure_numeric(vertex_coordinates, Float)
73
74        #Build underlying mesh
75        if verbose: print 'Building mesh'
76        #self.mesh = General_mesh(vertex_coordinates, triangles,
77        #FIXME: Trying the normal mesh while testing precrop,
78        #       The functionality of boundary_polygon is needed for that
79
80        #FIXME - geo ref does not have to go into mesh.
81        # Change the point co-ords to conform to the
82        # mesh co-ords early in the code
83
84        #FIXME: geo_ref can also be a geo_ref object
85        if mesh_origin is None:
86            geo = None
87        else:
88            geo = Geo_reference(mesh_origin[0],mesh_origin[1],mesh_origin[2])
89        self.mesh = Mesh(vertex_coordinates, triangles,
90                         geo_reference = geo)
91       
92        self.mesh.check_integrity()
93
94        self.root = build_quadtree(self.mesh,
95                              max_points_per_cell = max_points_per_cell)
96       
97       
98    def _build_interpolation_matrix_A(self,
99                                     point_coordinates,
100                                     verbose = False):
101        """Build n x m interpolation matrix, where
102        n is the number of data points and
103        m is the number of basis functions phi_k (one per vertex)
104
105        This algorithm uses a quad tree data structure for fast binning of data points
106        origin is a 3-tuple consisting of UTM zone, easting and northing.
107        If specified coordinates are assumed to be relative to this origin.
108
109        This one will override any data_origin that may be specified in
110        instance interpolation
111
112        Preconditions
113        Point_coordindates and mesh vertices have the same origin.
114        """
115
116
117       
118        #Convert point_coordinates to Numeric arrays, in case it was a list.
119        point_coordinates = ensure_numeric(point_coordinates, Float)
120       
121        #Remove points falling outside mesh boundary
122        # do this bit later - that sorta means this becomes an object
123        # get a list of what indices are outside the boundary
124        # maybe fill these rows with n/a?
125
126       
127        #Build n x m interpolation matrix
128        m = self.mesh.coordinates.shape[0] #Nbr of basis functions (1/vertex)
129        n = point_coordinates.shape[0]     #Nbr of data points
130
131        if verbose: print 'Number of datapoints: %d' %n
132        if verbose: print 'Number of basis functions: %d' %m
133
134        A = Sparse(n,m)
135
136
137        # I think this (along with other expanded_quad_searches stuff) can be removed
138        #self.expanded_quad_searches = []
139
140        #Compute matrix elements
141        for i in range(n):
142            #For each data_coordinate point
143
144            if verbose and i%((n+10)/10)==0: print 'Doing %d of %d' %(i, n)
145            x = point_coordinates[i]
146
147            #Find vertices near x
148            candidate_vertices = self.root.search(x[0], x[1])
149            is_more_elements = True
150
151            element_found, sigma0, sigma1, sigma2, k = \
152                self.search_triangles_of_vertices(candidate_vertices, x)
153            first_expansion = True
154            while not element_found and is_more_elements:
155                #if verbose: print 'Expanding search'
156                if first_expansion == True:
157                    #self.expanded_quad_searches.append(1)
158                    first_expansion = False
159                #else:
160                    #end = len(self.expanded_quad_searches) - 1
161                    #assert end >= 0
162                    #self.expanded_quad_searches[end] += 1
163                candidate_vertices, branch = self.root.expand_search()
164                if branch == []:
165                    # Searching all the verts from the root cell that haven't
166                    # been searched.  This is the last try
167                    element_found, sigma0, sigma1, sigma2, k = \
168                      self.search_triangles_of_vertices(candidate_vertices, x)
169                    is_more_elements = False
170                else:
171                    element_found, sigma0, sigma1, sigma2, k = \
172                      self.search_triangles_of_vertices(candidate_vertices, x)
173
174               
175            #Update interpolation matrix A if necessary
176            if element_found is True:
177                #Assign values to matrix A
178
179                j0 = self.mesh.triangles[k,0] #Global vertex id for sigma0
180                j1 = self.mesh.triangles[k,1] #Global vertex id for sigma1
181                j2 = self.mesh.triangles[k,2] #Global vertex id for sigma2
182
183                sigmas = {j0:sigma0, j1:sigma1, j2:sigma2}
184                js     = [j0,j1,j2]
185
186                for j in js:
187                    A[i,j] = sigmas[j]
188                   
189            else:
190                print 'Could not find triangle for point', x
191
192
193
194        return A
195
196
197    def search_triangles_of_vertices(self,
198                                 candidate_vertices, x):
199            #Find triangle containing x:
200            element_found = False
201
202            # This will be returned if element_found = False
203            sigma2 = -10.0
204            sigma0 = -10.0
205            sigma1 = -10.0
206            k = -10.0
207            #print "*$* candidate_vertices", candidate_vertices
208            #For all vertices in same cell as point x
209            for v in candidate_vertices:
210                #FIXME (DSG-DSG): this catches verts with no triangle.
211                #Currently pmesh is producing these.
212                #this should be stopped,
213                if self.mesh.vertexlist[v] is None:
214                    continue
215                #for each triangle id (k) which has v as a vertex
216                for k, _ in self.mesh.vertexlist[v]:
217
218                    #Get the three vertex_points of candidate triangle
219                    xi0 = self.mesh.get_vertex_coordinate(k, 0)
220                    xi1 = self.mesh.get_vertex_coordinate(k, 1)
221                    xi2 = self.mesh.get_vertex_coordinate(k, 2)
222
223                    #print "PDSG - k", k
224                    #print "PDSG - xi0", xi0
225                    #print "PDSG - xi1", xi1
226                    #print "PDSG - xi2", xi2
227                    #print "PDSG element %i verts((%f, %f),(%f, %f),(%f, %f))"\
228                    #   % (k, xi0[0], xi0[1], xi1[0], xi1[1], xi2[0], xi2[1])
229
230                    #Get the three normals
231                    n0 = self.mesh.get_normal(k, 0)
232                    n1 = self.mesh.get_normal(k, 1)
233                    n2 = self.mesh.get_normal(k, 2)
234
235
236                    #Compute interpolation
237                    sigma2 = dot((x-xi0), n2)/dot((xi2-xi0), n2)
238                    sigma0 = dot((x-xi1), n0)/dot((xi0-xi1), n0)
239                    sigma1 = dot((x-xi2), n1)/dot((xi1-xi2), n1)
240
241                    #print "PDSG - sigma0", sigma0
242                    #print "PDSG - sigma1", sigma1
243                    #print "PDSG - sigma2", sigma2
244
245                    #FIXME: Maybe move out to test or something
246                    epsilon = 1.0e-6
247                    assert abs(sigma0 + sigma1 + sigma2 - 1.0) < epsilon
248
249                    #Check that this triangle contains the data point
250
251                    #Sigmas can get negative within
252                    #machine precision on some machines (e.g nautilus)
253                    #Hence the small eps
254                    eps = 1.0e-15
255                    if sigma0 >= -eps and sigma1 >= -eps and sigma2 >= -eps:
256                        element_found = True
257                        break
258
259                if element_found is True:
260                    #Don't look for any other triangle
261                    break
262            return element_found, sigma0, sigma1, sigma2, k
263
264
265
266    #def get_A(self):
267     #       return self.A.todense()
268   
269    def interpolate(self, f):
270        """Interpolate mesh data f to points.
271
272        f is the data on the mesh vertices.
273
274
275        The mesh values representing a smooth surface are
276        assumed to be specified in f. This argument could,
277        for example have been obtained from the method self.fit()
278
279        Pre Condition:
280          self.A has been initialised
281
282        Inputs:
283          f: Vector or array of data at the mesh vertices.
284          If f is an array, interpolation will be done for each column as
285          per underlying matrix-matrix multiplication
286
287        Output:
288          Interpolated values at data points implied in self.A
289
290        """
291
292        return self.A * f
293#-------------------------------------------------------------
294if __name__ == "__main__":
295    a = [0.0, 0.0]
296    b = [0.0, 2.0]
297    c = [2.0,0.0]
298    points = [a, b, c]
299    vertices = [ [1,0,2] ]   #bac
300   
301    data = [ [2.0/3, 2.0/3] ] #Use centroid as one data point
302   
303    interp = Interpolate(points, vertices) #, data)
304    A  = interp._build_interpolation_matrix_A(data, verbose=True)
305    A = A.todense()
306    print "A",A
307    assert allclose(A, [[1./3, 1./3, 1./3]])
308    print "finished"
Note: See TracBrowser for help on using the repository browser.