Context Navigation

source: inundation/ga/storm_surge/pyvolution/least_squares.py @ 1454

Last change on this file since 1454 was 1423, checked in by duncan, 20 years ago
I've removed the function mesh_file_to_mesh_dictionary. Please use the function import_mesh_file instead.
File size: 32.2 KB

Line
1	"""Least squares smooting and interpolation.
2
3	Implements a penalised least-squares fit and associated interpolations.
4
5	The penalty term (or smoothing term) is controlled by the smoothing
6	parameter alpha.
7	With a value of alpha=0, the fit function will attempt
8	to interpolate as closely as possible in the least-squares sense.
9	With values alpha > 0, a certain amount of smoothing will be applied.
10	A positive alpha is essential in cases where there are too few
11	data points.
12	A negative alpha is not allowed.
13	A typical value of alpha is 1.0e-6
14
15
16	Ole Nielsen, Stephen Roberts, Duncan Gray, Christopher Zoppou
17	Geoscience Australia, 2004.
18	"""
19
20	import exceptions
21	class ShapeError(exceptions.Exception): pass
22
23	#from general_mesh import General_mesh
24	from Numeric import zeros, array, Float, Int, dot, transpose, concatenate, ArrayType
25	from mesh import Mesh
26
27	from Numeric import zeros, take, array, Float, Int, dot, transpose, concatenate, ArrayType
28	from sparse import Sparse, Sparse_CSR
29	from cg_solve import conjugate_gradient, VectorShapeError
30
31	from coordinate_transforms.geo_reference import Geo_reference
32
33	import time
34
35
36	try:
37	from util import gradient
38	except ImportError, e:
39	#FIXME reduce the dependency of modules in pyvolution
40	# Have util in a dir, working like load_mesh, and get rid of this
41	def gradient(x0, y0, x1, y1, x2, y2, q0, q1, q2):
42	"""
43	"""
44
45	det = (y2-y0)(x1-x0) - (y1-y0)(x2-x0)
46	a = (y2-y0)(q1-q0) - (y1-y0)(q2-q0)
47	a /= det
48
49	b = (x1-x0)(q2-q0) - (x2-x0)(q1-q0)
50	b /= det
51
52	return a, b
53
54
55	DEFAULT_ALPHA = 0.001
56
57	def fit_to_mesh_file(mesh_file, point_file, mesh_output_file,
58	alpha=DEFAULT_ALPHA, verbose= False,
59	expand_search = False,
60	data_origin = None,
61	mesh_origin = None,
62	precrop = False):
63	"""
64	Given a mesh file (tsh) and a point attribute file (xya), fit
65	point attributes to the mesh and write a mesh file with the
66	results.
67
68
69	If data_origin is not None it is assumed to be
70	a 3-tuple with geo referenced
71	UTM coordinates (zone, easting, northing)
72
73	mesh_origin is the same but refers to the input tsh file.
74	FIXME: When the tsh format contains it own origin, these parameters can go.
75	FIXME: And both origins should be obtained from the specified files.
76	"""
77
78	from load_mesh.loadASCII import import_mesh_file, \
79	import_points_file, export_mesh_file, \
80	concatinate_attributelist
81
82	mesh_dict = import_mesh_file(mesh_file)
83	vertex_coordinates = mesh_dict['vertices']
84	triangles = mesh_dict['triangles']
85	if type(mesh_dict['vertex_attributes']) == ArrayType:
86	old_point_attributes = mesh_dict['vertex_attributes'].tolist()
87	else:
88	old_point_attributes = mesh_dict['vertex_attributes']
89
90	if type(mesh_dict['vertex_attribute_titles']) == ArrayType:
91	old_title_list = mesh_dict['vertex_attribute_titles'].tolist()
92	else:
93	old_title_list = mesh_dict['vertex_attribute_titles']
94
95	if verbose: print 'tsh file %s loaded' %mesh_file
96
97	# load in the .pts file
98	try:
99	point_dict = import_points_file(point_file,
100	delimiter = ',',
101	verbose=verbose)
102	except SyntaxError,e:
103	point_dict = import_points_file(point_file,
104	delimiter = ' ',
105	verbose=verbose)
106
107	point_coordinates = point_dict['pointlist']
108	title_list,point_attributes = concatinate_attributelist(point_dict['attributelist'])
109
110	if point_dict.has_key('geo_reference') and not point_dict['geo_reference'] is None:
111	data_origin = point_dict['geo_reference'].get_origin()
112	else:
113	data_origin = (56, 0, 0) #FIXME(DSG-DSG)
114
115	if mesh_dict.has_key('geo_reference') and not mesh_dict['geo_reference'] is None:
116	mesh_origin = mesh_dict['geo_reference'].get_origin()
117	else:
118	mesh_origin = (56, 0, 0) #FIXME(DSG-DSG)
119
120	if verbose: print "points file loaded"
121	if verbose:print "fitting to mesh"
122	f = fit_to_mesh(vertex_coordinates,
123	triangles,
124	point_coordinates,
125	point_attributes,
126	alpha = alpha,
127	verbose = verbose,
128	expand_search = expand_search,
129	data_origin = data_origin,
130	mesh_origin = mesh_origin,
131	precrop = precrop)
132	if verbose: print "finished fitting to mesh"
133
134	# convert array to list of lists
135	new_point_attributes = f.tolist()
136	#FIXME have this overwrite attributes with the same title - DSG
137	#Put the newer attributes last
138	if old_title_list <> []:
139	old_title_list.extend(title_list)
140	#FIXME can this be done a faster way? - DSG
141	for i in range(len(old_point_attributes)):
142	old_point_attributes[i].extend(new_point_attributes[i])
143	mesh_dict['vertex_attributes'] = old_point_attributes
144	mesh_dict['vertex_attribute_titles'] = old_title_list
145	else:
146	mesh_dict['vertex_attributes'] = new_point_attributes
147	mesh_dict['vertex_attribute_titles'] = title_list
148
149	#FIXME (Ole): Remember to output mesh_origin as well
150	if verbose: print "exporting to file ",mesh_output_file
151	export_mesh_file(mesh_output_file, mesh_dict)
152
153
154	def fit_to_mesh(vertex_coordinates,
155	triangles,
156	point_coordinates,
157	point_attributes,
158	alpha = DEFAULT_ALPHA,
159	verbose = False,
160	expand_search = False,
161	data_origin = None,
162	mesh_origin = None,
163	precrop = False):
164	"""
165	Fit a smooth surface to a triangulation,
166	given data points with attributes.
167
168
169	Inputs:
170
171	vertex_coordinates: List of coordinate pairs [xi, eta] of points
172	constituting mesh (or a an m x 2 Numeric array)
173
174	triangles: List of 3-tuples (or a Numeric array) of
175	integers representing indices of all vertices in the mesh.
176
177	point_coordinates: List of coordinate pairs [x, y] of data points
178	(or an nx2 Numeric array)
179
180	alpha: Smoothing parameter.
181
182	point_attributes: Vector or array of data at the point_coordinates.
183
184	data_origin and mesh_origin are 3-tuples consisting of
185	UTM zone, easting and northing. If specified
186	point coordinates and vertex coordinates are assumed to be
187	relative to their respective origins.
188
189	"""
190	interp = Interpolation(vertex_coordinates,
191	triangles,
192	point_coordinates,
193	alpha = alpha,
194	verbose = verbose,
195	expand_search = expand_search,
196	data_origin = data_origin,
197	mesh_origin = mesh_origin,
198	precrop = precrop)
199
200	vertex_attributes = interp.fit_points(point_attributes, verbose = verbose)
201	return vertex_attributes
202
203
204
205	def pts2rectangular(pts_name, M, N, alpha = DEFAULT_ALPHA,
206	verbose = False, reduction = 1, format = 'netcdf'):
207	"""Fits attributes from pts file to MxN rectangular mesh
208
209	Read pts file and create rectangular mesh of resolution MxN such that
210	it covers all points specified in pts file.
211
212	FIXME: This may be a temporary function until we decide on
213	netcdf formats etc
214
215	FIXME: Uses elevation hardwired
216	"""
217
218	import util, mesh_factory
219
220	if verbose: print 'Read pts'
221	points, attributes = util.read_xya(pts_name, format)
222
223	#Reduce number of points a bit
224	points = points[::reduction]
225	elevation = attributes['elevation'] #Must be elevation
226	elevation = elevation[::reduction]
227
228	if verbose: print 'Got %d data points' %len(points)
229
230	if verbose: print 'Create mesh'
231	#Find extent
232	max_x = min_x = points[0][0]
233	max_y = min_y = points[0][1]
234	for point in points[1:]:
235	x = point[0]
236	if x > max_x: max_x = x
237	if x < min_x: min_x = x
238	y = point[1]
239	if y > max_y: max_y = y
240	if y < min_y: min_y = y
241
242	#Create appropriate mesh
243	vertex_coordinates, triangles, boundary =\
244	mesh_factory.rectangular(M, N, max_x-min_x, max_y-min_y,
245	(min_x, min_y))
246
247	#Fit attributes to mesh
248	vertex_attributes = fit_to_mesh(vertex_coordinates,
249	triangles,
250	points,
251	elevation, alpha=alpha, verbose=verbose)
252
253
254
255	return vertex_coordinates, triangles, boundary, vertex_attributes
256
257
258
259	class Interpolation:
260
261	def __init__(self,
262	vertex_coordinates,
263	triangles,
264	point_coordinates = None,
265	alpha = DEFAULT_ALPHA,
266	verbose = False,
267	expand_search = True,
268	max_points_per_cell = 30,
269	mesh_origin = None,
270	data_origin = None,
271	precrop = False):
272
273
274	""" Build interpolation matrix mapping from
275	function values at vertices to function values at data points
276
277	Inputs:
278
279	vertex_coordinates: List of coordinate pairs [xi, eta] of
280	points constituting mesh (or a an m x 2 Numeric array)
281
282	triangles: List of 3-tuples (or a Numeric array) of
283	integers representing indices of all vertices in the mesh.
284
285	point_coordinates: List of coordinate pairs [x, y] of
286	data points (or an nx2 Numeric array)
287	If point_coordinates is absent, only smoothing matrix will
288	be built
289
290	alpha: Smoothing parameter
291
292	data_origin and mesh_origin are 3-tuples consisting of
293	UTM zone, easting and northing. If specified
294	point coordinates and vertex coordinates are assumed to be
295	relative to their respective origins.
296
297	"""
298	from util import ensure_numeric
299
300	#Convert input to Numeric arrays
301	triangles = ensure_numeric(triangles, Int)
302	vertex_coordinates = ensure_numeric(vertex_coordinates, Float)
303
304	#Build underlying mesh
305	if verbose: print 'Building mesh'
306	#self.mesh = General_mesh(vertex_coordinates, triangles,
307	#FIXME: Trying the normal mesh while testing precrop,
308	# The functionality of boundary_polygon is needed for that
309
310	#FIXME - geo ref does not have to go into mesh.
311	# Change the point co-ords to conform to the
312	# mesh co-ords early in the code
313	if mesh_origin == None:
314	geo = None
315	else:
316	geo = Geo_reference(mesh_origin[0],mesh_origin[1],mesh_origin[2])
317	self.mesh = Mesh(vertex_coordinates, triangles,
318	geo_reference = geo)
319	#FIXME, remove if mesh checks it.
320	self.mesh.check_integrity()
321	self.data_origin = data_origin
322
323	self.point_indices = None
324
325	#Smoothing parameter
326	self.alpha = alpha
327
328	#Build coefficient matrices
329	self.build_coefficient_matrix_B(point_coordinates,
330	verbose = verbose,
331	expand_search = expand_search,
332	max_points_per_cell =\
333	max_points_per_cell,
334	data_origin = data_origin,
335	precrop = precrop)
336
337
338	def set_point_coordinates(self, point_coordinates,
339	data_origin = None):
340	"""
341	A public interface to setting the point co-ordinates.
342	"""
343	self.build_coefficient_matrix_B(point_coordinates, data_origin)
344
345	def build_coefficient_matrix_B(self, point_coordinates=None,
346	verbose = False, expand_search = True,
347	max_points_per_cell=30,
348	data_origin = None,
349	precrop = False):
350	"""Build final coefficient matrix"""
351
352
353	if self.alpha <> 0:
354	if verbose: print 'Building smoothing matrix'
355	self.build_smoothing_matrix_D()
356
357	if point_coordinates is not None:
358
359	if verbose: print 'Building interpolation matrix'
360	self.build_interpolation_matrix_A(point_coordinates,
361	verbose = verbose,
362	expand_search = expand_search,
363	max_points_per_cell =\
364	max_points_per_cell,
365	data_origin = data_origin,
366	precrop = precrop)
367
368	if self.alpha <> 0:
369	self.B = self.AtA + self.alpha*self.D
370	else:
371	self.B = self.AtA
372
373	#Convert self.B matrix to CSR format for faster matrix vector
374	self.B = Sparse_CSR(self.B)
375
376	def build_interpolation_matrix_A(self, point_coordinates,
377	verbose = False, expand_search = True,
378	max_points_per_cell=30,
379	data_origin = None,
380	precrop = False):
381	"""Build n x m interpolation matrix, where
382	n is the number of data points and
383	m is the number of basis functions phi_k (one per vertex)
384
385	This algorithm uses a quad tree data structure for fast binning of data points
386	origin is a 3-tuple consisting of UTM zone, easting and northing.
387	If specified coordinates are assumed to be relative to this origin.
388
389	This one will override any data_origin that may be specified in
390	interpolation instance
391
392	"""
393
394	from quad import build_quadtree
395	from util import ensure_numeric
396
397	if data_origin is None:
398	data_origin = self.data_origin #Use the one from
399	#interpolation instance
400
401	#Convert input to Numeric arrays just in case.
402	point_coordinates = ensure_numeric(point_coordinates, Float)
403
404
405	#Shift data points to same origin as mesh (if specified)
406
407	#FIXME this will shift if there was no geo_ref.
408	#But all this should be removed amyhow.
409	#change coords before this point
410	mesh_origin = self.mesh.geo_reference.get_origin()
411	if point_coordinates is not None:
412	if data_origin is not None:
413	if mesh_origin is not None:
414
415	#Transformation:
416	#
417	#Let x_0 be the reference point of the point coordinates
418	#and xi_0 the reference point of the mesh.
419	#
420	#A point coordinate (x + x_0) is then made relative
421	#to xi_0 by
422	#
423	# x_new = x + x_0 - xi_0
424	#
425	#and similarly for eta
426
427	x_offset = data_origin[1] - mesh_origin[1]
428	y_offset = data_origin[2] - mesh_origin[2]
429	else: #Shift back to a zero origin
430	x_offset = data_origin[1]
431	y_offset = data_origin[2]
432
433	point_coordinates[:,0] += x_offset
434	point_coordinates[:,1] += y_offset
435	else:
436	if mesh_origin is not None:
437	#Use mesh origin for data points
438	point_coordinates[:,0] -= mesh_origin[1]
439	point_coordinates[:,1] -= mesh_origin[2]
440
441
442
443	#Remove points falling outside mesh boundary
444	#This reduced one example from 1356 seconds to 825 seconds
445	#And more could be had by writing util.inside_polygon in C
446	if precrop is True:
447	from Numeric import take
448	from util import inside_polygon
449
450	if verbose: print 'Getting boundary polygon'
451	P = self.mesh.get_boundary_polygon()
452
453	if verbose: print 'Getting indices inside mesh boundary'
454	indices = inside_polygon(point_coordinates, P, verbose = verbose)
455
456	if verbose:
457	print 'Done'
458	if len(indices) != point_coordinates.shape[0]:
459	print '%d points outside mesh have been cropped.'\
460	%(point_coordinates.shape[0] - len(indices))
461	point_coordinates = take(point_coordinates, indices)
462	self.point_indices = indices
463
464
465
466
467	#Build n x m interpolation matrix
468	m = self.mesh.coordinates.shape[0] #Nbr of basis functions (1/vertex)
469	n = point_coordinates.shape[0] #Nbr of data points
470
471	if verbose: print 'Number of datapoints: %d' %n
472	if verbose: print 'Number of basis functions: %d' %m
473
474	#FIXME (Ole): We should use CSR here since mat-mat mult is now OK.
475	#However, Sparse_CSR does not have the same methods as Sparse yet
476	#The tests will reveal what needs to be done
477	self.A = Sparse(n,m)
478	self.AtA = Sparse(m,m)
479
480	#Build quad tree of vertices (FIXME: Is this the right spot for that?)
481	root = build_quadtree(self.mesh,
482	max_points_per_cell = max_points_per_cell)
483
484	#Compute matrix elements
485	for i in range(n):
486	#For each data_coordinate point
487
488	if verbose and i%((n+10)/10)==0: print 'Doing %d of %d' %(i, n)
489
490	x = point_coordinates[i]
491
492	#Find vertices near x
493	candidate_vertices = root.search(x[0], x[1])
494	is_more_elements = True
495
496	element_found, sigma0, sigma1, sigma2, k = \
497	self.search_triangles_of_vertices(candidate_vertices, x)
498	while not element_found and is_more_elements and expand_search:
499	#if verbose: print 'Expanding search'
500	candidate_vertices, branch = root.expand_search()
501	if branch == []:
502	# Searching all the verts from the root cell that haven't
503	# been searched. This is the last try
504	element_found, sigma0, sigma1, sigma2, k = \
505	self.search_triangles_of_vertices(candidate_vertices, x)
506	is_more_elements = False
507	else:
508	element_found, sigma0, sigma1, sigma2, k = \
509	self.search_triangles_of_vertices(candidate_vertices, x)
510
511
512	#Update interpolation matrix A if necessary
513	if element_found is True:
514	#Assign values to matrix A
515
516	j0 = self.mesh.triangles[k,0] #Global vertex id for sigma0
517	j1 = self.mesh.triangles[k,1] #Global vertex id for sigma1
518	j2 = self.mesh.triangles[k,2] #Global vertex id for sigma2
519
520	sigmas = {j0:sigma0, j1:sigma1, j2:sigma2}
521	js = [j0,j1,j2]
522
523	for j in js:
524	self.A[i,j] = sigmas[j]
525	for k in js:
526	self.AtA[j,k] += sigmas[j]*sigmas[k]
527	else:
528	pass
529	#Ok if there is no triangle for datapoint
530	#(as in brute force version)
531	#raise 'Could not find triangle for point', x
532
533
534
535	def search_triangles_of_vertices(self, candidate_vertices, x):
536	#Find triangle containing x:
537	element_found = False
538
539	# This will be returned if element_found = False
540	sigma2 = -10.0
541	sigma0 = -10.0
542	sigma1 = -10.0
543	k = -10.0
544
545	#For all vertices in same cell as point x
546	for v in candidate_vertices:
547
548	#for each triangle id (k) which has v as a vertex
549	for k, _ in self.mesh.vertexlist[v]:
550
551	#Get the three vertex_points of candidate triangle
552	xi0 = self.mesh.get_vertex_coordinate(k, 0)
553	xi1 = self.mesh.get_vertex_coordinate(k, 1)
554	xi2 = self.mesh.get_vertex_coordinate(k, 2)
555
556	#print "PDSG - k", k
557	#print "PDSG - xi0", xi0
558	#print "PDSG - xi1", xi1
559	#print "PDSG - xi2", xi2
560	#print "PDSG element %i verts((%f, %f),(%f, %f),(%f, %f))"\
561	# % (k, xi0[0], xi0[1], xi1[0], xi1[1], xi2[0], xi2[1])
562
563	#Get the three normals
564	n0 = self.mesh.get_normal(k, 0)
565	n1 = self.mesh.get_normal(k, 1)
566	n2 = self.mesh.get_normal(k, 2)
567
568
569	#Compute interpolation
570	sigma2 = dot((x-xi0), n2)/dot((xi2-xi0), n2)
571	sigma0 = dot((x-xi1), n0)/dot((xi0-xi1), n0)
572	sigma1 = dot((x-xi2), n1)/dot((xi1-xi2), n1)
573
574	#print "PDSG - sigma0", sigma0
575	#print "PDSG - sigma1", sigma1
576	#print "PDSG - sigma2", sigma2
577
578	#FIXME: Maybe move out to test or something
579	epsilon = 1.0e-6
580	assert abs(sigma0 + sigma1 + sigma2 - 1.0) < epsilon
581
582	#Check that this triangle contains the data point
583
584	#Sigmas can get negative within
585	#machine precision on some machines (e.g nautilus)
586	#Hence the small eps
587	eps = 1.0e-15
588	if sigma0 >= -eps and sigma1 >= -eps and sigma2 >= -eps:
589	element_found = True
590	break
591
592	if element_found is True:
593	#Don't look for any other triangle
594	break
595	return element_found, sigma0, sigma1, sigma2, k
596
597
598
599	def build_interpolation_matrix_A_brute(self, point_coordinates):
600	"""Build n x m interpolation matrix, where
601	n is the number of data points and
602	m is the number of basis functions phi_k (one per vertex)
603
604	This is the brute force which is too slow for large problems,
605	but could be used for testing
606	"""
607
608	from util import ensure_numeric
609
610	#Convert input to Numeric arrays
611	point_coordinates = ensure_numeric(point_coordinates, Float)
612
613	#Build n x m interpolation matrix
614	m = self.mesh.coordinates.shape[0] #Nbr of basis functions (1/vertex)
615	n = point_coordinates.shape[0] #Nbr of data points
616
617	self.A = Sparse(n,m)
618	self.AtA = Sparse(m,m)
619
620	#Compute matrix elements
621	for i in range(n):
622	#For each data_coordinate point
623
624	x = point_coordinates[i]
625	element_found = False
626	k = 0
627	while not element_found and k < len(self.mesh):
628	#For each triangle (brute force)
629	#FIXME: Real algorithm should only visit relevant triangles
630
631	#Get the three vertex_points
632	xi0 = self.mesh.get_vertex_coordinate(k, 0)
633	xi1 = self.mesh.get_vertex_coordinate(k, 1)
634	xi2 = self.mesh.get_vertex_coordinate(k, 2)
635
636	#Get the three normals
637	n0 = self.mesh.get_normal(k, 0)
638	n1 = self.mesh.get_normal(k, 1)
639	n2 = self.mesh.get_normal(k, 2)
640
641	#Compute interpolation
642	sigma2 = dot((x-xi0), n2)/dot((xi2-xi0), n2)
643	sigma0 = dot((x-xi1), n0)/dot((xi0-xi1), n0)
644	sigma1 = dot((x-xi2), n1)/dot((xi1-xi2), n1)
645
646	#FIXME: Maybe move out to test or something
647	epsilon = 1.0e-6
648	assert abs(sigma0 + sigma1 + sigma2 - 1.0) < epsilon
649
650	#Check that this triangle contains data point
651	if sigma0 >= 0 and sigma1 >= 0 and sigma2 >= 0:
652	element_found = True
653	#Assign values to matrix A
654
655	j0 = self.mesh.triangles[k,0] #Global vertex id
656	#self.A[i, j0] = sigma0
657
658	j1 = self.mesh.triangles[k,1] #Global vertex id
659	#self.A[i, j1] = sigma1
660
661	j2 = self.mesh.triangles[k,2] #Global vertex id
662	#self.A[i, j2] = sigma2
663
664	sigmas = {j0:sigma0, j1:sigma1, j2:sigma2}
665	js = [j0,j1,j2]
666
667	for j in js:
668	self.A[i,j] = sigmas[j]
669	for k in js:
670	self.AtA[j,k] += sigmas[j]*sigmas[k]
671	k = k+1
672
673
674
675	def get_A(self):
676	return self.A.todense()
677
678	def get_B(self):
679	return self.B.todense()
680
681	def get_D(self):
682	return self.D.todense()
683
684	#FIXME: Remember to re-introduce the 1/n factor in the
685	#interpolation term
686
687	def build_smoothing_matrix_D(self):
688	"""Build m x m smoothing matrix, where
689	m is the number of basis functions phi_k (one per vertex)
690
691	The smoothing matrix is defined as
692
693	D = D1 + D2
694
695	where
696
697	[D1]_{k,l} = \int_\Omega
698	\frac{\partial \phi_k}{\partial x}
699	\frac{\partial \phi_l}{\partial x}\,
700	dx dy
701
702	[D2]_{k,l} = \int_\Omega
703	\frac{\partial \phi_k}{\partial y}
704	\frac{\partial \phi_l}{\partial y}\,
705	dx dy
706
707
708	The derivatives \frac{\partial \phi_k}{\partial x},
709	\frac{\partial \phi_k}{\partial x} for a particular triangle
710	are obtained by computing the gradient a_k, b_k for basis function k
711	"""
712
713	#FIXME: algorithm might be optimised by computing local 9x9
714	#"element stiffness matrices:
715
716	m = self.mesh.coordinates.shape[0] #Nbr of basis functions (1/vertex)
717
718	self.D = Sparse(m,m)
719
720	#For each triangle compute contributions to D = D1+D2
721	for i in range(len(self.mesh)):
722
723	#Get area
724	area = self.mesh.areas[i]
725
726	#Get global vertex indices
727	v0 = self.mesh.triangles[i,0]
728	v1 = self.mesh.triangles[i,1]
729	v2 = self.mesh.triangles[i,2]
730
731	#Get the three vertex_points
732	xi0 = self.mesh.get_vertex_coordinate(i, 0)
733	xi1 = self.mesh.get_vertex_coordinate(i, 1)
734	xi2 = self.mesh.get_vertex_coordinate(i, 2)
735
736	#Compute gradients for each vertex
737	a0, b0 = gradient(xi0[0], xi0[1], xi1[0], xi1[1], xi2[0], xi2[1],
738	1, 0, 0)
739
740	a1, b1 = gradient(xi0[0], xi0[1], xi1[0], xi1[1], xi2[0], xi2[1],
741	0, 1, 0)
742
743	a2, b2 = gradient(xi0[0], xi0[1], xi1[0], xi1[1], xi2[0], xi2[1],
744	0, 0, 1)
745
746	#Compute diagonal contributions
747	self.D[v0,v0] += (a0a0 + b0b0)*area
748	self.D[v1,v1] += (a1a1 + b1b1)*area
749	self.D[v2,v2] += (a2a2 + b2b2)*area
750
751	#Compute contributions for basis functions sharing edges
752	e01 = (a0a1 + b0b1)*area
753	self.D[v0,v1] += e01
754	self.D[v1,v0] += e01
755
756	e12 = (a1a2 + b1b2)*area
757	self.D[v1,v2] += e12
758	self.D[v2,v1] += e12
759
760	e20 = (a2a0 + b2b0)*area
761	self.D[v2,v0] += e20
762	self.D[v0,v2] += e20
763
764
765	def fit(self, z):
766	"""Fit a smooth surface to given 1d array of data points z.
767
768	The smooth surface is computed at each vertex in the underlying
769	mesh using the formula given in the module doc string.
770
771	Pre Condition:
772	self.A, self.At and self.B have been initialised
773
774	Inputs:
775	z: Single 1d vector or array of data at the point_coordinates.
776	"""
777
778	#Convert input to Numeric arrays
779	from util import ensure_numeric
780	z = ensure_numeric(z, Float)
781
782	if len(z.shape) > 1 :
783	raise VectorShapeError, 'Can only deal with 1d data vector'
784
785	if self.point_indices is not None:
786	#Remove values for any points that were outside mesh
787	z = take(z, self.point_indices)
788
789	#Compute right hand side based on data
790	Atz = self.A.trans_mult(z)
791
792
793	#Check sanity
794	n, m = self.A.shape
795	if n<m and self.alpha == 0.0:
796	msg = 'ERROR (least_squares): Too few data points\n'
797	msg += 'There are only %d data points and alpha == 0. ' %n
798	msg += 'Need at least %d\n' %m
799	msg += 'Alternatively, set smoothing parameter alpha to a small '
800	msg += 'positive value,\ne.g. 1.0e-3.'
801	raise msg
802
803
804
805	return conjugate_gradient(self.B, Atz, Atz,imax=2*len(Atz) )
806	#FIXME: Should we store the result here for later use? (ON)
807
808
809	def fit_points(self, z, verbose=False):
810	"""Like fit, but more robust when each point has two or more attributes
811	FIXME (Ole): The name fit_points doesn't carry any meaning
812	for me. How about something like fit_multiple or fit_columns?
813	"""
814
815	try:
816	if verbose: print 'Solving penalised least_squares problem'
817	return self.fit(z)
818	except VectorShapeError, e:
819	# broadcasting is not supported.
820
821	#Convert input to Numeric arrays
822	from util import ensure_numeric
823	z = ensure_numeric(z, Float)
824
825	#Build n x m interpolation matrix
826	m = self.mesh.coordinates.shape[0] #Number of vertices
827	n = z.shape[1] #Number of data points
828
829	f = zeros((m,n), Float) #Resulting columns
830
831	for i in range(z.shape[1]):
832	f[:,i] = self.fit(z[:,i])
833
834	return f
835
836
837	def interpolate(self, f):
838	"""Evaluate smooth surface f at data points implied in self.A.
839
840	The mesh values representing a smooth surface are
841	assumed to be specified in f. This argument could,
842	for example have been obtained from the method self.fit()
843
844	Pre Condition:
845	self.A has been initialised
846
847	Inputs:
848	f: Vector or array of data at the mesh vertices.
849	If f is an array, interpolation will be done for each column as
850	per underlying matrix-matrix multiplication
851
852	Output:
853	Interpolated values at data points implied in self.A
854
855	"""
856
857	return self.A * f
858
859	def cull_outsiders(self, f):
860	pass
861
862
863	#-------------------------------------------------------------
864	if __name__ == "__main__":
865	"""
866	Load in a mesh and data points with attributes.
867	Fit the attributes to the mesh.
868	Save a new mesh file.
869	"""
870	import os, sys
871	usage = "usage: %s mesh_input.tsh point.xya mesh_output.tsh [expand\|no_expand][vervose\|non_verbose] [alpha]"\
872	%os.path.basename(sys.argv[0])
873
874	if len(sys.argv) < 4:
875	print usage
876	else:
877	mesh_file = sys.argv[1]
878	point_file = sys.argv[2]
879	mesh_output_file = sys.argv[3]
880
881	expand_search = False
882	if len(sys.argv) > 4:
883	if sys.argv[4][0] == "e" or sys.argv[4][0] == "E":
884	expand_search = True
885	else:
886	expand_search = False
887
888	verbose = False
889	if len(sys.argv) > 5:
890	if sys.argv[5][0] == "n" or sys.argv[5][0] == "N":
891	verbose = False
892	else:
893	verbose = True
894
895	if len(sys.argv) > 6:
896	alpha = sys.argv[6]
897	else:
898	alpha = DEFAULT_ALPHA
899
900	t0 = time.time()
901	fit_to_mesh_file(mesh_file,
902	point_file,
903	mesh_output_file,
904	alpha,
905	verbose= verbose,
906	expand_search = expand_search)
907
908	print 'That took %.2f seconds' %(time.time()-t0)
909

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

Download in other formats: