[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 | |
---|
| 8 | DESIGN 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 | |
---|
| 14 | TO 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 | |
---|
| 20 | import time |
---|
| 21 | |
---|
| 22 | from Numeric import zeros, array, Float, Int, dot, transpose, concatenate, \ |
---|
| 23 | ArrayType, allclose, take |
---|
| 24 | |
---|
| 25 | from pyvolution.mesh import Mesh |
---|
| 26 | from pyvolution.sparse import Sparse, Sparse_CSR |
---|
| 27 | from pyvolution.cg_solve import conjugate_gradient, VectorShapeError |
---|
| 28 | from coordinate_transforms.geo_reference import Geo_reference |
---|
| 29 | from pyvolution.quad import build_quadtree |
---|
| 30 | from utilities.numerical_tools import ensure_numeric |
---|
| 31 | from utilities.polygon import inside_polygon |
---|
| 32 | |
---|
[2190] | 33 | from search_functions import search_tree_of_vertices |
---|
[2187] | 34 | |
---|
| 35 | class 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 | #------------------------------------------------------------- |
---|
| 342 | if __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" |
---|