source: inundation/fit_interpolate/interpolate.py @ 2191

Last change on this file since 2191 was 2191, checked in by duncan, 18 years ago

functionn name change

File size: 13.6 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
14TO DO
15* remove points outside the mesh
16* geo-ref
17"""
18
19import time
20
21from Numeric import zeros, array, Float, Int, dot, transpose, concatenate, \
22     ArrayType, allclose, take
23
24from pyvolution.mesh import Mesh
25from pyvolution.sparse import Sparse, Sparse_CSR
26from pyvolution.cg_solve import conjugate_gradient, VectorShapeError
27from coordinate_transforms.geo_reference import Geo_reference
28from pyvolution.quad import build_quadtree
29from utilities.numerical_tools import ensure_numeric
30from utilities.polygon import inside_polygon
31
32from search_functions import search_tree_of_vertices
33
34class Interpolate:
35   
36    def __init__(self,
37                 vertex_coordinates,
38                 triangles,
39                 #point_coordinates = None,
40                 mesh_origin = None,
41                 verbose=False,
42                 max_points_per_cell=30):
43
44
45        """ Build interpolation matrix mapping from
46        function values at vertices to function values at data points
47
48        Inputs:
49
50          vertex_coordinates: List of coordinate pairs [xi, eta] of
51          points constituting mesh (or a an m x 2 Numeric array)
52          Points may appear multiple times
53          (e.g. if vertices have discontinuities)
54
55          triangles: List of 3-tuples (or a Numeric array) of
56          integers representing indices of all vertices in the mesh.
57
58          point_coordinates: List of coordinate pairs [x, y] of
59          data points (or an nx2 Numeric array)
60          If point_coordinates is absent, only smoothing matrix will
61          be built
62
63          alpha: Smoothing parameter
64
65          data_origin and mesh_origin are 3-tuples consisting of
66          UTM zone, easting and northing. If specified
67          point coordinates and vertex coordinates are assumed to be
68          relative to their respective origins.
69
70        """
71       
72        from pyvolution.util import ensure_numeric
73
74        # Initialise variabels
75        self.A_can_be_reused = False
76        self.point_coordinates = None
77       
78        #Convert input to Numeric arrays
79        triangles = ensure_numeric(triangles, Int)
80        vertex_coordinates = ensure_numeric(vertex_coordinates, Float)
81
82        #Build underlying mesh
83        if verbose: print 'Building mesh'
84        #self.mesh = General_mesh(vertex_coordinates, triangles,
85        #FIXME: Trying the normal mesh while testing precrop,
86        #       The functionality of boundary_polygon is needed for that
87
88        #FIXME - geo ref does not have to go into mesh.
89        # Change the point co-ords to conform to the
90        # mesh co-ords early in the code
91
92        #FIXME: geo_ref can also be a geo_ref object
93        if mesh_origin is None:
94            geo = None
95        else:
96            geo = Geo_reference(mesh_origin[0],mesh_origin[1],mesh_origin[2])
97        self.mesh = Mesh(vertex_coordinates, triangles,
98                         geo_reference = geo)
99       
100        self.mesh.check_integrity()
101
102        self.root = build_quadtree(self.mesh,
103                              max_points_per_cell = max_points_per_cell)
104       
105       
106    def _build_interpolation_matrix_A(self,
107                                     point_coordinates,
108                                     verbose = False):
109        """Build n x m interpolation matrix, where
110        n is the number of data points and
111        m is the number of basis functions phi_k (one per vertex)
112
113        This algorithm uses a quad tree data structure for fast binning of data points
114        origin is a 3-tuple consisting of UTM zone, easting and northing.
115        If specified coordinates are assumed to be relative to this origin.
116
117        This one will override any data_origin that may be specified in
118        instance interpolation
119
120        Preconditions
121        Point_coordindates and mesh vertices have the same origin.
122        """
123
124
125       
126        #Convert point_coordinates to Numeric arrays, in case it was a list.
127        point_coordinates = ensure_numeric(point_coordinates, Float)
128       
129        #Remove points falling outside mesh boundary
130        # do this bit later - that sorta means this becomes an object
131        # get a list of what indices are outside the boundary
132        # maybe fill these rows with n/a?
133
134       
135        #Build n x m interpolation matrix
136        m = self.mesh.coordinates.shape[0] #Nbr of basis functions (1/vertex)
137        n = point_coordinates.shape[0]     #Nbr of data points
138
139        if verbose: print 'Number of datapoints: %d' %n
140        if verbose: print 'Number of basis functions: %d' %m
141
142        A = Sparse(n,m)
143
144        #Compute matrix elements
145        for i in range(n):
146            #For each data_coordinate point
147            if verbose and i%((n+10)/10)==0: print 'Doing %d of %d' %(i, n)
148
149            x = point_coordinates[i]
150           
151            element_found, sigma0, sigma1, sigma2, k = \
152                           search_tree_of_vertices(self.root, self.mesh, x)
153            #Update interpolation matrix A if necessary
154            if element_found is True:
155                #Assign values to matrix A
156
157                j0 = self.mesh.triangles[k,0] #Global vertex id for sigma0
158                j1 = self.mesh.triangles[k,1] #Global vertex id for sigma1
159                j2 = self.mesh.triangles[k,2] #Global vertex id for sigma2
160
161                sigmas = {j0:sigma0, j1:sigma1, j2:sigma2}
162                js     = [j0,j1,j2]
163
164                for j in js:
165                    A[i,j] = sigmas[j]
166            else:
167                print 'Could not find triangle for point', x
168               
169        return A
170
171    def _search_tree_of_vertices_OBSOLETE(self, root, mesh, x):
172        """
173        Find the triangle (element) that the point x is in.
174
175        root: A quad tree of the vertices
176        Return the associated sigma and k values
177           (and if the element was found) .
178        """
179        #Find triangle containing x:
180        element_found = False
181
182        # This will be returned if element_found = False
183        sigma2 = -10.0
184        sigma0 = -10.0
185        sigma1 = -10.0
186        k = -10.0
187           
188        #Find vertices near x
189        candidate_vertices = root.search(x[0], x[1])
190        is_more_elements = True
191
192        element_found, sigma0, sigma1, sigma2, k = \
193                       self._search_triangles_of_vertices(mesh,
194                                                          candidate_vertices, x)
195        while not element_found and is_more_elements:
196            candidate_vertices, branch = root.expand_search()
197            if branch == []:
198                # Searching all the verts from the root cell that haven't
199                # been searched.  This is the last try
200                element_found, sigma0, sigma1, sigma2, k = \
201                      self._search_triangles_of_vertices(mesh,
202                                                         candidate_vertices, x)
203                is_more_elements = False
204            else:
205                element_found, sigma0, sigma1, sigma2, k = \
206                      self._search_triangles_of_vertices(mesh,
207                                                         candidate_vertices, x)
208
209        return element_found, sigma0, sigma1, sigma2, k
210
211    def _search_triangles_of_vertices_OBSOLETE(self, mesh, candidate_vertices, x):
212        #Find triangle containing x:
213        element_found = False
214
215        # This will be returned if element_found = False
216        sigma2 = -10.0
217        sigma0 = -10.0
218        sigma1 = -10.0
219        k = -10.0
220        #print "*$* candidate_vertices", candidate_vertices
221        #For all vertices in same cell as point x
222        for v in candidate_vertices:
223            #FIXME (DSG-DSG): this catches verts with no triangle.
224            #Currently pmesh is producing these.
225            #this should be stopped,
226            if mesh.vertexlist[v] is None:
227                continue
228            #for each triangle id (k) which has v as a vertex
229            for k, _ in mesh.vertexlist[v]:
230                #Get the three vertex_points of candidate triangle
231                xi0 = mesh.get_vertex_coordinate(k, 0)
232                xi1 = mesh.get_vertex_coordinate(k, 1)
233                xi2 = mesh.get_vertex_coordinate(k, 2)
234
235                #Get the three normals
236                n0 = mesh.get_normal(k, 0)
237                n1 = mesh.get_normal(k, 1)
238                n2 = mesh.get_normal(k, 2)
239
240                #Compute interpolation
241                sigma2 = dot((x-xi0), n2)/dot((xi2-xi0), n2)
242                sigma0 = dot((x-xi1), n0)/dot((xi0-xi1), n0)
243                sigma1 = dot((x-xi2), n1)/dot((xi1-xi2), n1)
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     # FIXME: What is a good start_blocking_count value?
267    def interpolate(self, f, point_coordinates = None,
268                    start_blocking_len = 500000, verbose=False):
269        """Interpolate mesh data f to determine values, z, at points.
270
271        f is the data on the mesh vertices.
272
273        The mesh values representing a smooth surface are
274        assumed to be specified in f.
275
276        Inputs:
277          f: Vector or array of data at the mesh vertices.
278              If f is an array, interpolation will be done for each column as
279              per underlying matrix-matrix multiplication
280          point_coordinates: Interpolate mesh data to these positions.
281          start_blocking_len: If the # of points is more or greater than this,
282              start blocking
283
284        Output:
285          Interpolated values at inputted points (z).
286        """
287       
288        # Can I interpolate, based on previous point_coordinates?
289        if point_coordinates is None:
290            if self.A_can_be_reused is True and \
291            len(self.point_coordinates) < start_blocking_len: 
292                z = self._get_point_data_z(f,
293                                           verbose=verbose)
294            elif self.point_coordinates is not None:
295                #     if verbose, give warning
296                if verbose:
297                    print 'WARNING: Recalculating A matrix, due to blocking.'
298                point_coordinates = self.point_coordinates
299            else:
300                #There are no good point_coordinates. import sys; sys.exit()
301                msg = 'ERROR (interpolate.py): No point_coordinates inputted'
302                raise msg
303           
304        if point_coordinates is not None:
305            self.point_coordinates = point_coordinates
306            if len(point_coordinates) < start_blocking_len or \
307                   start_blocking_len == 0:
308                self.A_can_be_reused = True
309                z = self.interpolate_block(f, point_coordinates,
310                                           verbose=verbose)
311            else:
312                #Handle blocking
313                self.A_can_be_reused = False
314                start=0
315                z = self.interpolate_block(f, point_coordinates[0:0])
316                for end in range(start_blocking_len
317                                 ,len(point_coordinates)
318                                 ,start_blocking_len):
319                    t = self.interpolate_block(f, point_coordinates[start:end],
320                                               verbose=verbose)
321                    z = concatenate((z,t))
322                    start = end
323                end = len(point_coordinates)
324                t = self.interpolate_block(f, point_coordinates[start:end],
325                                           verbose=verbose)
326                z = concatenate((z,t))
327        return z
328
329    def interpolate_block(self, f, point_coordinates = None, verbose=False):
330        """
331        Call this if you want to control the blocking or make sure blocking
332        doesn't occur.
333
334        See interpolate for doc info.
335        """
336        if point_coordinates is not None:
337            self.A =self._build_interpolation_matrix_A(point_coordinates,
338                                                       verbose=verbose)
339        return self._get_point_data_z(f)
340
341    def _get_point_data_z(self, f, verbose=False):
342        return self.A * f
343       
344#-------------------------------------------------------------
345if __name__ == "__main__":
346    a = [0.0, 0.0]
347    b = [0.0, 2.0]
348    c = [2.0,0.0]
349    points = [a, b, c]
350    vertices = [ [1,0,2] ]   #bac
351   
352    data = [ [2.0/3, 2.0/3] ] #Use centroid as one data point
353   
354    interp = Interpolate(points, vertices) #, data)
355    A  = interp._build_interpolation_matrix_A(data, verbose=True)
356    A = A.todense()
357    print "A",A
358    assert allclose(A, [[1./3, 1./3, 1./3]])
359    print "finished"
Note: See TracBrowser for help on using the repository browser.