source: inundation/fit_interpolate/interpolate.py @ 2263

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

comments

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