Context Navigation

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

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
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?
13
14	TO DO
15	* remove points outside the mesh ?(in interpolate_block)?
16	* geo-ref (in interpolate_block)
17	* add functional interpolate interface - in mesh and points, out interp data
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
33	from search_functions import search_tree_of_vertices
34
35	class Interpolate:
36
37	def __init__(self,
38	vertex_coordinates,
39	triangles,
40	mesh_origin = None,
41	verbose=False,
42	max_vertices_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 a mesh (or 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	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.
62
63	max_vertices_per_cell: Number of vertices in a quad tree cell
64	at which the cell is split into 4.
65
66	"""
67
68	# Initialise variabels
69	self._A_can_be_reused = False
70	self._point_coordinates = None
71
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
87	#FIXME: move this to interpolate_block
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,
98	max_points_per_cell = max_vertices_per_cell)
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
108	This algorithm uses a quad tree data structure for fast binning
109	of data points
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 = \
146	search_tree_of_vertices(self.root, self.mesh, x)
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:
161	print 'Could not find triangle for point', x
162	return A
163
164	def _search_tree_of_vertices_OBSOLETE(self, root, mesh, x):
165	"""
166	Find the triangle (element) that the point x is in.
167
168	root: A quad tree of the vertices
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
174
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
182	candidate_vertices = root.search(x[0], x[1])
183	is_more_elements = True
184
185	element_found, sigma0, sigma1, sigma2, k = \
186	self._search_triangles_of_vertices(mesh,
187	candidate_vertices, x)
188	while not element_found and is_more_elements:
189	candidate_vertices, branch = root.expand_search()
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 = \
194	self._search_triangles_of_vertices(mesh,
195	candidate_vertices, x)
196	is_more_elements = False
197	else:
198	element_found, sigma0, sigma1, sigma2, k = \
199	self._search_triangles_of_vertices(mesh,
200	candidate_vertices, x)
201
202	return element_found, sigma0, sigma1, sigma2, k
203
204	def _search_triangles_of_vertices_OBSOLETE(self, mesh, candidate_vertices, x):
205	#Find triangle containing x:
206	element_found = False
207
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)
227
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)
232
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
241
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
251
252	if element_found is True:
253	#Don't look for any other triangle
254	break
255	return element_found, sigma0, sigma1, sigma2, k
256
257
258
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.
263
264	f is the data on the mesh vertices.
265
266	The mesh values representing a smooth surface are
267	assumed to be specified in f.
268
269	Inputs:
270	f: Vector or array of data at the mesh vertices.
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.
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.
278	start_blocking_len: If the # of points is more or greater than this,
279	start blocking
280
281	Output:
282	Interpolated values at inputted points (z).
283	"""
284
285	# Can I interpolate, based on previous point_coordinates?
286	if point_coordinates is None:
287	if self._A_can_be_reused is True and \
288	len(self._point_coordinates) < start_blocking_len:
289	z = self._get_point_data_z(f,
290	verbose=verbose)
291	elif self._point_coordinates is not None:
292	# if verbose, give warning
293	if verbose:
294	print 'WARNING: Recalculating A matrix, due to blocking.'
295	point_coordinates = self._point_coordinates
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:
302	self._point_coordinates = point_coordinates
303	if len(point_coordinates) < start_blocking_len or \
304	start_blocking_len == 0:
305	self._A_can_be_reused = True
306	z = self.interpolate_block(f, point_coordinates,
307	verbose=verbose)
308	else:
309	#Handle blocking
310	self._A_can_be_reused = False
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
325
326	def interpolate_block(self, f, point_coordinates = None, verbose=False):
327	"""
328	Call this if you want to control the blocking or make sure blocking
329	doesn't occur.
330
331	See interpolate for doc info.
332	"""
333	if point_coordinates is not None:
334	self._A =self._build_interpolation_matrix_A(point_coordinates,
335	verbose=verbose)
336	return self._get_point_data_z(f)
337
338	def _get_point_data_z(self, f, verbose=False):
339	return self._A * f
340
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"

Note: See TracBrowser for help on using the repository browser.

Download in other formats: