Context Navigation

source: anuga_core/source/anuga/pmesh/mesh.py @ 4668

Last change on this file since 4668 was 4668, checked in by duncan, 17 years ago
ticket#189 - hacky fix to unproven memory leak
File size: 129.7 KB

Line
1	#!/usr/bin/env python
2	#
3	"""General 2D triangular classes for triangular mesh generation.
4
5	Note: A .index attribute is added to objects such as vertices and
6	segments, often when reading and writing to files. This information
7	should not be used as persistant information. It is not the 'index' of
8	an element in a list.
9
10
11	Copyright 2003/2004
12	Ole Nielsen, Stephen Roberts, Duncan Gray, Christopher Zoppou
13	Geoscience Australia
14	"""
15
16	import sys
17	import math
18	import re
19	import os
20	import pickle
21
22	import types
23	import exceptions
24	from Numeric import array, Float, Int
25
26
27	class NoTrianglesError(exceptions.Exception): pass
28
29	#import load_mesh
30	from anuga.coordinate_transforms.geo_reference import Geo_reference, \
31	DEFAULT_ZONE
32	from anuga.utilities.polygon import point_in_polygon
33	from anuga.load_mesh.loadASCII import NOMAXAREA, export_mesh_file, \
34	import_mesh_file
35	import anuga.alpha_shape.alpha_shape
36	from anuga.geospatial_data.geospatial_data import Geospatial_data, \
37	ensure_geospatial, ensure_absolute, ensure_numeric
38	from anuga.mesh_engine.mesh_engine import generate_mesh
39
40	try:
41	import kinds
42	except ImportError:
43	# Hand-built mockup of the things we need from the kinds package, since it
44	# was recently removed from the standard Numeric distro. Some users may
45	# not have it by default.
46	class _bunch:
47	pass
48
49	class _kinds(_bunch):
50	default_float_kind = _bunch()
51	default_float_kind.MIN = 2.2250738585072014e-308 #smallest +ve number
52	default_float_kind.MAX = 1.7976931348623157e+308
53
54	kinds = _kinds()
55
56	SET_COLOUR='red'
57
58	#FIXME: this is not tested.
59	from anuga.utilities.numerical_tools import gradient
60
61
62
63	# 1st and third values must be the same
64	# FIXME: maybe make this a switch that the user can change? - DSG
65	initialconversions = ['', 'exterior', '']
66
67	#from os import sep
68	#sys.path.append('..'+sep+'pmesh')
69	#print "sys.path",sys.path
70
71	class MeshObject:
72	"""
73	An abstract superclass for the basic mesh objects, eg vertex, segment,
74	triangle.
75	"""
76	def __init__(self):
77	pass
78
79	class Point(MeshObject):
80	"""
81	Define a point in a 2D space.
82	"""
83	def __init__(self,X,Y):
84	__slots__ = ['x','y']
85	self.x=X
86	self.y=Y
87
88	def DistanceToPoint(self, OtherPoint):
89	"""
90	Returns the distance from this point to another
91	"""
92	SumOfSquares = ((self.x - OtherPoint.x)**2) + \
93	((self.y - OtherPoint.y)**2)
94	return math.sqrt(SumOfSquares)
95
96	def IsInsideCircle(self, Center, Radius):
97	"""
98	Return 1 if this point is inside the circle,
99	0 otherwise
100	"""
101
102	if (self.DistanceToPoint(Center)<Radius):
103	return 1
104	else:
105	return 0
106
107	def __repr__(self):
108	return "(%f,%f)" % (self.x,self.y)
109
110	def cmp_xy(self, point):
111	if self.x < point.x:
112	return -1
113	elif self.x > point.x:
114	return 1
115	else:
116	if self.y < point.y:
117	return -1
118	elif self.y > point.y:
119	return 1
120	else:
121	return 0
122
123	def same_x_y(self, point):
124	if self.x == point.x and self.y == point.y:
125	return True
126	else:
127	return False
128
129
130
131	class Vertex(Point):
132	"""
133	A point on the mesh.
134	Object attributes based on the Triangle program
135	"""
136	def __init__(self,X,Y, attributes = None):
137	__slots__ = ['x','y','attributes']
138
139	assert (type(X) == types.FloatType or type(X) == types.IntType)
140	assert (type(Y) == types.FloatType or type(Y) == types.IntType)
141	self.x=X
142	self.y=Y
143	self.attributes=[]
144
145	if attributes is None:
146	self.attributes=[]
147	else:
148	self.attributes=attributes
149
150
151	def setAttributes(self,attributes):
152	"""
153	attributes is a list.
154	"""
155	self.attributes = attributes
156
157	VERTEXSQUARESIDELENGTH = 6
158	def draw(self, canvas, tags, colour = 'black',scale = 1, xoffset = 0,
159	yoffset =0, ):
160	x = scale*(self.x + xoffset)
161	y = -1scale(self.y + yoffset) # - since for a canvas - is up
162	#print "draw x:", x
163	#print "draw y:", y
164	cornerOffset= self.VERTEXSQUARESIDELENGTH/2
165
166	# A hack to see the vert tags
167	# note: there will be many tags, since tags will not be removed
168	#when zooming
169	#canvas.create_text(x+ 2*cornerOffset,
170	# y+ 2*cornerOffset,
171	# text=tags)
172
173	return canvas.create_rectangle(x-cornerOffset,
174	y-cornerOffset,
175	x+cornerOffset,
176	y+cornerOffset,
177	tags = tags,
178	outline=colour,
179	fill = 'white')
180
181	#return tags
182
183	def __repr__(self):
184	return "[(%f,%f),%r]" % (self.x,self.y,self.attributes)
185
186	class Hole(Point):
187	"""
188	A region of the mesh were no triangles are generated.
189	Defined by a point in the hole enclosed by segments.
190	"""
191
192	HOLECORNERLENGTH = 6
193
194	def draw(self, canvas, tags, colour = 'purple',scale = 1, xoffset = 0,
195	yoffset =0, ):
196	x = scale*(self.x + xoffset)
197	y = -1scale(self.y + yoffset) # - since for a canvas - is up
198	#print "draw x:", x
199	#print "draw y:", y
200	cornerOffset= self.HOLECORNERLENGTH/2
201	return canvas.create_oval(x-cornerOffset,
202	y-cornerOffset,
203	x+cornerOffset,
204	y+cornerOffset,
205	tags = tags,
206	outline=colour,
207	fill = 'white')
208
209	class Region(Point):
210	"""
211	A region of the mesh, defined by a point in the region
212	enclosed by segments. Used to tag areas.
213	"""
214	CROSSLENGTH = 6
215	TAG = 0
216	MAXAREA = 1
217
218	def __init__(self,X,Y, tag = None, maxArea = None):
219	"""Precondition: tag is a string and maxArea is a real
220	"""
221	# This didn't work.
222	#super(Region,self)._init_(self,X,Y)
223	self.x=X
224	self.y=Y
225	self.attributes =[] # index 0 is the tag string
226	#optoinal index 1 is the max triangle area
227	#NOTE the size of this attribute is assumed
228	# to be 1 or 2 in regionstrings2int
229	if tag is None:
230	self.attributes.append("")
231	else:
232	self.attributes.append(tag) #this is a string
233
234	if maxArea is not None:
235	self.setMaxArea(maxArea) # maxArea is a number
236
237	def getTag(self,):
238	return self.attributes[self.TAG]
239
240	def setTag(self,tag):
241	self.attributes[self.TAG] = tag
242
243	def getMaxArea(self):
244	""" Returns the Max Area of a Triangle or
245	None, if the Max Area has not been set.
246	"""
247	if self.isMaxArea():
248	return self.attributes[1]
249	else:
250	return None
251
252	def setMaxArea(self,MaxArea):
253	if MaxArea is not None:
254	if self.isMaxArea():
255	self.attributes[self.MAXAREA] = float(MaxArea)
256	else:
257	self.attributes.append( float(MaxArea) )
258
259	def deleteMaxArea(self):
260	if self.isMaxArea():
261	self.attributes.pop(self.MAXAREA)
262
263	def isMaxArea(self):
264	return len(self.attributes)> 1
265
266	def draw(self, canvas, tags, scale=1, xoffset = 0, yoffset =0,
267	colour = "black"):
268	"""
269	Draw a black cross, returning the objectID
270	"""
271	x = scale*(self.x + xoffset)
272	y = -1scale(self.y + yoffset)
273	cornerOffset= self.CROSSLENGTH/2
274	return canvas.create_polygon(x,
275	y-cornerOffset,
276	x,
277	y,
278	x+cornerOffset,
279	y,
280	x,
281	y,
282	x,
283	y+cornerOffset,
284	x,
285	y,
286	x-cornerOffset,
287	y,
288	x,
289	y,
290	tags = tags,
291	outline = colour,fill = '')
292
293	def __repr__(self):
294	if self.isMaxArea():
295	area = self.getMaxArea()
296	return "(%f,%f,%s,%f)" % (self.x,self.y,
297	self.getTag(), self.getMaxArea())
298	else:
299	return "(%f,%f,%s)" % (self.x,self.y,
300	self.getTag())
301
302	class Triangle(MeshObject):
303	"""
304	A triangle element, defined by 3 vertices.
305	Attributes based on the Triangle program.
306	"""
307
308	def __init__(self, vertex1, vertex2, vertex3, attribute = None,
309	neighbors = None ):
310	"""
311	Vertices, the initial arguments, are listed in counterclockwise order.
312	"""
313	self.vertices= [vertex1,vertex2, vertex3 ]
314
315	if attribute is None:
316	self.attribute =""
317	else:
318	self.attribute = attribute #this is a string
319
320	if neighbors is None:
321	self.neighbors=[]
322	else:
323	self.neighbors=neighbors
324
325	def replace(self,new_triangle):
326	self = new_triangle
327
328	def longestSideID(self):
329	ax = self.vertices[0].x
330	ay = self.vertices[0].y
331
332	bx = self.vertices[1].x
333	by = self.vertices[1].y
334
335	cx = self.vertices[2].x
336	cy = self.vertices[2].y
337
338	lenA = ((cx-bx)2+(cy-by)2)**0.5
339	lenB = ((ax-cx)2+(ay-cy)2)**0.5
340	lenC = ((bx-ax)2+(by-ay)2)**0.5
341
342	len = [lenA,lenB,lenC]
343	return len.index(max(len))
344
345	def rotate(self,offset):
346	"""
347	permute the order of the sides of the triangle
348	offset must be 0,1 or 2
349	"""
350
351	if offset == 0:
352	pass
353	else:
354	if offset == 1:
355	self.vertices = [self.vertices[1],self.vertices[2],
356	self.vertices[0]]
357	self.neighbors = [self.neighbors[1],self.neighbors[2],
358	self.neighbors[0]]
359	if offset == 2:
360	self.vertices = [self.vertices[2],self.vertices[0],
361	self.vertices[1]]
362	self.neighbors = [self.neighbors[2],self.neighbors[0],
363	self.neighbors[1]]
364
365	def rotate_longest_side(self):
366	self.rotate(self.longestSideID())
367
368	def getVertices(self):
369	return self.vertices
370
371	def get_vertices(self):
372	"""
373	Return a list of the vertices. The x and y values will be relative
374	Easting and Northings for the zone of the current geo_ref.
375	"""
376	return self.vertices
377
378	def calcArea(self):
379	ax = self.vertices[0].x
380	ay = self.vertices[0].y
381
382	bx = self.vertices[1].x
383	by = self.vertices[1].y
384
385	cx = self.vertices[2].x
386	cy = self.vertices[2].y
387
388	return abs((bxay-axby)+(cxby-bxcy)+(axcy-cxay))/2
389
390	def calcP(self):
391	#calculate the perimeter
392	ax = self.vertices[0].x
393	ay = self.vertices[0].y
394
395	bx = self.vertices[1].x
396	by = self.vertices[1].y
397
398	cx = self.vertices[2].x
399	cy = self.vertices[2].y
400
401	a = ((cx-bx)2+(cy-by)2)**0.5
402	b = ((ax-cx)2+(ay-cy)2)**0.5
403	c = ((bx-ax)2+(by-ay)2)**0.5
404
405	return a+b+c
406
407	def setNeighbors(self,neighbor1=None, neighbor2=None, neighbor3=None):
408	"""
409	neighbor1 is the triangle opposite vertex1 and so on.
410	Null represents no neighbor
411	"""
412	self.neighbors = [neighbor1, neighbor2, neighbor3]
413
414	def setAttribute(self,attribute):
415	"""
416	neighbor1 is the triangle opposite vertex1 and so on.
417	Null represents no neighbor
418	"""
419	self.attribute = attribute #this is a string
420
421	def __repr__(self):
422	return "[%s,%s]" % (self.vertices,self.attribute)
423
424
425	def draw(self, canvas, tags, scale=1, xoffset=0, yoffset=0,
426	colour="green"):
427	"""
428	Draw a triangle, returning the objectID
429	"""
430	return canvas.create_polygon(scale*(self.vertices[1].x + xoffset),
431	scale-1(self.vertices[1].y + yoffset),
432	scale*(self.vertices[0].x + xoffset),
433	scale-1(self.vertices[0].y + yoffset),
434	scale*(self.vertices[2].x + xoffset),
435	scale-1(self.vertices[2].y + yoffset),
436	tags = tags,
437	outline = colour,fill = '')
438
439	class Segment(MeshObject):
440	"""
441	Segments are edges whose presence in the triangulation is enforced.
442
443	"""
444	def __init__(self, vertex1, vertex2, tag = None ):
445	"""
446	Each segment is specified by listing the vertices of its endpoints
447	The vertices are Vertex objects
448	"""
449	assert(vertex1 != vertex2)
450	self.vertices = [vertex1,vertex2 ]
451
452	if tag is None:
453	self.tag = self.__class__.default
454	else:
455	self.tag = tag #this is a string
456
457	def __repr__(self):
458	return "[%s,%s]" % (self.vertices,self.tag)
459
460
461	def draw(self, canvas, tags,scale=1, xoffset=0, yoffset=0,colour='blue'):
462	x=[]
463	y=[]
464	for end in self.vertices:
465	#end.draw(canvas,scale, xoffset, yoffset ) # draw the vertices
466	x.append(scale*(end.x + xoffset))
467	y.append(-1scale(end.y + yoffset)) # - since for a canvas - is up
468
469	return canvas.create_line(x[0], y[0], x[1], y[1],
470	tags = tags,fill=colour)
471	def set_tag(self,tag):
472	self.tag = tag
473
474	# Class methods
475	def set_default_tag(cls, default):
476	cls.default = default
477
478	def get_default_tag(cls):
479	return cls.default
480
481	set_default_tag = classmethod(set_default_tag)
482	get_default_tag = classmethod(get_default_tag)
483
484	Segment.set_default_tag("")
485
486	class Rigid_triangulation:
487	"""
488	This is a triangulation that can't have triangles added or taken away.
489
490	It just represents the triangulation, not the mesh outline needed to
491	build the triangulation.
492
493	This is the guts of the data structure;
494	generated vertex list: [(x1,y1),(x2,y2),...] (Tuples of doubles)
495	generated segment list: [(point1,point2),(p3,p4),...]
496	(Tuples of integers)
497	generated segment tag list: [tag,tag,...] list of strings
498
499	generated triangle list: [(p1,p2,p3), (p4,p5,p6),....] tuple of points
500
501	generated triangle attribute list: [s1,s2,...] list of strings
502
503	generated triangle neighbor list: [(t1,t2,t3), (t4,t5,t6),....]
504	tuple of triangles
505
506	Should I do the basic structure like triangle, general
507	mesh or something else? How about like triangle, since
508	that's where the info is from, and it'll fit easier into
509	this file..
510
511	Removing the loners is difficult, since all the vert's
512	after it must be removed.
513
514	This happens in set_triangulation.
515
516	"""
517
518	def __init__(self,
519	triangles,
520	segments,
521	vertices,
522	triangle_tags,
523	triangle_neighbor,
524	segment_tags,
525	):
526
527	self.triangles=ensure_numeric(triangles)
528	self.triangle_neighbor=ensure_numeric(triangle_neighbor)
529	self.triangle_tags=triangle_tags
530	self.segments=ensure_numeric(segments)
531
532	self.segment_tags=segment_tags # string
533	self.vertices=ensure_numeric(vertices)
534
535
536
537	class Mesh:
538	"""
539	Representation of a 2D triangular mesh.
540	User attributes describe the mesh region/segments/vertices/attributes
541
542	mesh attributes describe the mesh that is produced eg triangles and
543	vertices.
544	All point information is relative to the geo_reference passed in
545
546
547	"""
548
549	def __repr__(self):
550	return """
551	mesh Triangles: %s
552	mesh Attribute Titles: %s
553	mesh Segments: %s
554	mesh Vertices: %s
555	user Segments: %s
556	user Vertices: %s
557	holes: %s
558	regions: %s""" % (self.meshTriangles,
559	self.attributeTitles,
560	self.meshSegments,
561	self.meshVertices,
562	self.getUserSegments(),
563	self.userVertices,
564	self.holes,
565	self.regions)
566
567	def __init__(self,
568	userSegments=None,
569	userVertices=None,
570	holes=None,
571	regions=None,
572	geo_reference=None):
573	self.meshTriangles=[]
574	self.attributeTitles=[]
575	self.meshSegments=[]
576	self.meshVertices=[]
577
578	self.setID={}
579	#a dictionary of names.
580	#multiple sets are allowed, but the gui does not yet
581	#support this
582
583	self.setID['None']=0
584	#contains the names of the sets pointing to the indexes
585	#in the list.
586
587	self.sets=[[]]
588	#Contains the lists of triangles (triangle sets)
589
590
591	self.visualise_graph = True
592
593	if userSegments is None:
594	self.userSegments=[]
595	else:
596	self.userSegments=userSegments
597	self.alphaUserSegments=[]
598
599	if userVertices is None:
600	self.userVertices=[]
601	else:
602	self.userVertices=userVertices
603
604	if holes is None:
605	self.holes=[]
606	else:
607	self.holes=holes
608
609	if regions is None:
610	self.regions=[]
611	else:
612	self.regions=regions
613
614	if geo_reference is None:
615	self.geo_reference = Geo_reference(DEFAULT_ZONE,0,0)
616	else:
617	self.geo_reference = geo_reference
618
619	self.shape = None
620
621	def __cmp__(self,other):
622
623	# A dic for the initial m
624	dic = self.Mesh2triangList()
625	dic_mesh = self.Mesh2MeshList()
626	for element in dic_mesh.keys():
627	dic[element] = dic_mesh[element]
628	for element in dic.keys():
629	dic[element].sort()
630
631	# A dic for the exported/imported m
632	dic_other = other.Mesh2triangList()
633	dic_mesh = other.Mesh2MeshList()
634	for element in dic_mesh.keys():
635	dic_other[element] = dic_mesh[element]
636	for element in dic.keys():
637	dic_other[element].sort()
638
639	#print "dsg************************8"
640	#print "dic ",dic
641	#print "*******8"
642	#print "mesh",dic_other
643	#print "dic.__cmp__(dic_o)",dic.__cmp__(dic_other)
644	#print "dsg************************8"
645
646	return (dic.__cmp__(dic_other))
647
648	def addUserPoint(self, pointType, x,y):
649	if pointType == Vertex:
650	point = self.addUserVertex(x,y)
651	if pointType == Hole:
652	point = self._addHole(x,y)
653	if pointType == Region:
654	point = self._addRegion(x,y)
655	return point
656
657	def addUserVertex(self, x,y):
658	v=Vertex(x, y)
659	self.userVertices.append(v)
660	return v
661
662	def _addHole(self, x,y):
663	h=Hole(x, y)
664	self.holes.append(h)
665	return h
666
667	def add_hole(self, x,y, geo_reference=None):
668	"""
669	adds a point, which represents a hole.
670
671	The point data can have it's own geo_refernece.
672	If geo_reference is None the data is asumed to be absolute
673	"""
674	[[x,y]] = self.geo_reference.change_points_geo_ref([x,y],
675	points_geo_ref=geo_reference)
676	return self._addHole(x, y)
677
678	def _addRegion(self, x,y):
679	h=Region(x, y)
680	self.regions.append(h)
681	return h
682
683	def add_region(self, x,y, geo_reference=None):
684	"""
685	adds a point, which represents a region.
686
687	The point data can have it's own geo_refernece.
688	If geo_reference is None the data is asumed to be absolute
689	"""
690	#FIXME: have the user set the tag and resolution here,
691	# but still return the instance, just in case.
692	[[x,y]] = self.geo_reference.change_points_geo_ref([x,y],
693	points_geo_ref=geo_reference)
694	return self._addRegion(x, y)
695
696	def build_grid(self, vert_rows, vert_columns):
697	"""
698	Build a grid with vert_rows number of vertex rows and
699	vert_columns number if vertex columns
700
701	Grid spacing of 1, the origin is the lower left hand corner.
702
703	FIXME(DSG-DSG) no test.
704	"""
705
706	for i in range(vert_rows):
707	for j in range(vert_columns):
708	self.addUserVertex(j,i)
709	self.auto_segment()
710	self.generateMesh(mode = "Q", minAngle=20.0)
711
712	# Depreciated
713	def addRegionEN(self, x,y):
714	print "depreciated, use add_region"
715	return self.add_region(x,y)
716
717
718	def add_vertices(self, point_data):
719	"""
720	Add user vertices.
721
722	The point_data can be a list of (x,y) values, a numeric
723	array or a geospatial_data instance.
724	"""
725	point_data = ensure_geospatial(point_data)
726	#print "point_data",point_data
727	# get points relative to the mesh geo_ref
728	points = point_data.get_data_points(geo_reference=self.geo_reference)
729
730	for point in points:
731	v=Vertex(point[0], point[1])
732	self.userVertices.append(v)
733
734	def add_hole_from_polygon(self, polygon, tags=None, geo_reference=None):
735	"""
736	Add a polygon with tags to the current mesh, as a region.
737	The maxArea of the region can be specified.
738
739	If a geo_reference of the polygon points is given, this is used.
740	If not;
741	The x,y info is assumed to be Easting and Northing, absolute,
742	for the meshes zone.
743
744	polygon a list of points, in meters that describe the polygon
745	(e.g. [[x1,y1],[x2,y2],...]
746	tags (e.g.{'wall':[0,1,3],'ocean':[2]})
747
748	This returns the region instance, so if the user whats to modify
749	it they can.
750	"""
751	return self._add_area_from_polygon(polygon, tags=tags,
752	hole=True,
753	geo_reference=geo_reference
754	)
755
756
757	def add_region_from_polygon(self, polygon, tags=None,
758	max_triangle_area=None, geo_reference=None):
759	"""
760	Add a polygon with tags to the current mesh, as a region.
761	The maxArea of the region can be specified.
762
763	If a geo_reference of the polygon points is given, this is used.
764	If not;
765	The x,y info is assumed to be Easting and Northing, absolute,
766	for the meshes zone.
767
768	polygon a list of points, in meters that describe the polygon
769	(e.g. [[x1,y1],[x2,y2],...]
770	tags (e.g.{'wall':[0,1,3],'ocean':[2]})
771
772	This returns the region instance (if a max_triangle_area is given),
773	so if the user whats to modify it they can.
774	"""
775	if max_triangle_area is None:
776	create_region = False
777	else:
778	create_region = True
779
780	region = self._add_area_from_polygon(polygon, tags=tags,
781	geo_reference=geo_reference,
782	region=create_region)
783	if max_triangle_area is not None:
784	region.setMaxArea(max_triangle_area)
785	return region
786
787
788	def _add_area_from_polygon(self, polygon, tags=None,
789	geo_reference=None,
790	hole=False,
791	region=False):
792	"""
793	Add a polygon with tags to the current mesh, as a region.
794	The maxArea of the region can be specified.
795
796	If a geo_reference of the polygon points is given, this is used.
797	If not;
798	The x,y info is assumed to be Easting and Northing, absolute,
799	for the meshes zone.
800
801	polygon a list of points, in meters that describe the polygon
802	(e.g. [[x1,y1],[x2,y2],...]
803	tags (e.g.{'wall':[0,1,3],'ocean':[2]})
804
805	This returns the region instance, so if the user whats to modify
806	it they can.
807
808	"""
809	#get absolute values
810	if geo_reference is not None:
811	polygon = geo_reference.get_absolute(polygon)
812	# polygon is now absolute
813	#print "polygon should be absolute",polygon
814
815	#create points, segs and tags
816	region_dict = {}
817	region_dict['points'] = polygon
818
819	#Create segments
820	#E.g. [[0,1], [1,2], [2,3], [3,0]]
821	#from polygon
822	#[0,1,2,3]
823	segments = []
824	N = len(polygon)
825	for i in range(N):
826	lo = i
827	hi = (lo + 1) % N
828	segments.append( [lo, hi] )
829	region_dict['segments'] = segments
830	region_dict['segment_tags'] = self._tag_dict2list(tags, N)
831
832
833	self.addVertsSegs(region_dict) #this is passing absolute values
834
835	if region is True:
836	#get inner point - absolute values
837	inner_point = point_in_polygon(polygon)
838	inner = self.add_region(inner_point[0], inner_point[1],
839	geo_reference=None)
840	elif hole is True:
841	#get inner point - absolute values
842	inner_point = point_in_polygon(polygon)
843	inner = self.add_hole(inner_point[0], inner_point[1],
844	geo_reference=None)
845	else:
846	inner = None
847
848	return inner
849
850	def _tag_dict2list(self, tags, number_of_segs):
851	"""
852	Convert a tag dictionary from this sort of format;
853	#{'wall':[0,3],'ocean':[2]}
854
855	To a list format
856	# ['wall', '', 'ocean', 'wall']
857
858	Note: '' is a default value of nothing
859	"""
860	# FIXME (DSG-DSG): Using '' as a default isn't good.
861	#Try None.
862	# Due to this default this method is too connected to
863	# _add_area_from_polygon
864
865	segment_tags = ['']*number_of_segs
866	if tags is not None:
867	for key in tags:
868	indices = tags[key]
869	for i in indices:
870	segment_tags[i] = key
871	return segment_tags
872
873	def add_circle(self, center, radius, segment_count=100,
874	center_geo_reference=None, tag = None,
875	region=False, hole=False):
876	"""
877	center is a point, in absulute or relative to center_geo_ref
878	radius is the radius of the circle
879	segment_count is the number of segments used to represent the circle
880	tag is a string name, given to the segments.
881	If region is True a region object is returned.
882	If hole is True a hole object is returned.
883	(Don't have them both True.)
884
885
886	"""
887	# convert center and radius to a polygon
888	cuts = []
889	factor = 2* math.pi/segment_count
890	for cut in range(segment_count):
891	cuts.append(cut*factor)
892
893	polygon = []
894	for cut in cuts:
895
896	x = center[0] + radius * math.cos(cut)
897	y = center[1] + radius * math.sin(cut)
898	polygon.append([x,y])
899	# build the tags
900	tags = {tag:range(segment_count)}
901
902	return self._add_area_from_polygon(polygon, tags=tags,
903	region=region, hole=hole,
904	geo_reference=center_geo_reference)
905
906	def auto_set_geo_reference(self):
907	"""
908	Automatically set the georeference, based in the minimum x and y
909	user vertex values.
910
911	Not to be used with the graphical interface
912
913	Not implemented.
914	Don't implement now. using the new georeferenced points class
915	will change this?
916	"""
917	#to do
918	# find the lower left hand corner
919	[xmin, ymin, xmax, ymax] = self.boxsize()
920
921	# set all points to that lower left hand corner.
922	# use change_geo_reference
923	new_geo = Geo_reference(self.geo_reference.get_zone(), xmin, ymin)
924	self.change_geo_reference(new_geo)
925
926	def change_geo_reference(self, new_geo_reference):
927	"""
928	Change from one geo_reference to another.
929	Not to be used with the graphical interface
930	"""
931	# FIXME
932	# change georeference of;
933	#self.userVertices = self.geo_reference.change_points_geo_ref( \
934	#self.userVertices)
935	#self.holes = self.geo_reference.change_points_geo_ref(self.holes)
936	#self.regions = self.geo_reference.change_points_geo_ref(self.regions)
937	# The above will not work.
938	# since userVertices (etc) is a list of point objects,
939	#not a list of lists.
940	# add a method to the points class to fix this up.
941
942	def add_segment(self, v1, v2, tag):
943	"""
944	Don't do this function.
945	what will the v's be objects? How is the user suppost to get them?
946
947	Indexes? If add vertstosegs or add_region_from_polygon is called
948	more than once then the actual index is not obvious. Making this
949	function confusing.
950	"""
951	pass
952
953
954	def add_points_and_segments(self, points,
955	segments, segment_tags = None):
956	"""
957	Add an outline of the mesh.
958	Vertices is a list of points/ a standard representation of points.
959	Segments is a list of tuples of integers. Each tuple defines the
960	start and end of the segment by it's vertex index, in relation to
961	the list of vertices.
962	segment_tags is an optional dictionary which is used to add tags to
963	the segments. The key is the tag name, value is the list of segment
964	indexes the tag will apply to.
965	eg. {'wall':[0,3],'ocean':[2]}
966
967	"""
968	#make sure the points are absolute
969	points = ensure_absolute(points)
970
971	#create points, segs and tags
972	region_dict = {}
973	region_dict['points'] = points
974	region_dict['segments'] = segments
975	region_dict['segment_tags'] = self._tag_dict2list(segment_tags,
976	len(segments))
977	self.addVertsSegs(region_dict)
978
979	def addVertsSegs(self, outlineDict):
980	"""
981	Add out-line (user Mesh) attributes given a dictionary of the lists
982	points: [(x1,y1),(x2,y2),...] (Tuples of doubles)
983	segments: [(point1,point2),(p3,p4),...] (Tuples of integers)
984	segment_tags: [S1Tag, S2Tag, ...] (list of strings)
985
986	Assume the values are in Eastings and Northings, with no reference
987	point. eg absolute
988	"""
989	if not outlineDict.has_key('segment_tags'):
990	outlineDict['segment_tags'] = []
991	for i in range(len(outlineDict['segments'])):
992	outlineDict['segment_tags'].append('')
993	#print "outlineDict['segment_tags']",outlineDict['segment_tags']
994	#print "outlineDict['points']",outlineDict['points']
995	#print "outlineDict['segments']",outlineDict['segments']
996
997	i_offset = len(self.userVertices)
998	#print "self.userVertices",self.userVertices
999	#print "index_offset",index_offset
1000	for point in outlineDict['points']:
1001	v=Vertex(point[0]-self.geo_reference.xllcorner,
1002	point[1]-self.geo_reference.yllcorner)
1003	self.userVertices.append(v)
1004
1005	for seg,seg_tag in map(None,outlineDict['segments'],
1006	outlineDict['segment_tags']):
1007	segObject = Segment(self.userVertices[int(seg[0])+i_offset],
1008	self.userVertices[int(seg[1])+i_offset] )
1009	if not seg_tag == '':
1010	segObject.set_tag(seg_tag)
1011	self.userSegments.append(segObject)
1012
1013
1014	def get_triangle_count(self):
1015	return len(self.getTriangulation())
1016
1017	def getUserVertices(self):
1018	"""
1019	Note: The x,y values will be relative to the mesh geo_reference
1020	This returns a list of vertex objects
1021	"""
1022	return self.userVertices
1023
1024	def get_user_vertices(self, absolute=True):
1025	"""
1026	This returns a list of (x,y) values
1027	(maybe it should return a geospatical object?
1028	It makes it a bit confusing though.)
1029	"""
1030	pointlist=[]
1031	for vertex in self.userVertices:
1032	pointlist.append([vertex.x,vertex.y])
1033	spat = Geospatial_data(pointlist, geo_reference=self.geo_reference)
1034	return spat.get_data_points(absolute=absolute)
1035
1036	def getUserSegments(self):
1037	allSegments = self.userSegments + self.alphaUserSegments
1038	#print "self.userSegments",self.userSegments
1039	#print "self.alphaUserSegments",self.alphaUserSegments
1040	#print "allSegments",allSegments
1041	return allSegments
1042
1043	def deleteUserSegments(self,seg):
1044	if self.userSegments.count(seg) == 0:
1045	self.alphaUserSegments.remove(seg)
1046	pass
1047	else:
1048	self.userSegments.remove(seg)
1049
1050	def clearUserSegments(self):
1051	self.userSegments = []
1052	self.alphaUserSegments = []
1053
1054	def getTriangulation(self):
1055	return self.meshTriangles
1056
1057	def getMeshVertices(self):
1058	return self.meshVertices
1059
1060	def getMeshSegments(self):
1061	return self.meshSegments
1062
1063	def getHoles(self):
1064	return self.holes
1065
1066	def getRegions(self):
1067	return self.regions
1068
1069	def isTriangulation(self):
1070	if self.meshVertices == []:
1071	return False
1072	else:
1073	return True
1074
1075	def addUserSegment(self, v1,v2):
1076	"""
1077	PRECON: A segment between the two vertices is not already present.
1078	Check by calling isUserSegmentNew before calling this function.
1079
1080	"""
1081	s=Segment( v1,v2)
1082	self.userSegments.append(s)
1083	return s
1084
1085	def generate_mesh(self,
1086	maximum_triangle_area="",
1087	minimum_triangle_angle=28.0,
1088	verbose=True):
1089	if verbose is True:
1090	silent = ''
1091	else:
1092	silent = 'Q'
1093	self.generateMesh(mode = silent +"pzq"+str(minimum_triangle_angle)
1094	+"a"+str(maximum_triangle_area)
1095	+"a")
1096	#The last a is so areas for regions will be used
1097
1098	def generateMesh(self, mode = None, maxArea = None, minAngle= None,
1099	isRegionalMaxAreas = True):
1100	"""
1101	Based on the current user vaules, holes and regions
1102	generate a new mesh
1103	mode is a string that sets conditions on the mesh generations
1104	see triangle_instructions.txt for a definition of the commands
1105	PreCondition: maxArea is a double
1106	"""
1107	#print "mode ",mode
1108	if mode == None:
1109	self.mode = ""
1110	else:
1111	self.mode = mode
1112
1113	if self.mode.find('p') < 0:
1114	self.mode += 'p' #p - read a planar straight line graph.
1115	#there must be segments to use this switch
1116	# TODO throw an aception if there aren't seg's
1117	# it's more comlex than this. eg holes
1118	if self.mode.find('z') < 0:
1119	self.mode += 'z' # z - Number all items starting from zero
1120	# (rather than one)
1121	if self.mode.find('n'):
1122	self.mode += 'n' # n - output a list of neighboring triangles
1123	if self.mode.find('A') < 0:
1124	self.mode += 'A' # A - output region attribute list for triangles
1125
1126	if not self.mode.find('V') and not self.mode.find('Q'):
1127	self.mode += 'V' # V - output info about what Triangle is doing
1128
1129	if self.mode.find('q') < 0 and minAngle is not None:
1130	# print "********8minAngle****** ",minAngle
1131	min_angle = 'q' + str(minAngle)
1132	self.mode += min_angle # z - Number all items starting from zero
1133	# (rather than one)
1134	if maxArea != None:
1135	self.mode += 'a' + str(maxArea)
1136	#FIXME (DSG-DSG) This isn't explained.
1137	if isRegionalMaxAreas:
1138	self.mode += 'a'
1139	#print "mesh#generateMesh# self.mode",self.mode
1140	meshDict = self.Mesh2triangList()
1141
1142	#FIXME (DSG-DSG) move below section into generate_mesh.py
1143	# & 4 functions eg segment_strings2ints
1144	#print "*************************!@!@ This is going to triangle !@!@"
1145	#print meshDict
1146	#print "************************!@!@ This is going to triangle !@!@"
1147
1148	#print "meshDict['segmenttaglist']", meshDict['segmenttaglist']
1149	[meshDict['segmenttaglist'],
1150	segconverter] = segment_strings2ints(meshDict['segmenttaglist'],
1151	initialconversions)
1152	#print "regionlist",meshDict['regionlist']
1153	[meshDict['regionlist'],
1154	regionconverter] = region_strings2ints(meshDict['regionlist'])
1155	#print "%%%%%%%%%%%%%%%%%%%%%%%%%%%regionlist",meshDict['regionlist']
1156	#print "meshDict['segmenttaglist']", meshDict['segmenttaglist'
1157	generatedMesh = generate_mesh(
1158	meshDict['pointlist'],
1159	meshDict['segmentlist'],
1160	meshDict['holelist'],
1161	meshDict['regionlist'],
1162	meshDict['pointattributelist'],
1163	meshDict['segmenttaglist'],
1164	self.mode,
1165	meshDict['pointlist'])
1166	#print "%%%%%%%%%%%%%%%%%%%%%%%%%%%generated",generatedMesh
1167	generatedMesh['qaa'] = 1
1168	generatedMesh['generatedsegmentmarkerlist'] = \
1169	segment_ints2strings(generatedMesh['generatedsegmentmarkerlist'],
1170	segconverter)
1171	#print "processed gen",generatedMesh['generatedsegmentmarkerlist']
1172	generatedMesh['generatedtriangleattributelist'] = \
1173	region_ints2strings(generatedMesh['generatedtriangleattributelist'],
1174	regionconverter)
1175
1176	#FIXME (DSG-DSG) move above section into generate_mesh.py
1177
1178	if len(generatedMesh['generatedpointattributelist'][0])==0:
1179	self.attributeTitles = []
1180	generatedMesh['generatedpointattributetitlelist']= \
1181	self.attributeTitles
1182	#print "################ FROM TRIANGLE"
1183	#print "generatedMesh",generatedMesh
1184	#print "##################"
1185
1186	self.setTriangulation(generatedMesh)
1187
1188	def clearTriangulation(self):
1189
1190	#Clear the current generated mesh values
1191	self.meshTriangles=[]
1192	self.meshSegments=[]
1193	self.meshVertices=[]
1194
1195	def removeDuplicatedUserVertices(self):
1196	"""Pre-condition: There are no user segments
1197	This function will keep the first duplicate
1198	"""
1199	assert self.getUserSegments() == []
1200	self.userVertices, counter = self.removeDuplicatedVertices(
1201	self.userVertices)
1202	return counter
1203
1204	def removeDuplicatedVertices(self, Vertices):
1205	"""
1206	This function will keep the first duplicate, remove all others
1207	Precondition: Each vertex has a dupindex, which is the list
1208	index.
1209
1210	Note: this removes vertices that have the same x,y values,
1211	not duplicate instances in the Vertices list.
1212	"""
1213	remove = []
1214	index = 0
1215	for v in Vertices:
1216	v.dupindex = index
1217	index += 1
1218	t = list(Vertices)
1219	t.sort(Point.cmp_xy)
1220
1221	length = len(t)
1222	behind = 0
1223	ahead = 1
1224	counter = 0
1225	while ahead < length:
1226	b = t[behind]
1227	ah = t[ahead]
1228	if (b.y == ah.y and b.x == ah.x):
1229	remove.append(ah.dupindex)
1230	behind += 1
1231	ahead += 1
1232
1233	# remove the duplicate vertices
1234	remove.sort()
1235	remove.reverse()
1236	for i in remove:
1237	Vertices.pop(i)
1238	pass
1239
1240	#Remove the attribute that this function added
1241	for v in Vertices:
1242	del v.dupindex
1243	return Vertices,counter
1244
1245	# FIXME (DSG-DSG) Move this to geospatial
1246	def thinoutVertices(self, delta):
1247	"""Pre-condition: There are no user segments
1248	This function will keep the first duplicate
1249	"""
1250	assert self.getUserSegments() == []
1251	#t = self.userVertices
1252	#self.userVertices =[]
1253	boxedVertices = {}
1254	thinnedUserVertices =[]
1255	delta = round(delta,1)
1256
1257	for v in self.userVertices :
1258	# tag is the center of the boxes
1259	tag = (round(v.x/delta,0)delta,round(v.y/delta,0)delta)
1260	#this creates a dict of lists of faces, indexed by tag
1261	boxedVertices.setdefault(tag,[]).append(v)
1262
1263	for [tag,verts] in boxedVertices.items():
1264	min = delta
1265	bestVert = None
1266	tagVert = Vertex(tag[0],tag[1])
1267	for v in verts:
1268	dist = v.DistanceToPoint(tagVert)
1269	if (dist<min):
1270	min = dist
1271	bestVert = v
1272	thinnedUserVertices.append(bestVert)
1273	self.userVertices =thinnedUserVertices
1274
1275
1276	def isUserSegmentNew(self, v1,v2):
1277	identicalSegs= [x for x in self.getUserSegments() \
1278	if (x.vertices[0] == v1 and x.vertices[1] == v2)
1279	or (x.vertices[0] == v2 and x.vertices[1] == v1) ]
1280
1281	return len(identicalSegs) == 0
1282
1283
1284	def deleteSegsOfVertex(self, delVertex):
1285	"""
1286	Delete this vertex and any segments that connect to it.
1287	"""
1288	#Find segments that connect to delVertex
1289	deletedSegments = []
1290	for seg in self.getUserSegments():
1291	if (delVertex in seg.vertices):
1292	deletedSegments.append(seg)
1293	# Delete segments that connect to delVertex
1294	for seg in deletedSegments:
1295	self.deleteUserSegments(seg)
1296	return deletedSegments
1297
1298
1299	def deleteMeshObject(self, MeshObject):
1300	"""
1301	Returns a list of all objects that were removed
1302	"""
1303	deletedObs = []
1304	if isinstance(MeshObject, Vertex ):
1305	deletedObs = self.deleteSegsOfVertex(MeshObject)
1306	deletedObs.append(MeshObject)
1307	self.userVertices.remove(MeshObject)
1308	elif isinstance(MeshObject, Segment):
1309	deletedObs.append(MeshObject)
1310	self.deleteUserSegments(MeshObject)
1311	elif isinstance(MeshObject, Hole):
1312	deletedObs.append(MeshObject)
1313	self.holes.remove(MeshObject)
1314	elif isinstance(MeshObject, Region):
1315	deletedObs.append(MeshObject)
1316	self.regions.remove(MeshObject)
1317	return deletedObs
1318
1319	def Mesh2triangList(self, userVertices=None,
1320	userSegments=None,
1321	holes=None,
1322	regions=None):
1323	"""
1324	Convert the Mesh to a dictionary of the lists needed for the
1325	triang module
1326	points list: [(x1,y1),(x2,y2),...] (Tuples of doubles)
1327	pointattributelist: [(a11,a12,...),(a21,a22),...] (Tuples of doubles)
1328	segment list: [(point1,point2),(p3,p4),...] (Tuples of integers)
1329	hole list: [(x1,y1),...](Tuples of doubles, one inside each hole)
1330	regionlist: [ (x1,y1,tag, max area),...] (Tuple of 3-4 doubles)
1331
1332	Note, this adds an index attribute to the user Vertex objects.
1333
1334	Used to produce output to triangle
1335	"""
1336	if userVertices is None:
1337	userVertices = self.getUserVertices()
1338	if userSegments is None:
1339	userSegments = self.getUserSegments()
1340	if holes is None:
1341	holes = self.getHoles()
1342	if regions is None:
1343	regions = self.getRegions()
1344
1345	meshDict = {}
1346
1347	pointlist=[]
1348	pointattributelist=[]
1349	index = 0
1350	for vertex in userVertices:
1351	vertex.index = index
1352	pointlist.append((vertex.x,vertex.y))
1353	pointattributelist.append((vertex.attributes))
1354
1355	index += 1
1356	meshDict['pointlist'] = pointlist
1357	meshDict['pointattributelist'] = pointattributelist
1358
1359	segmentlist=[]
1360	segmenttaglist=[]
1361	for seg in userSegments:
1362	segmentlist.append((seg.vertices[0].index,seg.vertices[1].index))
1363	segmenttaglist.append(seg.tag)
1364	meshDict['segmentlist'] =segmentlist
1365	meshDict['segmenttaglist'] =segmenttaglist
1366
1367	holelist=[]
1368	for hole in holes:
1369	holelist.append((hole.x,hole.y))
1370	meshDict['holelist'] = holelist
1371
1372	regionlist=[]
1373	for region in regions:
1374	if (region.getMaxArea() != None):
1375	regionlist.append((region.x,region.y,region.getTag(),
1376	region.getMaxArea()))
1377	else:
1378	regionlist.append((region.x,region.y,region.getTag()))
1379	meshDict['regionlist'] = regionlist
1380	#print "(("
1381	#print meshDict
1382	#print meshDict['regionlist']
1383	#print "(("
1384	return meshDict
1385
1386	def Mesh2MeshList(self):
1387	"""
1388	Convert the Mesh to a dictionary of lists describing the
1389	triangulation variables;
1390	generated point list: [(x1,y1),(x2,y2),...] (Tuples of doubles)
1391	generated point attribute list: [(a11,a12,...),(a21,a22),...]
1392	(Tuples of doubles)
1393	generated point attribute title list:[A1Title, A2Title ...]
1394	(list of strings)
1395	generated segment list: [(point1,point2),(p3,p4),...]
1396	(Tuples of integers)
1397	generated segment tag list: [tag,tag,...] list of strings
1398
1399	generated triangle list: [(p1,p2,p3), (p4,p5,p6),....] tuple of points
1400
1401	generated triangle attribute list: [s1,s2,...] list of strings
1402
1403	generated triangle neighbor list: [(t1,t2,t3), (t4,t5,t6),....]
1404	tuple of triangles
1405
1406	Used to produce .tsh file
1407	"""
1408
1409	meshDict = {}
1410	pointlist=[]
1411	pointattributelist=[]
1412
1413
1414	self.maxVertexIndex=0
1415	for vertex in self.meshVertices:
1416	vertex.index = self.maxVertexIndex
1417	pointlist.append((vertex.x,vertex.y))
1418	pointattributelist.append((vertex.attributes))
1419	self.maxVertexIndex += 1
1420
1421	meshDict['generatedpointlist'] = pointlist
1422	meshDict['generatedpointattributelist'] = pointattributelist
1423	meshDict['generatedpointattributetitlelist'] = self.attributeTitles
1424	#segments
1425	segmentlist=[]
1426	segmenttaglist=[]
1427	for seg in self.meshSegments:
1428	segmentlist.append((seg.vertices[0].index,seg.vertices[1].index))
1429	segmenttaglist.append(seg.tag)
1430	meshDict['generatedsegmentlist'] =segmentlist
1431	meshDict['generatedsegmenttaglist'] =segmenttaglist
1432
1433	# Make sure that the indexation is correct
1434	index = 0
1435	for tri in self.meshTriangles:
1436	tri.index = index
1437	index += 1
1438
1439	trianglelist = []
1440	triangleattributelist = []
1441	triangleneighborlist = []
1442	for tri in self.meshTriangles:
1443	trianglelist.append((tri.vertices[0].index,tri.vertices[1].index,
1444	tri.vertices[2].index))
1445	triangleattributelist.append([tri.attribute])
1446	neighborlist = [-1,-1,-1]
1447	for neighbor,index in map(None,tri.neighbors,
1448	range(len(tri.neighbors))):
1449	if neighbor:
1450	neighborlist[index] = neighbor.index
1451	triangleneighborlist.append(neighborlist)
1452
1453	meshDict['generatedtrianglelist'] = trianglelist
1454	meshDict['generatedtriangleattributelist'] = triangleattributelist
1455	meshDict['generatedtriangleneighborlist'] = triangleneighborlist
1456
1457	#print "mesh.Mesh2MeshList))"
1458	#print meshDict
1459	#print "mesh.Mesh2MeshList))"
1460
1461	return meshDict
1462
1463
1464	def Mesh2MeshDic(self):
1465	"""
1466	Convert the user and generated info of a mesh to a dictionary
1467	structure
1468	"""
1469	dic = self.Mesh2triangList()
1470	dic_mesh = self.Mesh2MeshList()
1471	for element in dic_mesh.keys():
1472	dic[element] = dic_mesh[element]
1473	return dic
1474
1475	def setTriangulation(self, genDict):
1476	"""
1477	Set the mesh attributes given a dictionary of the lists
1478	returned from the triang module
1479	generated point list: [(x1,y1),(x2,y2),...] (Tuples of doubles)
1480	generated point attribute list:[(P1att1,P1attt2, ...),
1481	(P2att1,P2attt2,...),...]
1482	generated point attribute title list:[A1Title, A2Title ...]
1483	(list of strings)
1484	generated segment list: [(point1,point2),(p3,p4),...]
1485	(Tuples of integers)
1486	generated segment marker list: [S1Tag, S2Tag, ...] (list of ints)
1487	triangle list: [(point1,point2, point3),(p5,p4, p1),...]
1488	(Tuples of integers)
1489	triangle neighbor list: [(triangle1,triangle2, triangle3),
1490	(t5,t4, t1),...] (Tuples of integers)
1491	-1 means there's no triangle neighbor
1492	triangle attribute list: [(T1att), (T2att), ...]
1493
1494	"""
1495	#Clear the current generated mesh values
1496	self.meshTriangles=[]
1497	self.attributeTitles=[]
1498	self.meshSegments=[]
1499	self.meshVertices=[]
1500
1501	#print "mesh.setTriangulation@#@#@#"
1502	#print genDict
1503	#print "@#@#@#"
1504
1505	self.maxVertexIndex = 0
1506	for point in genDict['generatedpointlist']:
1507	v=Vertex(point[0], point[1])
1508	v.index = self.maxVertexIndex
1509	self.maxVertexIndex +=1
1510	self.meshVertices.append(v)
1511
1512	self.attributeTitles = genDict['generatedpointattributetitlelist']
1513
1514	index = 0
1515	for seg,marker in map(None,genDict['generatedsegmentlist'],
1516	genDict['generatedsegmentmarkerlist']):
1517	segObject = Segment( self.meshVertices[int(seg[0])],
1518	self.meshVertices[int(seg[1])], tag = marker )
1519	segObject.index = index
1520	index +=1
1521	self.meshSegments.append(segObject)
1522
1523	index = 0
1524	for triangle in genDict['generatedtrianglelist']:
1525	tObject =Triangle( self.meshVertices[int(triangle[0])],
1526	self.meshVertices[int(triangle[1])],
1527	self.meshVertices[int(triangle[2])] )
1528	tObject.index = index
1529	index +=1
1530	self.meshTriangles.append(tObject)
1531
1532	index = 0
1533	for att in genDict['generatedtriangleattributelist']:
1534	if att == []:
1535	self.meshTriangles[index].setAttribute("")
1536	else:
1537	self.meshTriangles[index].setAttribute(att[0])
1538	index += 1
1539
1540	index = 0
1541	for att in genDict['generatedpointattributelist']:
1542	if att == None:
1543	self.meshVertices[index].setAttributes([])
1544	else:
1545	self.meshVertices[index].setAttributes(att)
1546	index += 1
1547
1548	index = 0
1549	for triangle in genDict['generatedtriangleneighborlist']:
1550	# Build a list of triangle object neighbors
1551	ObjectNeighbor = []
1552	for neighbor in triangle:
1553	if ( neighbor != -1):
1554	ObjectNeighbor.append(self.meshTriangles[int(neighbor)])
1555	else:
1556	ObjectNeighbor.append(None)
1557	self.meshTriangles[index].setNeighbors(ObjectNeighbor[0],
1558	ObjectNeighbor[1],
1559	ObjectNeighbor[2])
1560	index += 1
1561	#genDict['lonepointlist'].sort()
1562	#genDict['lonepointlist'].reverse()
1563	#for loner in genDict['lonepointlist']:
1564	# Remove the loner vertex
1565	#print "Removing the loner", loner
1566	#self.meshVertices.pop(loner)
1567	# delete the info from this data structure
1568	del genDict['generatedtriangleneighborlist']
1569	del genDict['generatedpointmarkerlist']
1570	del genDict['generatedpointlist']
1571	del genDict['generatedpointattributetitlelist']
1572	del genDict['generatedsegmentlist']
1573	del genDict['generatedsegmentmarkerlist']
1574	del genDict['generatedtrianglelist']
1575	del genDict['generatedtriangleattributelist']
1576	del genDict['generatedpointattributelist']
1577	genDict['zpp'] = 1 # just an identifier
1578
1579	def setMesh(self, genDict):
1580	"""
1581	Set the user Mesh attributes given a dictionary of the lists
1582	point list: [(x1,y1),(x2,y2),...] (Tuples of doubles)
1583	point attribute list:[(P1att1,P1attt2, ...),(P2att1,P2attt2,...),...]
1584	segment list: [(point1,point2),(p3,p4),...] (Tuples of integers)
1585	segment tag list: [S1Tag, S2Tag, ...] (list of ints)
1586	region list: [(x1,y1),(x2,y2),...] (Tuples of doubles)
1587	region attribute list: ["","reservoir",""] list of strings
1588	region max area list:[real, None, Real,...] list of None and reals
1589
1590	mesh is an instance of a mesh object
1591	"""
1592	#Clear the current user mesh values
1593	self.clearUserSegments()
1594	self.userVertices=[]
1595	self.Holes=[]
1596	self.Regions=[]
1597
1598	#print "mesh.setMesh@#@#@#"
1599	#print genDict
1600	#print "@#@#@#"
1601
1602	#index = 0
1603	for point in genDict['pointlist']:
1604	v=Vertex(point[0], point[1])
1605	#v.index = index
1606	#index +=1
1607	self.userVertices.append(v)
1608
1609	#index = 0
1610	for seg,tag in map(None,genDict['segmentlist'],
1611	genDict['segmenttaglist']):
1612	segObject = Segment( self.userVertices[int(seg[0])],
1613	self.userVertices[int(seg[1])], tag = tag )
1614	#segObject.index = index
1615	#index +=1
1616	self.userSegments.append(segObject)
1617
1618	# Remove the loading of attribute info.
1619	# Have attribute info added using least_squares in pyvolution
1620	# index = 0
1621	# for att in genDict['pointattributelist']:
1622	# if att == None:
1623	# self.userVertices[index].setAttributes([])
1624	# else:
1625	# self.userVertices[index].setAttributes(att)
1626	# index += 1
1627
1628	#index = 0
1629	for point in genDict['holelist']:
1630	h=Hole(point[0], point[1])
1631	#h.index = index
1632	#index +=1
1633	self.holes.append(h)
1634
1635	#index = 0
1636	for reg,att,maxArea in map(None,
1637	genDict['regionlist'],
1638	genDict['regionattributelist'],
1639	genDict['regionmaxarealist']):
1640	Object = Region( reg[0],
1641	reg[1],
1642	tag = att,
1643	maxArea = maxArea)
1644	#Object.index = index
1645	#index +=1
1646	self.regions.append(Object)
1647
1648	def Testauto_segment(self):
1649	newsegs = []
1650	s1 = Segment(self.userVertices[0],
1651	self.userVertices[1])
1652	s2 = Segment(self.userVertices[0],
1653	self.userVertices[2])
1654	s3 = Segment(self.userVertices[2],
1655	self.userVertices[1])
1656	if self.isUserSegmentNew(s1.vertices[0],s1.vertices[1]):
1657	newsegs.append(s1)
1658	if self.isUserSegmentNew(s2.vertices[0],s2.vertices[1]):
1659	newsegs.append(s2)
1660	if self.isUserSegmentNew(s3.vertices[0],s3.vertices[1]):
1661	newsegs.append(s3)
1662	#DSG!!!
1663	self.userSegments.extend(newsegs)
1664	return newsegs
1665
1666
1667	def savePickle(self, currentName):
1668	fd = open(currentName, 'w')
1669	pickle.dump(self,fd)
1670	fd.close()
1671
1672	def auto_segmentFilter(self,raw_boundary=True,
1673	remove_holes=False,
1674	smooth_indents=False,
1675	expand_pinch=False):
1676	"""
1677	Change the filters applied on the alpha shape boundary
1678	"""
1679	if self.shape is None:
1680	return [],[],0.0
1681	return self._boundary2mesh(raw_boundary=raw_boundary,
1682	remove_holes=remove_holes,
1683	smooth_indents=smooth_indents,
1684	expand_pinch=expand_pinch)
1685
1686
1687
1688	def auto_segment(self, alpha = None,
1689	raw_boundary=True,
1690	remove_holes=False,
1691	smooth_indents=False,
1692	expand_pinch=False):
1693	"""
1694	Precon: There must be 3 or more vertices in the userVertices structure
1695	"""
1696	self._createBoundary(alpha=alpha)
1697	return self._boundary2mesh(raw_boundary=raw_boundary,
1698	remove_holes=remove_holes,
1699	smooth_indents=smooth_indents,
1700	expand_pinch=expand_pinch)
1701
1702	def _createBoundary(self,alpha=None):
1703	"""
1704	"""
1705	points=[]
1706	for vertex in self.getUserVertices():
1707	points.append((vertex.x,vertex.y))
1708	self.shape = anuga.alpha_shape.alpha_shape.Alpha_Shape(points,
1709	alpha=alpha)
1710
1711
1712	def _boundary2mesh(self, raw_boundary=True,
1713	remove_holes=False,
1714	smooth_indents=False,
1715	expand_pinch=False):
1716	"""
1717	Precon there must be a shape object.
1718	"""
1719	self.shape.set_boundary_type(raw_boundary=raw_boundary,
1720	remove_holes=remove_holes,
1721	smooth_indents=smooth_indents,
1722	expand_pinch=expand_pinch)
1723	boundary_segs = self.shape.get_boundary()
1724	#print "boundary_segs",boundary_segs
1725	segs2delete = self.alphaUserSegments
1726	#FIXME(DSG-DSG) this algorithm needs comments
1727	new_segs = {}
1728	#alpha_segs = []
1729	#user_segs = []
1730	for seg in boundary_segs:
1731	v1 = self.userVertices[int(seg[0])]
1732	v2 = self.userVertices[int(seg[1])]
1733	boundary_seg = Segment(v1, v2)
1734	new_segs[(v1,v2)] = boundary_seg
1735
1736	for user_seg in self.userSegments:
1737	if new_segs.has_key((user_seg.vertices[0],
1738	user_seg.vertices[1])):
1739	del new_segs[user_seg.vertices[0],
1740	user_seg.vertices[1]]
1741	elif new_segs.has_key((user_seg.vertices[1],
1742	user_seg.vertices[0])):
1743	del new_segs[user_seg.vertices[1],
1744	user_seg.vertices[0]]
1745
1746	optimum_alpha = self.shape.get_alpha()
1747	alpha_segs_no_user_segs = new_segs.values()
1748	self.alphaUserSegments = alpha_segs_no_user_segs
1749	return alpha_segs_no_user_segs, segs2delete, optimum_alpha
1750
1751	def _boundary2mesh_old(self, raw_boundary=True,
1752	remove_holes=False,
1753	smooth_indents=False,
1754	expand_pinch=False):
1755	"""
1756	Precon there must be a shape object.
1757	"""
1758	self.shape.set_boundary_type(raw_boundary=raw_boundary,
1759	remove_holes=remove_holes,
1760	smooth_indents=smooth_indents,
1761	expand_pinch=expand_pinch)
1762	boundary_segs = self.shape.get_boundary()
1763	#print "boundary_segs",boundary_segs
1764	segs2delete = self.alphaUserSegments
1765
1766	#FIXME(DSG-DSG) this algorithm needs comments
1767	#FIXME(DSG-DSG) can it be sped up? It's slow
1768	new_segs = []
1769	alpha_segs = []
1770	user_segs = []
1771	for seg in boundary_segs:
1772	v1 = self.userVertices[int(seg[0])]
1773	v2 = self.userVertices[int(seg[1])]
1774	alpha_seg = self.representedAlphaUserSegment(v1, v2)
1775	user_seg = self.representedUserSegment(v1, v2)
1776	#DSG!!!
1777	assert not(not (alpha_seg == None) and not (user_seg == None))
1778	if not alpha_seg == None:
1779	alpha_segs.append(alpha_seg)
1780	elif not user_seg == None:
1781	user_segs.append(user_seg)
1782	else:
1783	unique_seg = Segment(v1, v2)
1784	new_segs.append(unique_seg)
1785
1786	for seg in alpha_segs:
1787	try:
1788	segs2delete.remove(seg)
1789	except:
1790	pass
1791
1792	self.alphaUserSegments = []
1793	self.alphaUserSegments.extend(new_segs)
1794	self.alphaUserSegments.extend(alpha_segs)
1795
1796	optimum_alpha = self.shape.get_alpha()
1797	# need to draw newsegs
1798	return new_segs, segs2delete, optimum_alpha
1799
1800	def representedAlphaUserSegment(self, v1,v2):
1801	identicalSegs= [x for x in self.alphaUserSegments \
1802	if (x.vertices[0] == v1 and x.vertices[1] == v2)
1803	or (x.vertices[0] == v2 and x.vertices[1] == v1) ]
1804
1805	if identicalSegs == []:
1806	return None
1807	else:
1808	# Only return the first one.
1809	return identicalSegs[0]
1810
1811	def representedUserSegment(self, v1,v2):
1812	identicalSegs= [x for x in self.userSegments \
1813	if (x.vertices[0] == v1 and x.vertices[1] == v2)
1814	or (x.vertices[0] == v2 and x.vertices[1] == v1) ]
1815
1816	if identicalSegs == []:
1817	return None
1818	else:
1819	# Only return the first one.
1820	return identicalSegs[0]
1821
1822	def joinVertices(self):
1823	"""
1824	Return list of segments connecting all userVertices
1825	in the order they were given
1826
1827	Precon: There must be 3 or more vertices in the userVertices structure
1828	"""
1829
1830	newsegs = []
1831
1832	v1 = self.userVertices[0]
1833	for v2 in self.userVertices[1:]:
1834	if self.isUserSegmentNew(v1,v2):
1835	newseg = Segment(v1, v2)
1836	newsegs.append(newseg)
1837	v1 = v2
1838
1839	#Connect last point to the first
1840	v2 = self.userVertices[0]
1841	if self.isUserSegmentNew(v1,v2):
1842	newseg = Segment(v1, v2)
1843	newsegs.append(newseg)
1844
1845
1846	#Update list of user segments
1847	#DSG!!!
1848	self.userSegments.extend(newsegs)
1849	return newsegs
1850
1851	def normaliseMesh(self,scale, offset, height_scale):
1852	[xmin, ymin, xmax, ymax] = self.boxsize()
1853	[attmin0, attmax0] = self.maxMinVertAtt(0)
1854	#print "[attmin0, attmax0]" ,[attmin0, attmax0]
1855	[attmin1, attmax1] = self.maxMinVertAtt(1)
1856	#print [xmin, ymin, xmax, ymax]
1857	xrange = xmax - xmin
1858	yrange = ymax - ymin
1859	if xrange > yrange:
1860	min,max = xmin, xmax
1861	else:
1862	min,max = ymin, ymax
1863
1864	for obj in self.getUserVertices():
1865	obj.x = (obj.x - xmin)/(max- min)*scale + offset
1866	obj.y = (obj.y - ymin)/(max- min)*scale + offset
1867	if len(obj.attributes) > 0 and attmin0 != attmax0:
1868	obj.attributes[0] = (obj.attributes[0]-attmin0)/ \
1869	(attmax0-attmin0)*height_scale
1870	if len(obj.attributes) > 1 and attmin1 != attmax1:
1871	obj.attributes[1] = (obj.attributes[1]-attmin1)/ \
1872	(attmax1-attmin1)*height_scale
1873
1874	for obj in self.getMeshVertices():
1875	obj.x = (obj.x - xmin)/(max- min)*scale + offset
1876	obj.y = (obj.y - ymin)/(max- min)*scale + offset
1877	if len(obj.attributes) > 0 and attmin0 != attmax0:
1878	obj.attributes[0] = (obj.attributes[0]-attmin0)/ \
1879	(attmax0-attmin0)*height_scale
1880	if len(obj.attributes) > 1 and attmin1 != attmax1:
1881	obj.attributes[1] = (obj.attributes[1]-attmin1)/ \
1882	(attmax1-attmin1)*height_scale
1883
1884	for obj in self.getHoles():
1885	obj.x = (obj.x - xmin)/(max- min)*scale + offset
1886	obj.y = (obj.y - ymin)/(max- min)*scale + offset
1887	for obj in self.getRegions():
1888	obj.x = (obj.x - xmin)/(max- min)*scale + offset
1889	obj.y = (obj.y - ymin)/(max- min)*scale + offset
1890	[xmin, ymin, xmax, ymax] = self.boxsize()
1891	#print [xmin, ymin, xmax, ymax]
1892
1893	def boxsizeVerts(self):
1894	"""
1895	Returns a list of verts denoting a box or triangle that contains
1896	verts on the xmin, ymin, xmax and ymax axis.
1897	Structure: list of verts
1898	"""
1899
1900	large = kinds.default_float_kind.MAX
1901	xmin= large
1902	xmax=-large
1903	ymin= large
1904	ymax=-large
1905	for vertex in self.userVertices:
1906	if vertex.x < xmin:
1907	xmin = vertex.x
1908	xminVert = vertex
1909	if vertex.x > xmax:
1910	xmax = vertex.x
1911	xmaxVert = vertex
1912
1913	if vertex.y < ymin:
1914	ymin = vertex.y
1915	yminVert = vertex
1916	if vertex.y > ymax:
1917	ymax = vertex.y
1918	ymaxVert = vertex
1919	verts, count = self.removeDuplicatedVertices([xminVert,
1920	xmaxVert,
1921	yminVert,
1922	ymaxVert])
1923
1924	return verts
1925
1926	def boxsize(self):
1927	"""
1928	Returns a list denoting a box that contains the entire
1929	structure of vertices
1930	Structure: [xmin, ymin, xmax, ymax]
1931	"""
1932
1933	large = kinds.default_float_kind.MAX
1934	xmin= large
1935	xmax=-large
1936	ymin= large
1937	ymax=-large
1938	for vertex in self.userVertices:
1939	if vertex.x < xmin:
1940	xmin = vertex.x
1941	if vertex.x > xmax:
1942	xmax = vertex.x
1943
1944	if vertex.y < ymin:
1945	ymin = vertex.y
1946	if vertex.y > ymax:
1947	ymax = vertex.y
1948	return [xmin, ymin, xmax, ymax]
1949
1950	def maxMinVertAtt(self, iatt):
1951	"""
1952	Returns a list denoting a box that contains the entire structure
1953	of vertices
1954	Structure: [xmin, ymin, xmax, ymax]
1955	"""
1956
1957	large = kinds.default_float_kind.MAX
1958	min= large
1959	max=-large
1960	for vertex in self.userVertices:
1961	if len(vertex.attributes) > iatt:
1962	if vertex.attributes[iatt] < min:
1963	min = vertex.attributes[iatt]
1964	if vertex.attributes[iatt] > max:
1965	max = vertex.attributes[iatt]
1966	for vertex in self.meshVertices:
1967	if len(vertex.attributes) > iatt:
1968	if vertex.attributes[iatt] < min:
1969	min = vertex.attributes[iatt]
1970	if vertex.attributes[iatt] > max:
1971	max = vertex.attributes[iatt]
1972	return [min, max]
1973
1974	def scaleoffset(self, WIDTH, HEIGHT):
1975	"""
1976	Returns a list denoting the scale and offset terms that need to be
1977	applied when converting mesh co-ordinates onto grid co-ordinates
1978	Structure: [scale, xoffset, yoffset]
1979	"""
1980	OFFSET = 0.05*min([WIDTH, HEIGHT])
1981	[xmin, ymin, xmax, ymax] = self.boxsize()
1982	SCALE = min([0.9WIDTH, 0.9HEIGHT])/max([xmax-xmin, ymax-ymin])
1983
1984	if SCALE*xmin < OFFSET:
1985	xoffset = abs(SCALE*xmin) + OFFSET
1986	if SCALE*xmax > WIDTH - OFFSET:
1987	xoffset= -(SCALE*xmax - WIDTH + OFFSET)
1988	if SCALE*ymin < OFFSET:
1989	b = abs(SCALE*ymin)+OFFSET
1990	if SCALE*ymax > HEIGHT-OFFSET:
1991	b = -(SCALE*ymax - HEIGHT + OFFSET)
1992	yoffset = HEIGHT - b
1993	return [SCALE, xoffset, yoffset]
1994
1995
1996	def exportASCIIobj(self,ofile):
1997	"""
1998	export a file, ofile, with the format
1999	lines: v <x> <y> <first attribute>
2000	f <vertex #> <vertex #> <vertex #> (of the triangles)
2001	"""
2002	fd = open(ofile,'w')
2003	self.writeASCIIobj(fd)
2004	fd.close()
2005
2006
2007	def writeASCIIobj(self,fd):
2008	fd.write(" # Triangulation as an obj file\n")
2009	numVert = str(len(self.meshVertices))
2010
2011	index1 = 1
2012	for vert in self.meshVertices:
2013	vert.index1 = index1
2014	index1 += 1
2015
2016	fd.write("v "
2017	+ str(vert.x) + " "
2018	+ str(vert.y) + " "
2019	+ str(vert.attributes[0]) + "\n")
2020
2021	for tri in self.meshTriangles:
2022	fd.write("f "
2023	+ str(tri.vertices[0].index1) + " "
2024	+ str(tri.vertices[1].index1) + " "
2025	+ str(tri.vertices[2].index1) + "\n")
2026
2027	def exportASCIIsegmentoutlinefile(self,ofile):
2028	"""Write the boundary user mesh info, eg
2029	vertices that are connected to segments,
2030	segments
2031	"""
2032
2033	verts = {}
2034	for seg in self.getUserSegments():
2035	verts[seg.vertices[0]] = seg.vertices[0]
2036	verts[seg.vertices[1]] = seg.vertices[1]
2037	meshDict = self.Mesh2IOOutlineDict(userVertices=verts.values())
2038	export_mesh_file(ofile,meshDict)
2039
2040	# exportASCIImeshfile - this function is used
2041	def export_mesh_file(self,ofile):
2042	"""
2043	export a file, ofile, with the format
2044	"""
2045
2046	dict = self.Mesh2IODict()
2047	export_mesh_file(ofile,dict)
2048
2049	# FIXME(DSG-DSG):Break this into two functions.
2050	#One for the outline points.
2051	#One for the mesh points.
2052	# Note: this function is not in the gui
2053	def exportPointsFile(self,ofile):
2054	"""
2055	export a points file, ofile.
2056
2057	"""
2058
2059	mesh_dict = self.Mesh2IODict()
2060	#point_dict = {}
2061	#point_dict['attributelist'] = {} #this will need to be expanded..
2062	# if attributes are brought back in.
2063	#point_dict['geo_reference'] = self.geo_reference
2064	if mesh_dict['vertices'] == []:
2065	#point_dict['pointlist'] = mesh_dict['points']
2066	geo = Geospatial_data(mesh_dict['points'],
2067	geo_reference=self.geo_reference)
2068	else:
2069	#point_dict['pointlist'] = mesh_dict['vertices']
2070	geo = Geospatial_data(mesh_dict['vertices'],
2071	geo_reference=self.geo_reference)
2072
2073	geo.export_points_file(ofile, absolute=True)
2074
2075
2076
2077	def import_ungenerate_file(self,ofile, tag=None):
2078	"""
2079	Imports an ungenerate file, from arcGIS into mesh.
2080
2081	ofile is the name of the ungenerated file.
2082	Tag is a string name to be taggged on each segment.
2083
2084	WARNING values are assumed to be absolute.
2085	geo-refs are not taken into account..
2086	"""
2087
2088	dict = importUngenerateFile(ofile)
2089	default_tag = Segment.get_default_tag()
2090	if tag is not None:
2091	Segment.set_default_tag(str(tag))
2092	self.addVertsSegs(dict)
2093	Segment.set_default_tag(default_tag)
2094
2095	# change the tag back to it's default
2096
2097
2098	########### IO CONVERTERS ##################
2099	"""
2100	The dict fromat for IO with .tsh files is;
2101	(the triangulation)
2102	vertices: [[x1,y1],[x2,y2],...] (lists of doubles)
2103	vertex_attributes: [[a11,a12,...],[a21,a22],...] (lists of doubles)
2104	vertex_attribute_titles:[A1Title, A2Title ...] (A list of strings)
2105	segments: [[v1,v2],[v3,v4],...] (lists of integers)
2106	segment_tags : [tag,tag,...] list of strings
2107	triangles : [(v1,v2,v3), (v4,v5,v6),....] lists of points
2108	triangle tags: [s1,s2,...] A list of strings
2109	triangle neighbors: [[t1,t2,t3], [t4,t5,t6],..] lists of triangles
2110
2111	(the outline)
2112	points: [[x1,y1],[x2,y2],...] (lists of doubles)
2113	point_attributes: [[a11,a12,...],[a21,a22],...] (lists of doubles)
2114	outline_segments: [[point1,point2],[p3,p4],...] (lists of integers)
2115	outline_segment_tags : [tag1,tag2,...] list of strings
2116	holes : [[x1,y1],...](List of doubles, one inside each hole region)
2117	regions : [ [x1,y1],...] (List of 4 doubles)
2118	region_tags : [tag1,tag2,...] (list of strings)
2119	region_max_areas: [ma1,ma2,...] (A list of doubles)
2120	{Convension: A -ve max area means no max area}
2121
2122	"""
2123
2124
2125
2126	def Mesh2IODict(self):
2127	"""
2128	Convert the triangulation and outline info of a mesh to a dictionary
2129	structure
2130	"""
2131	dict = self.Mesh2IOTriangulationDict()
2132	dict_mesh = self.Mesh2IOOutlineDict()
2133	for element in dict_mesh.keys():
2134	dict[element] = dict_mesh[element]
2135
2136	# add the geo reference
2137	dict['geo_reference'] = self.geo_reference
2138	return dict
2139
2140	def Mesh2IOTriangulationDict(self):
2141	"""
2142	Convert the Mesh to a dictionary of lists describing the
2143	triangulation variables;
2144
2145	Used to produce .tsh file
2146	"""
2147
2148	meshDict = {}
2149	vertices=[]
2150	vertex_attributes=[]
2151
2152	self.maxVertexIndex=0
2153	for vertex in self.meshVertices:
2154	vertex.index = self.maxVertexIndex
2155	vertices.append([vertex.x,vertex.y])
2156	vertex_attributes.append(vertex.attributes)
2157	self.maxVertexIndex += 1
2158
2159	meshDict['vertices'] = vertices
2160	meshDict['vertex_attributes'] = vertex_attributes
2161	meshDict['vertex_attribute_titles'] = self.attributeTitles
2162	#segments
2163	segments=[]
2164	segment_tags=[]
2165	for seg in self.meshSegments:
2166	segments.append([seg.vertices[0].index,seg.vertices[1].index])
2167	segment_tags.append(seg.tag)
2168	meshDict['segments'] =segments
2169	meshDict['segment_tags'] =segment_tags
2170
2171	# Make sure that the indexation is correct
2172	index = 0
2173	for tri in self.meshTriangles:
2174	tri.index = index
2175	index += 1
2176
2177	triangles = []
2178	triangle_tags = []
2179	triangle_neighbors = []
2180	for tri in self.meshTriangles:
2181	triangles.append([tri.vertices[0].index,
2182	tri.vertices[1].index,
2183	tri.vertices[2].index])
2184	triangle_tags.append(tri.attribute)
2185	neighborlist = [-1,-1,-1]
2186	for neighbor,index in map(None,tri.neighbors,
2187	range(len(tri.neighbors))):
2188	if neighbor:
2189	neighborlist[index] = neighbor.index
2190	triangle_neighbors.append(neighborlist)
2191
2192	meshDict['triangles'] = triangles
2193	meshDict['triangle_tags'] = triangle_tags
2194	meshDict['triangle_neighbors'] = triangle_neighbors
2195
2196	#print "mesh.Mesh2IOTriangulationDict))"
2197	#print meshDict
2198	#print "mesh.Mesh2IOTriangulationDict))"
2199
2200	return meshDict
2201
2202
2203	def Mesh2IOOutlineDict(self, userVertices=None,
2204	userSegments=None,
2205	holes=None,
2206	regions=None):
2207	"""
2208	Convert the mesh outline to a dictionary of the lists needed for the
2209	triang module;
2210
2211	Note, this adds an index attribute to the user Vertex objects.
2212
2213	Used to produce .tsh file and output to triangle
2214	"""
2215
2216	if userVertices is None:
2217	userVertices = self.getUserVertices()
2218	if userSegments is None:
2219	userSegments = self.getUserSegments()
2220	if holes is None:
2221	holes = self.getHoles()
2222	if regions is None:
2223	regions = self.getRegions()
2224
2225	meshDict = {}
2226	#print "userVertices",userVertices
2227	#print "userSegments",userSegments
2228	pointlist=[]
2229	pointattributelist=[]
2230	index = 0
2231	for vertex in userVertices:
2232	vertex.index = index
2233	pointlist.append([vertex.x,vertex.y])
2234	pointattributelist.append(vertex.attributes)
2235
2236	index += 1
2237	meshDict['points'] = pointlist
2238	meshDict['point_attributes'] = pointattributelist
2239
2240	segmentlist=[]
2241	segmenttaglist=[]
2242	for seg in userSegments:
2243	segmentlist.append([seg.vertices[0].index,seg.vertices[1].index])
2244	segmenttaglist.append(seg.tag)
2245	meshDict['outline_segments'] =segmentlist
2246	meshDict['outline_segment_tags'] =segmenttaglist
2247
2248	holelist=[]
2249	for hole in holes:
2250	holelist.append([hole.x,hole.y])
2251	meshDict['holes'] = holelist
2252
2253	regionlist=[]
2254	regiontaglist = []
2255	regionmaxarealist = []
2256	for region in regions:
2257	regionlist.append([region.x,region.y])
2258	regiontaglist.append(region.getTag())
2259
2260	if (region.getMaxArea() != None):
2261	regionmaxarealist.append(region.getMaxArea())
2262	else:
2263	regionmaxarealist.append(NOMAXAREA)
2264	meshDict['regions'] = regionlist
2265	meshDict['region_tags'] = regiontaglist
2266	meshDict['region_max_areas'] = regionmaxarealist
2267	#print "(("
2268	#print meshDict
2269	#print meshDict['regionlist']
2270	#print "(("
2271	return meshDict
2272
2273	def IOTriangulation2Mesh(self, genDict):
2274	"""
2275	Set the mesh attributes given an tsh IO dictionary
2276	"""
2277	#Clear the current generated mesh values
2278	self.meshTriangles=[]
2279	self.attributeTitles=[]
2280	self.meshSegments=[]
2281	self.meshVertices=[]
2282
2283	#print "mesh.setTriangulation@#@#@#"
2284	#print genDict
2285	#print "@#@#@#"
2286
2287	self.maxVertexIndex = 0
2288	for point in genDict['vertices']:
2289	v=Vertex(point[0], point[1])
2290	v.index = self.maxVertexIndex
2291	self.maxVertexIndex +=1
2292	self.meshVertices.append(v)
2293
2294	self.attributeTitles = genDict['vertex_attribute_titles']
2295
2296	index = 0
2297	for seg,tag in map(None,genDict['segments'],genDict['segment_tags']):
2298	segObject = Segment( self.meshVertices[int(seg[0])],
2299	self.meshVertices[int(seg[1])], tag = tag )
2300	segObject.index = index
2301	index +=1
2302	self.meshSegments.append(segObject)
2303
2304	index = 0
2305	for triangle in genDict['triangles']:
2306	tObject =Triangle( self.meshVertices[int(triangle[0])],
2307	self.meshVertices[int(triangle[1])],
2308	self.meshVertices[int(triangle[2])] )
2309	tObject.index = index
2310	index +=1
2311	self.meshTriangles.append(tObject)
2312
2313	index = 0
2314	for att in genDict['triangle_tags']:
2315	if att == []:
2316	self.meshTriangles[index].setAttribute("")
2317	else:
2318	self.meshTriangles[index].setAttribute(att)
2319	index += 1
2320
2321	index = 0
2322	for att in genDict['vertex_attributes']:
2323	if att == None:
2324	self.meshVertices[index].setAttributes([])
2325	else:
2326	self.meshVertices[index].setAttributes(att)
2327	index += 1
2328
2329	index = 0
2330	for triangle in genDict['triangle_neighbors']:
2331	# Build a list of triangle object neighbors
2332	ObjectNeighbor = []
2333	for neighbor in triangle:
2334	if ( neighbor != -1):
2335	ObjectNeighbor.append(self.meshTriangles[int(neighbor)])
2336	else:
2337	ObjectNeighbor.append(None)
2338	self.meshTriangles[index].setNeighbors(ObjectNeighbor[0],
2339	ObjectNeighbor[1],
2340	ObjectNeighbor[2])
2341	index += 1
2342
2343
2344	def IOOutline2Mesh(self, genDict):
2345	"""
2346	Set the outline (user Mesh attributes) given a IO tsh dictionary
2347
2348	mesh is an instance of a mesh object
2349	"""
2350	#Clear the current user mesh values
2351	self.clearUserSegments()
2352	self.userVertices=[]
2353	self.Holes=[]
2354	self.Regions=[]
2355
2356	#print "mesh.IOOutline2Mesh@#@#@#"
2357	#print "genDict",genDict
2358	#print "@#@#@#"
2359
2360	#index = 0
2361	for point in genDict['points']:
2362	v=Vertex(point[0], point[1])
2363	#v.index = index
2364	#index +=1
2365	self.userVertices.append(v)
2366
2367	#index = 0
2368	for seg,tag in map(None,genDict['outline_segments'],
2369	genDict['outline_segment_tags']):
2370
2371	segObject = Segment( self.userVertices[int(seg[0])],
2372	self.userVertices[int(seg[1])], tag = tag )
2373	#segObject.index = index
2374	#index +=1
2375	self.userSegments.append(segObject)
2376
2377	# Remove the loading of attribute info.
2378	# Have attribute info added using least_squares in pyvolution
2379	# index = 0
2380	# for att in genDict['point_attributes']:
2381	# if att == None:
2382	# self.userVertices[index].setAttributes([])
2383	# else:
2384	# self.userVertices[index].setAttributes(att)
2385	# index += 1
2386
2387	#index = 0
2388	for point in genDict['holes']:
2389	h=Hole(point[0], point[1])
2390	#h.index = index
2391	#index +=1
2392	self.holes.append(h)
2393
2394	#index = 0
2395	for reg,att,maxArea in map(None,
2396	genDict['regions'],
2397	genDict['region_tags'],
2398	genDict['region_max_areas']):
2399	if maxArea > 0: # maybe I should ref NOMAXAREA? Prob' not though
2400	Object = Region( reg[0],
2401	reg[1],
2402	tag = att,
2403	maxArea = maxArea)
2404	else:
2405	Object = Region( reg[0],
2406	reg[1],
2407	tag = att)
2408
2409	#Object.index = index
2410	#index +=1
2411	self.regions.append(Object)
2412
2413	############################################
2414
2415
2416	def refineSet(self,setName):
2417	Triangles = self.sets[self.setID[setName]]
2418	Refine(self,Triangles)
2419
2420	def selectAllTriangles(self):
2421	A=[]
2422	A.extend(self.meshTriangles)
2423	if not('All' in self.setID.keys()):
2424	self.setID['All']=len(self.sets)
2425	self.sets.append(A)
2426	else:
2427	self.sets[self.setID['All']]=A
2428	return 'All'
2429	# and objectIDs
2430
2431
2432	def clearSelection(self):
2433	A = []
2434	if not('None' in self.setID.keys()):
2435	self.setID['None']=len(self.sets)
2436	self.sets.append(A)
2437	return 'None'
2438
2439	def drawSet(self,canvas,setName,SCALE,colour=SET_COLOUR):
2440	#FIXME Draws over previous triangles - may bloat canvas
2441	Triangles = self.sets[self.setID[setName]]
2442	for triangle in Triangles:
2443	triangle.draw(canvas,1,
2444	scale = SCALE,
2445	colour = colour)
2446
2447	def undrawSet(self,canvas,setName,SCALE,colour='green'):
2448	#FIXME Draws over previous lines - may bloat canvas
2449	Triangles = self.sets[self.setID[setName]]
2450	for triangle in Triangles:
2451	triangle.draw(canvas,1,
2452	scale = SCALE,
2453	colour = colour)
2454
2455	def weed(self,Vertices,Segments):
2456	#Depreciated
2457	#weed out existing duplicates
2458	print 'len(self.getUserSegments())'
2459	print len(self.getUserSegments())
2460	print 'len(self.getUserVertices())'
2461	print len(self.getUserVertices())
2462
2463	point_keys = {}
2464	for vertex in Vertices:
2465	point = (vertex.x,vertex.y)
2466	point_keys[point]=vertex
2467	#inlined would looks very ugly
2468
2469	line_keys = {}
2470	for segment in Segments:
2471	vertex1 = segment.vertices[0]
2472	vertex2 = segment.vertices[1]
2473	point1 = (vertex1.x,vertex1.y)
2474	point2 = (vertex2.x,vertex2.y)
2475	segment.vertices[0]=point_keys[point1]
2476	segment.vertices[1]=point_keys[point2]
2477	vertex1 = segment.vertices[0]
2478	vertex2 = segment.vertices[1]
2479	point1 = (vertex1.x,vertex1.y)
2480	point2 = (vertex2.x,vertex2.y)
2481	line1 = (point1,point2)
2482	line2 = (point2,point1)
2483	if not (line_keys.has_key(line1) \
2484	or line_keys.has_key(line2)):
2485	line_keys[line1]=segment
2486	Vertices=point_keys.values()
2487	Segments=line_keys.values()
2488	return Vertices,Segments
2489
2490	def segs_to_dict(self,segments):
2491	dict={}
2492	for segment in segments:
2493	vertex1 = segment.vertices[0]
2494	vertex2 = segment.vertices[1]
2495	point1 = (vertex1.x,vertex1.y)
2496	point2 = (vertex2.x,vertex2.y)
2497	line = (point1,point2)
2498	dict[line]=segment
2499	return dict
2500
2501	def seg2line(self,s):
2502	return ((s.vertices[0].x,s.vertices[0].y,)\
2503	(s.vertices[1].x,s.vertices[1].y))
2504
2505	def line2seg(self,line,tag=None):
2506	point0 = self.point2ver(line[0])
2507	point1 = self.point2ver(line[1])
2508	return Segment(point0,point1,tag=tag)
2509
2510	def ver2point(self,vertex):
2511	return (vertex.x,vertex.y)
2512
2513	def point2ver(self,point):
2514	return Vertex(point[0],point[1])
2515
2516	def smooth_polySet(self,min_radius=0.05):
2517	#for all pairs of connecting segments:
2518	# propose a new segment that replaces the 2
2519
2520	# If the difference between the new segment
2521	# and the old lines is small: replace the
2522	# old lines.
2523
2524	seg2line = self.seg2line
2525	ver2point= self.ver2point
2526	line2seg = self.line2seg
2527	point2ver= self.point2ver
2528
2529	#create dictionaries of lines -> segments
2530	userSegments = self.segs_to_dict(self.userSegments)
2531	alphaSegments = self.segs_to_dict(self.alphaUserSegments)
2532
2533	#lump user and alpha segments
2534	for key in alphaSegments.keys():
2535	userSegments[key]=alphaSegments[key]
2536
2537	#point_keys = tuple -> vertex
2538	#userVertices = vertex -> [line,line] - lines from that node
2539	point_keys = {}
2540	userVertices={}
2541	for vertex in self.getUserVertices():
2542	point = ver2point(vertex)
2543	if not point_keys.has_key(point):
2544	point_keys[point]=vertex
2545	userVertices[vertex]=[]
2546	for key in userSegments.keys():
2547	line = key
2548	point_0 = key[0]
2549	point_1 = key[1]
2550	userVertices[point_keys[point_0]].append(line)
2551	userVertices[point_keys[point_1]].append(line)
2552
2553	for point in point_keys.keys():
2554	try:
2555	#removed keys can cause keyerrors
2556	vertex = point_keys[point]
2557	lines = userVertices[vertex]
2558
2559	#if there are 2 lines on the node
2560	if len(lines)==2:
2561	line_0 = lines[0]
2562	line_1 = lines[1]
2563
2564	#if the tags are the the same on the 2 lines
2565	if userSegments[line_0].tag == userSegments[line_1].tag:
2566	tag = userSegments[line_0].tag
2567
2568	#point_a is one of the next nodes, point_b is the other
2569	if point==line_0[0]:
2570	point_a = line_0[1]
2571	if point==line_0[1]:
2572	point_a = line_0[0]
2573	if point==line_1[0]:
2574	point_b = line_1[1]
2575	if point==line_1[1]:
2576	point_b = line_1[0]
2577
2578
2579	#line_2 is proposed
2580	line_2 = (point_a,point_b)
2581
2582	#calculate the area of the triangle between
2583	#the two existing segments and the proposed
2584	#new segment
2585	ax = point_a[0]
2586	ay = point_a[1]
2587	bx = point_b[0]
2588	by = point_b[1]
2589	cx = point[0]
2590	cy = point[1]
2591	area=abs((bxay-axby)+(cxby-bxcy)+(axcy-cxay))/2
2592
2593	#calculate the perimeter
2594	len_a = ((cx-bx)2+(cy-by)2)**0.5
2595	len_b = ((ax-cx)2+(ay-cy)2)**0.5
2596	len_c = ((bx-ax)2+(by-ay)2)**0.5
2597	perimeter = len_a+len_b+len_c
2598
2599	#calculate the radius
2600	r = area/(2*perimeter)
2601
2602	#if the radius is small: then replace the existing
2603	#segments with the new one
2604	if r < min_radius:
2605	if len_c < min_radius: append = False
2606	else: append = True
2607	#if the new seg is also time, don't add it
2608	if append:
2609	segment = self.line2seg(line_2,tag=tag)
2610
2611	list_a=userVertices[point_keys[point_a]]
2612	list_b=userVertices[point_keys[point_b]]
2613
2614	if line_0 in list_a:
2615	list_a.remove(line_0)
2616	else:
2617	list_a.remove(line_1)
2618
2619	if line_0 in list_b:
2620	list_b.remove(line_0)
2621	else:
2622	list_b.remove(line_1)
2623
2624	if append:
2625	list_a.append(line_2)
2626	list_b.append(line_2)
2627	else:
2628	if len(list_a)==0:
2629	userVertices.pop(point_keys[point_a])
2630	point_keys.pop(point_a)
2631	if len(list_b)==0:
2632	userVertices.pop(point_keys[point_b])
2633	point_keys.pop(point_b)
2634
2635	userVertices.pop(point_keys[point])
2636	point_keys.pop(point)
2637	userSegments.pop(line_0)
2638	userSegments.pop(line_1)
2639
2640	if append:
2641	userSegments[line_2]=segment
2642	except:
2643	pass
2644
2645	#self.userVerticies = userVertices.keys()
2646	#self.userSegments = []
2647	#for key in userSegments.keys():
2648	# self.userSegments.append(userSegments[key])
2649	#self.alphaUserSegments = []
2650
2651	self.userVerticies = []
2652	self.userSegments = []
2653	self.alphaUserSegments = []
2654
2655	return userVertices,userSegments,alphaSegments
2656
2657	def triangles_to_polySet(self,setName):
2658	#self.smooth_polySet()
2659
2660	seg2line = self.seg2line
2661	ver2point= self.ver2point
2662	line2seg = self.line2seg
2663	point2ver= self.point2ver
2664
2665	from Numeric import array,allclose
2666	#turn the triangles into a set
2667	Triangles = self.sets[self.setID[setName]]
2668	Triangles_dict = {}
2669	for triangle in Triangles:
2670	Triangles_dict[triangle]=None
2671
2672
2673	#create a dict of points to vertexes (tuple -> object)
2674	#also create a set of vertexes (object -> True)
2675	point_keys = {}
2676	userVertices={}
2677	for vertex in self.getUserVertices():
2678	point = ver2point(vertex)
2679	if not point_keys.has_key(point):
2680	point_keys[point]=vertex
2681	userVertices[vertex]=True
2682
2683	#create a dict of lines to segments (tuple -> object)
2684	userSegments = self.segs_to_dict(self.userSegments)
2685	#append the userlines in an affine linespace
2686	affine_lines = Affine_Linespace()
2687	for line in userSegments.keys():
2688	affine_lines.append(line)
2689	alphaSegments = self.segs_to_dict(self.alphaUserSegments)
2690	for line in alphaSegments.keys():
2691	affine_lines.append(line)
2692
2693	for triangle in Triangles:
2694	for i in (0,1,2):
2695	#for every triangles neighbour:
2696	if not Triangles_dict.has_key(triangle.neighbors[i]):
2697	#if the neighbour is not in the set:
2698	a = triangle.vertices[i-1]
2699	b = triangle.vertices[i-2]
2700	#Get possible matches:
2701	point_a = ver2point(a)
2702	point_b = ver2point(b)
2703	midpoint = ((a.x+b.x)/2,(a.y+b.y)/2)
2704	line = (point_a,point_b)
2705	tag = None
2706
2707
2708	#this bit checks for matching lines
2709	possible_lines = affine_lines[line]
2710	possible_lines = unique(possible_lines)
2711	found = 0
2712	for user_line in possible_lines:
2713	if self.point_on_line(midpoint,user_line):
2714	found+=1
2715	assert found<2
2716	if userSegments.has_key(user_line):
2717	parent_segment = userSegments.pop(user_line)
2718	if alphaSegments.has_key(user_line):
2719	parent_segment = alphaSegments.pop(user_line)
2720	tag = parent_segment.tag
2721	offspring = [line]
2722	offspring.extend(self.subtract_line(user_line,
2723	line))
2724	affine_lines.remove(user_line)
2725	for newline in offspring:
2726	line_vertices = []
2727	for point in newline:
2728	if point_keys.has_key(point):
2729	vert = point_keys[point]
2730	else:
2731	vert = Vertex(point[0],point[1])
2732	userVertices[vert]=True
2733	point_keys[point]=vert
2734	line_vertices.append(vert)
2735	segment = Segment(line_vertices[0],
2736	line_vertices[1],tag)
2737	userSegments[newline]=segment
2738	affine_lines.append(newline)
2739	#break
2740	assert found<2
2741
2742
2743
2744	#if no matching lines
2745	if not found:
2746	line_vertices = []
2747	for point in line:
2748	if point_keys.has_key(point):
2749	vert = point_keys[point]
2750	else:
2751	vert = Vertex(point[0],point[1])
2752	userVertices[vert]=True
2753	point_keys[point]=vert
2754	line_vertices.append(vert)
2755	segment = Segment(line_vertices[0],
2756	line_vertices[1],tag)
2757	userSegments[line]=segment
2758	affine_lines.append(line)
2759
2760	self.userVerticies = []
2761	self.userSegments = []
2762	self.alphaUserSegments = []
2763
2764	return userVertices,userSegments,alphaSegments
2765
2766	def subtract_line(self,parent,child):
2767	#Subtracts child from parent
2768	#Requires that the child is a
2769	#subline of parent to work.
2770
2771	from Numeric import allclose,dot,array
2772	A= parent[0]
2773	B= parent[1]
2774	a = child[0]
2775	b = child[1]
2776
2777	A_array = array(parent[0])
2778	B_array = array(parent[1])
2779	a_array = array(child[0])
2780	b_array = array(child[1])
2781
2782	assert not A == B
2783	assert not a == b
2784
2785	answer = []
2786
2787	#if the new line does not share a
2788	#vertex with the old one
2789	if not (allclose(A_array,a_array)\
2790	or allclose(B_array,b_array)\
2791	or allclose(A_array,b_array)\
2792	or allclose(a_array,B_array)):
2793	if dot(A_array-a_array,A_array-a_array) \
2794	< dot(A_array-b_array,A_array-b_array):
2795	sibling1 = (A,a)
2796	sibling2 = (B,b)
2797	return [sibling1,sibling2]
2798	else:
2799	sibling1 = (A,b)
2800	sibling2 = (B,a)
2801	return [sibling1,sibling2]
2802
2803	elif allclose(A_array,a_array):
2804	if allclose(B_array,b_array):
2805	return []
2806	else:
2807	sibling = (b,B)
2808	return [sibling]
2809	elif allclose(B_array,b_array):
2810	sibling = (a,A)
2811	return [sibling]
2812
2813	elif allclose(A_array,b_array):
2814	if allclose(B,a):
2815	return []
2816	else:
2817	sibling = (a,B)
2818	return [sibling]
2819	elif allclose(a_array,B_array):
2820	sibling = (b,A)
2821	return [sibling]
2822
2823	def point_on_line(self,point,line):
2824	#returns true within a tolerance of 3 degrees
2825	x=point[0]
2826	y=point[1]
2827	x0=line[0][0]
2828	x1=line[1][0]
2829	y0=line[0][1]
2830	y1=line[1][1]
2831	from Numeric import array, dot, allclose
2832	from math import sqrt
2833	tol = 3. #DEGREES
2834	tol = tol*3.1415/180
2835
2836	a = array([x - x0, y - y0])
2837	a_normal = array([a[1], -a[0]])
2838	len_a_normal = sqrt(sum(a_normal**2))
2839
2840	b = array([x1 - x0, y1 - y0])
2841	len_b = sqrt(sum(b**2))
2842
2843	if abs(dot(a_normal, b)/(len_b*len_a_normal))< tol:
2844	#Point is somewhere on the infinite extension of the line
2845
2846	len_a = sqrt(sum(a**2))
2847	if dot(a, b) >= 0 and len_a <= len_b:
2848	return True
2849	else:
2850	return False
2851	else:
2852	return False
2853
2854	def line_length(self,line):
2855	x0=line[0][0]
2856	x1=line[1][0]
2857	y0=line[0][1]
2858	y1=line[1][1]
2859	return ((x1-x0)2-(y1-y0)2)**0.5
2860
2861	def threshold(self,setName,min=None,max=None,attribute_name='elevation'):
2862	"""
2863	threshold using d
2864	"""
2865	triangles = self.sets[self.setID[setName]]
2866	A = []
2867
2868	if attribute_name in self.attributeTitles:
2869	i = self.attributeTitles.index(attribute_name)
2870	else: i = -1#no attribute
2871	if not max == None:
2872	for t in triangles:
2873	if (min<self.av_att(t,i)<max):
2874	A.append(t)
2875	else:
2876	for t in triangles:
2877	if (min<self.av_att(t,i)):
2878	A.append(t)
2879	self.sets[self.setID[setName]] = A
2880
2881	def general_threshold(self,setName,min=None,max=None\
2882	,attribute_name = 'elevation',function=None):
2883	"""
2884	Thresholds the triangles
2885	"""
2886	from visual.graph import arange,ghistogram,color as colour
2887	triangles = self.sets[self.setID[setName]]
2888	A = []
2889	data=[]
2890	#data is for the graph
2891
2892	if attribute_name in self.attributeTitles:
2893	i = self.attributeTitles.index(attribute_name)
2894	else: i = -1
2895	if not max == None:
2896	for t in triangles:
2897	value=function(t,i)
2898	if (min<value<max):
2899	A.append(t)
2900	data.append(value)
2901	else:
2902	for t in triangles:
2903	value=function(t,i)
2904	if (min<value):
2905	A.append(t)
2906	data.append(value)
2907	self.sets[self.setID[setName]] = A
2908
2909	if self.visualise_graph:
2910	if len(data)>0:
2911	max=data[0]
2912	min=data[0]
2913	for value in data:
2914	if value > max:
2915	max = value
2916	if value < min:
2917	min = value
2918
2919	inc = (max-min)/100
2920
2921	histogram = ghistogram(bins=arange(min,max,inc),\
2922	color = colour.red)
2923	histogram.plot(data=data)
2924
2925	def av_att(self,triangle,i):
2926	if i==-1: return 1
2927	else:
2928	#evaluates the average attribute of the vertices of a triangle.
2929	V = triangle.getVertices()
2930	a0 = (V[0].attributes[i])
2931	a1 = (V[1].attributes[i])
2932	a2 = (V[2].attributes[i])
2933	return (a0+a1+a2)/3
2934
2935	def Courant_ratio(self,triangle,index):
2936	"""
2937	Uses the courant threshold
2938	"""
2939	e = self.av_att(triangle,index)
2940	A = triangle.calcArea()
2941	P = triangle.calcP()
2942	r = A/(2*P)
2943	e = max(0.1,abs(e))
2944	return r/e**0.5
2945
2946	def Gradient(self,triangle,index):
2947	V = triangle.vertices
2948	x0, y0, x1, y1, x2, y2, q0, q1, q2 = V[0].x,V[0].y,V[1].x,V[1].y,V[2].x,V[2].y,V[0].attributes[index],V[1].attributes[index],V[2].attributes[index]
2949	grad_x,grad_y = gradient(x0, y0, x1, y1, x2, y2, q0, q1, q2)
2950	if ((grad_x2)+(grad_y2))**(0.5)<0:
2951	print ((grad_x2)+(grad_y2))**(0.5)
2952	return ((grad_x2)+(grad_y2))**(0.5)
2953
2954
2955	def append_triangle(self,triangle):
2956	self.meshTriangles.append(triangle)
2957
2958	def replace_triangle(self,triangle,replacement):
2959	i = self.meshTriangles.index(triangle)
2960	self.meshTriangles[i]=replacement
2961	assert replacement in self.meshTriangles
2962
2963	def importUngenerateFile(ofile):
2964	"""
2965	import a file, ofile, with the format
2966	[poly]
2967	poly format:
2968	First line: <# of vertices> <x centroid> <y centroid>
2969	Following lines: <x> <y>
2970	last line: "END"
2971
2972	Note: These are clockwise.
2973	"""
2974	fd = open(ofile,'r')
2975	Dict = readUngenerateFile(fd)
2976	fd.close()
2977	return Dict
2978
2979	def readUngenerateFile(fd):
2980	"""
2981	import a file, ofile, with the format
2982	[poly]
2983	poly format:
2984	First line: <# of polynomial> <x centroid> <y centroid>
2985	Following lines: <x> <y>
2986	last line: "END"
2987	"""
2988	END_DELIMITER = 'END'
2989
2990	points = []
2991	segments = []
2992
2993	isEnd = False
2994	line = fd.readline() #not used <# of polynomial> <x> <y>
2995	while not isEnd:
2996	line = fd.readline()
2997	fragments = line.split()
2998	vert = [float(fragments.pop(0)),float(fragments.pop(0))]
2999	points.append(vert)
3000	PreviousVertIndex = len(points)-1
3001	firstVertIndex = PreviousVertIndex
3002
3003	line = fd.readline() #Read the next line
3004	while not line.startswith(END_DELIMITER):
3005	#print "line >" + line + "<"
3006	fragments = line.split()
3007	vert = [float(fragments.pop(0)),float(fragments.pop(0))]
3008	points.append(vert)
3009	thisVertIndex = len(points)-1
3010	segment = [PreviousVertIndex,thisVertIndex]
3011	segments.append(segment)
3012	PreviousVertIndex = thisVertIndex
3013	line = fd.readline() #Read the next line
3014	i =+ 1
3015	# If the last and first segments are the same,
3016	# Remove the last segment and the last vertex
3017	# then add a segment from the second last vert to the 1st vert
3018	thisVertIndex = len(points)-1
3019	firstVert = points[firstVertIndex]
3020	thisVert = points[thisVertIndex]
3021	#print "firstVert",firstVert
3022	#print "thisVert",thisVert
3023	if (firstVert[0] == thisVert[0] and firstVert[1] == thisVert[1]):
3024	points.pop()
3025	segments.pop()
3026	thisVertIndex = len(points)-1
3027	segments.append([thisVertIndex, firstVertIndex])
3028
3029	line = fd.readline() # read <# of polynomial> <x> <y> OR END
3030	#print "line >>" + line + "<<"
3031	if line.startswith(END_DELIMITER):
3032	isEnd = True
3033
3034	#print "points", points
3035	#print "segments", segments
3036	ungenerated_dict = {}
3037	ungenerated_dict['points'] = points
3038	ungenerated_dict['segments'] = segments
3039	return ungenerated_dict
3040
3041	def importMeshFromFile(ofile):
3042	"""returns a mesh object, made from a points file or .tsh/.msh file
3043	Often raises IOError,RuntimeError
3044	"""
3045	newmesh = None
3046	if (ofile[-4:]== ".pts" or ofile[-4:]== ".txt" or \
3047	ofile[-4:]== ".csv"):
3048	geospatial = Geospatial_data(ofile)
3049	dict = {}
3050	dict['points'] = geospatial.get_data_points(absolute=False)
3051	dict['outline_segments'] = []
3052	dict['outline_segment_tags'] = []
3053	dict['regions'] = []
3054	dict['region_tags'] = []
3055	dict['region_max_areas'] = []
3056	dict['holes'] = []
3057	newmesh= Mesh(geo_reference = geospatial.geo_reference)
3058	newmesh.IOOutline2Mesh(dict)
3059	counter = newmesh.removeDuplicatedUserVertices()
3060	if (counter >0):
3061	print "%i duplicate vertices removed from dataset" % (counter)
3062	elif (ofile[-4:]== ".tsh" or ofile[-4:]== ".msh"):
3063	dict = import_mesh_file(ofile)
3064	#print "********"
3065	#print "zq mesh.dict",dict
3066	#print "********"
3067	newmesh= Mesh()
3068	newmesh.IOOutline2Mesh(dict)
3069	newmesh.IOTriangulation2Mesh(dict)
3070	else:
3071	raise RuntimeError
3072
3073	if dict.has_key('geo_reference') and not dict['geo_reference'] == None:
3074	newmesh.geo_reference = dict['geo_reference']
3075	return newmesh
3076
3077	def loadPickle(currentName):
3078	fd = open(currentName)
3079	mesh = pickle.load(fd)
3080	fd.close()
3081	return mesh
3082
3083	def square_outline(side_length = 1,up = "top", left = "left", right = "right",
3084	down = "bottom", regions = False):
3085
3086	a = Vertex (0,0)
3087	b = Vertex (0,side_length)
3088	c = Vertex (side_length,0)
3089	d = Vertex (side_length,side_length)
3090
3091	s2 = Segment(b,d, tag = up)
3092	s3 = Segment(b,a, tag = left)
3093	s4 = Segment(d,c, tag = right)
3094	s5 = Segment(a,c, tag = down)
3095
3096	if regions:
3097	e = Vertex (side_length/2,side_length/2)
3098	s6 = Segment(a,e, tag = down + left)
3099	s7 = Segment(b,e, tag = up + left)
3100	s8 = Segment(c,e, tag = down + right)
3101	s9 = Segment(d,e, tag = up + right)
3102	r1 = Region(side_length/2,3.*side_length/4, tag = up)
3103	r2 = Region(1.*side_length/4,side_length/2, tag = left)
3104	r3 = Region(3.*side_length/4,side_length/2, tag = right)
3105	r4 = Region(side_length/2,1.*side_length/4, tag = down)
3106	mesh = Mesh(userVertices=[a,b,c,d,e],
3107	userSegments=[s2,s3,s4,s5,s6,s7,s8,s9],
3108	regions = [r1,r2,r3,r4])
3109	else:
3110	mesh = Mesh(userVertices=[a,b,c,d],
3111	userSegments=[s2,s3,s4,s5])
3112
3113	return mesh
3114
3115
3116
3117	def region_strings2ints(region_list):
3118	"""Given a list of (x_int,y_int,tag_string) lists it returns a list of
3119	(x_int,y_int,tag_int) and a list to convert the tag_int's back to
3120	the tag_strings
3121	"""
3122	# Make sure "" has an index of 0
3123	region_list.reverse()
3124	region_list.append((1.0,2.0,""))
3125	region_list.reverse()
3126	convertint2string = []
3127	for i in xrange(len(region_list)):
3128	convertint2string.append(region_list[i][2])
3129	if len(region_list[i]) == 4: # there's an area value
3130	region_list[i] = (region_list[i][0],
3131	region_list[i][1],i,region_list[i][3])
3132	elif len(region_list[i]) == 3: # no area value
3133	region_list[i] = (region_list[i][0],region_list[i][1],i)
3134	else:
3135	print "The region list has a bad size"
3136	# raise an error ..
3137	raise Error
3138
3139	#remove "" from the region_list
3140	region_list.pop(0)
3141
3142	return [region_list, convertint2string]
3143
3144	def region_ints2strings(region_list,convertint2string):
3145	"""Reverses the transformation of region_strings2ints
3146	"""
3147	if region_list == []:
3148	raise NoTrianglesError
3149	if region_list[0] != []:
3150	for i in xrange(len(region_list)):
3151	region_list[i] = [convertint2string[int(region_list[i][0])]]
3152	return region_list
3153
3154	def segment_ints2strings(intlist, convertint2string):
3155	"""Reverses the transformation of segment_strings2ints """
3156	stringlist = []
3157	for x in intlist:
3158	stringlist.append(convertint2string[x])
3159	return stringlist
3160
3161	def segment_strings2ints(stringlist, preset):
3162	"""Given a list of strings return a list of 0 to n ints which represent
3163	the strings and a converting list of the strings, indexed by 0 to n ints.
3164	Also, input an initial converting list of the strings
3165	Note, the converting list starts off with
3166	["internal boundary", "external boundary", "internal boundary"]
3167	example input and output
3168	input ["a","b","a","c"],["c"]
3169	output [[2, 1, 2, 0], ['c', 'b', 'a']]
3170
3171	the first element in the converting list will be
3172	overwritten with "".
3173	?This will become the third item in the converting list?
3174
3175	# note, order the initial converting list is important,
3176	since the index = the int tag
3177	"""
3178	nodups = unique(stringlist)
3179	# note, order is important, the index = the int tag
3180	#preset = ["internal boundary", "external boundary"]
3181	#Remove the preset tags from the list with no duplicates
3182	nodups = [x for x in nodups if x not in preset]
3183
3184	try:
3185	nodups.remove("") # this has to go to zero
3186	except ValueError:
3187	pass
3188
3189	# Add the preset list at the beginning of no duplicates
3190	preset.reverse()
3191	nodups.extend(preset)
3192	nodups.reverse()
3193
3194
3195	convertstring2int = {}
3196	convertint2string = []
3197	index = 0
3198	for x in nodups:
3199	convertstring2int[x] = index
3200	convertint2string.append(x)
3201	index += 1
3202	convertstring2int[""] = 0
3203
3204	intlist = []
3205	for x in stringlist:
3206	intlist.append(convertstring2int[x])
3207	return [intlist, convertint2string]
3208
3209
3210
3211
3212	def unique(s):
3213	"""Return a list of the elements in s, but without duplicates.
3214
3215	For example, unique([1,2,3,1,2,3]) is some permutation of [1,2,3],
3216	unique("abcabc") some permutation of ["a", "b", "c"], and
3217	unique(([1, 2], [2, 3], [1, 2])) some permutation of
3218	[[2, 3], [1, 2]].
3219
3220	For best speed, all sequence elements should be hashable. Then
3221	unique() will usually work in linear time.
3222
3223	If not possible, the sequence elements should enjoy a total
3224	ordering, and if list(s).sort() doesn't raise TypeError it's
3225	assumed that they do enjoy a total ordering. Then unique() will
3226	usually work in O(N*log2(N)) time.
3227
3228	If that's not possible either, the sequence elements must support
3229	equality-testing. Then unique() will usually work in quadratic
3230	time.
3231	"""
3232
3233	n = len(s)
3234	if n == 0:
3235	return []
3236
3237	# Try using a dict first, as that's the fastest and will usually
3238	# work. If it doesn't work, it will usually fail quickly, so it
3239	# usually doesn't cost much to try it. It requires that all the
3240	# sequence elements be hashable, and support equality comparison.
3241	u = {}
3242	try:
3243	for x in s:
3244	u[x] = 1
3245	except TypeError:
3246	del u # move on to the next method
3247	else:
3248	return u.keys()
3249
3250	# We can't hash all the elements. Second fastest is to sort,
3251	# which brings the equal elements together; then duplicates are
3252	# easy to weed out in a single pass.
3253	# NOTE: Python's list.sort() was designed to be efficient in the
3254	# presence of many duplicate elements. This isn't true of all
3255	# sort functions in all languages or libraries, so this approach
3256	# is more effective in Python than it may be elsewhere.
3257	try:
3258	t = list(s)
3259	t.sort()
3260	except TypeError:
3261	del t # move on to the next method
3262	else:
3263	assert n > 0
3264	last = t[0]
3265	lasti = i = 1
3266	while i < n:
3267	if t[i] != last:
3268	t[lasti] = last = t[i]
3269	lasti += 1
3270	i += 1
3271	return t[:lasti]
3272
3273	# Brute force is all that's left.
3274	u = []
3275	for x in s:
3276	if x not in u:
3277	u.append(x)
3278	return u
3279
3280	"""Refines triangles
3281
3282	Implements the #triangular bisection?# algorithm.
3283
3284
3285	"""
3286
3287	def Refine(mesh, triangles):
3288	"""
3289	Given a general_mesh, and a triangle number, split
3290	that triangle in the mesh in half. Then to prevent
3291	vertices and edges from meeting, keep refining
3292	neighbouring triangles until the mesh is clean.
3293	"""
3294	state = BisectionState(mesh)
3295	for triangle in triangles:
3296	if not state.refined_triangles.has_key(triangle):
3297	triangle.rotate_longest_side()
3298	state.start(triangle)
3299	Refine_mesh(mesh, state)
3300
3301	def Refine_mesh(mesh, state):
3302	"""
3303	"""
3304	state.getState(mesh)
3305	refine_triangle(mesh,state)
3306	state.evolve()
3307	if not state.end:
3308	Refine_mesh(mesh,state)
3309
3310	def refine_triangle(mesh,state):
3311	split(mesh,state.current_triangle,state.new_point)
3312	if state.case == 'one':
3313	state.r[3]=state.current_triangle#triangle 2
3314
3315	new_triangle_id = len(mesh.meshTriangles)-1
3316	new_triangle = mesh.meshTriangles[new_triangle_id]
3317
3318	split(mesh,new_triangle,state.old_point)
3319	state.r[2]=new_triangle#triangle 1.2
3320	state.r[4]=mesh.meshTriangles[len(mesh.meshTriangles)-1]#triangle 1.1
3321	r = state.r
3322	state.repairCaseOne()
3323
3324	if state.case == 'two':
3325	state.r[2]=mesh.meshTriangles[len(mesh.meshTriangles)-1]#triangle 1
3326
3327	new_triangle = state.current_triangle
3328
3329	split(mesh,new_triangle,state.old_point)
3330
3331	state.r[3]=mesh.meshTriangles[len(mesh.meshTriangles)-1]#triangle 2.1
3332	state.r[4]=new_triangle#triangle 2.2
3333	r = state.r
3334
3335	state.repairCaseTwo()
3336
3337	if state.case == 'vertex':
3338	state.r[2]=state.current_triangle#triangle 2
3339	state.r[3]=mesh.meshTriangles[len(mesh.meshTriangles)-1]#triangle 1
3340	r = state.r
3341	state.repairCaseVertex()
3342
3343	if state.case == 'start':
3344	state.r[2]=mesh.meshTriangles[len(mesh.meshTriangles)-1]#triangle 1
3345	state.r[3]=state.current_triangle#triangle 2
3346
3347	if state.next_case == 'boundary':
3348	state.repairCaseBoundary()
3349
3350
3351	def split(mesh, triangle, new_point):
3352	"""
3353	Given a mesh, triangle_id and a new point,
3354	split the corrosponding triangle into two
3355	new triangles and update the mesh.
3356	"""
3357
3358	new_triangle1 = Triangle(new_point,triangle.vertices[0],
3359	triangle.vertices[1],
3360	attribute = triangle.attribute,
3361	neighbors = None)
3362	new_triangle2 = Triangle(new_point,triangle.vertices[2],
3363	triangle.vertices[0],
3364	attribute = triangle.attribute,
3365	neighbors = None)
3366
3367	new_triangle1.setNeighbors(triangle.neighbors[2],None,new_triangle2)
3368	new_triangle2.setNeighbors(triangle.neighbors[1],new_triangle1,None)
3369
3370	mesh.meshTriangles.append(new_triangle1)
3371
3372	triangle.vertices = new_triangle2.vertices
3373	triangle.neighbors = new_triangle2.neighbors
3374
3375
3376	class State:
3377
3378	def __init__(self):
3379	pass
3380
3381	class BisectionState(State):
3382
3383
3384	def __init__(self,mesh):
3385	self.len = len(mesh.meshTriangles)
3386	self.refined_triangles = {}
3387	self.mesh = mesh
3388	self.current_triangle = None
3389	self.case = 'start'
3390	self.end = False
3391	self.r = [None,None,None,None,None]
3392
3393	def start(self, triangle):
3394	self.current_triangle = triangle
3395	self.case = 'start'
3396	self.end = False
3397	self.r = [None,None,None,None,None]
3398
3399	def getState(self,mesh):
3400	if not self.case == 'vertex':
3401	self.new_point=self.getNewVertex(mesh, self.current_triangle)
3402	#self.neighbour=self.getNeighbour(mesh, self.current_triangle)
3403	self.neighbour = self.current_triangle.neighbors[0]
3404	if not self.neighbour is None:
3405	self.neighbour.rotate_longest_side()
3406	self.next_case = self.get_next_case(mesh,self.neighbour,
3407	self.current_triangle)
3408	if self.case == 'vertex':
3409	self.new_point=self.old_point
3410
3411
3412	def evolve(self):
3413	if self.case == 'vertex':
3414	self.end = True
3415
3416	self.last_case = self.case
3417	self.case = self.next_case
3418
3419	self.old_point = self.new_point
3420	self.current_triangle = self.neighbour
3421
3422	if self.case == 'boundary':
3423	self.end = True
3424	self.refined_triangles[self.r[2]]=1
3425	self.refined_triangles[self.r[3]]=1
3426	if not self.r[4] is None:
3427	self.refined_triangles[self.r[4]]=1
3428	self.r[0]=self.r[2]
3429	self.r[1]=self.r[3]
3430
3431
3432	def getNewVertex(self,mesh,triangle):
3433	coordinate1 = triangle.vertices[1]
3434	coordinate2 = triangle.vertices[2]
3435	a = ([coordinate1.x1.,coordinate1.y1.])
3436	b = ([coordinate2.x1.,coordinate2.y1.])
3437	attributes = []
3438	for i in range(len(coordinate1.attributes)):
3439	att = (coordinate1.attributes[i]+coordinate2.attributes[i])/2
3440	attributes.append(att)
3441	new_coordinate = [((a[0]-b[0])/2+b[0]),((a[1]-b[1])/2+b[1])]
3442	newVertex = Vertex(new_coordinate[0],new_coordinate[1],
3443	attributes = attributes)
3444	mesh.maxVertexIndex+=1
3445	newVertex.index = mesh.maxVertexIndex
3446	mesh.meshVertices.append(newVertex)
3447	return newVertex
3448
3449	def get_next_case(self, mesh,neighbour,triangle):
3450	"""
3451	Given the locations of two neighbouring triangles,
3452	examine the interior indices of their vertices (i.e.
3453	0,1 or 2) to determine what how the neighbour needs
3454	to be refined.
3455	"""
3456	if (neighbour is None):
3457	next_case = 'boundary'
3458	else:
3459	if triangle.vertices[1].x==neighbour.vertices[2].x:
3460	if triangle.vertices[1].y==neighbour.vertices[2].y:
3461	next_case = 'vertex'
3462	if triangle.vertices[1].x==neighbour.vertices[0].x:
3463	if triangle.vertices[1].y==neighbour.vertices[0].y:
3464	next_case = 'two'
3465	if triangle.vertices[1].x==neighbour.vertices[1].x:
3466	if triangle.vertices[1].y==neighbour.vertices[1].y:
3467	next_case = 'one'
3468	return next_case
3469
3470
3471
3472	def repairCaseVertex(self):
3473
3474	r = self.r
3475
3476
3477	self.link(r[0],r[2])
3478	self.repair(r[0])
3479
3480	self.link(r[1],r[3])
3481	self.repair(r[1])
3482
3483	self.repair(r[2])
3484
3485	self.repair(r[3])
3486
3487
3488	def repairCaseOne(self):
3489	r = self.rkey
3490
3491
3492	self.link(r[0],r[2])
3493	self.repair(r[0])
3494
3495	self.link(r[1],r[4])
3496	self.repair(r[1])
3497
3498	self.repair(r[4])
3499
3500	def repairCaseTwo(self):
3501	r = self.r
3502
3503	self.link(r[0],r[4])
3504	self.repair(r[0])
3505
3506	self.link(r[1],r[3])
3507	self.repair(r[1])
3508
3509	self.repair(r[4])
3510
3511	def repairCaseBoundary(self):
3512	r = self.r
3513	self.repair(r[2])
3514	self.repair(r[3])
3515
3516
3517
3518	def repair(self,triangle):
3519	"""
3520	Given a triangle that knows its neighbours, this will
3521	force the neighbours to comply.
3522
3523	However, it needs to compare the vertices of triangles
3524	for this implementation
3525
3526	But it doesn't work for invalid neighbour structures
3527	"""
3528	n=triangle.neighbors
3529	for i in (0,1,2):
3530	if not n[i] is None:
3531	for j in (0,1,2):#to find which side of the list is broken
3532	if not (n[i].vertices[j] in triangle.vertices):
3533	#ie if j is the side of n that needs fixing
3534	k = j
3535	n[i].neighbors[k]=triangle
3536
3537
3538
3539	def link(self,triangle1,triangle2):
3540	"""
3541	make triangle1 neighbors point to t
3542	#count = 0riangle2
3543	"""
3544	count = 0
3545	for i in (0,1,2):#to find which side of the list is broken
3546	if not (triangle1.vertices[i] in triangle2.vertices):
3547	j = i
3548	count+=1
3549	assert count == 1
3550	triangle1.neighbors[j]=triangle2
3551
3552	class Discretised_Tuple_Set:
3553	"""
3554	if a={(0.0):[(0.01),(0.02)],(0.2):[(0.17)]}
3555	a[(0.01)]=a[(0.0)]=[(0.01),(0.02)]
3556	a[(10000)]=[] #NOT KEYERROR
3557
3558	a.append[(0.01)]
3559	=> {0.0:[(0.01),(0.02),(0.01)],0.2:[(0.17)]}
3560
3561	#NOT IMPLIMENTED
3562	a.remove[(0.01)]
3563	=> {(0.0):[(0.02),(0.01)],0.2:[(0.17)]}
3564
3565	a.remove[(0.17)]
3566	=> {(0.0):[(0.02),(0.01)],0.2:[]}
3567	#NOT IMPLIMENTED
3568	at a.precision = 2:
3569	a.round_up_rel[0.0]=
3570	a.round_flat[0.0]=
3571	a.round_down_rel[0.0]=
3572
3573	a.up((0.1,2.04))=
3574
3575	If t_rel = 0, nothing gets rounded into
3576	two bins. If t_rel = 0.5, everything does.
3577
3578	Ideally, precision can be set high, so that multiple
3579	entries are rarely in the same bin. And t_rel should
3580	be low (<0.1 for 1 dimension!,<(0.1/n) for small n!!)
3581	so that it is rare to put items in mutiple bins.
3582
3583	Ex bins per entry = product(a,b,c...,n)
3584	a = 1 or 2 s.t. Ex(a) = 1+2*t_rel
3585	b = 1 or 2 ...
3586
3587	BUT!!! to avoid missing the right bin:
3588	(-10)*(precision+1)t_rel must be greater than the
3589	greatest possible variation that an identical element
3590	can display.
3591
3592
3593	Note that if tol = 0.5 (the max allowed) 0.6 will round to .7 and .5
3594	but not .6 - this looks wrong, but note that everything will round,
3595	so .6 wont be missed as everything close to it will check in .7 and .5.
3596	"""
3597	def __init__(self,p_rel = 6,t_rel = 0.01):
3598	self.__p_rel__ = p_rel
3599	self.__t_rel__ = t_rel
3600
3601	self.__p_abs__ = p_rel+1
3602	self.__t_abs__ = t_rel
3603
3604	assert t_rel <= 0.5
3605	self.__items__ = {}
3606	from math import frexp
3607	self.frexp = frexp
3608	roundings = [self.round_up_rel,\
3609	self.round_down_rel,self.round_flat_rel,\
3610	self.round_down_abs,self.round_up_abs,\
3611	self.round_flat_abs]#
3612
3613	self.roundings = roundings
3614
3615	def __repr__(self):
3616	return '%s'%self.__items__
3617
3618	def rounded_keys(self,key):
3619	key = tuple(key)
3620	keys = [key]
3621	keys = self.__rounded_keys__(key)
3622	return (keys)
3623
3624	def __rounded_keys__(self,key):
3625	keys = []
3626	rounded_key=list(key)
3627	rounded_values=list(key)
3628
3629	roundings = list(self.roundings)
3630
3631	#initialise rounded_values
3632	round = roundings.pop(0)
3633	for i in range(len(rounded_values)):
3634	rounded_key[i]=round(key[i])
3635	rounded_values[i]={}
3636	rounded_values[i][rounded_key[i]]=None
3637	keys.append(tuple(rounded_key))
3638
3639	for round in roundings:
3640	for i in range(len(rounded_key)):
3641	rounded_value=round(key[i])
3642	if not rounded_values[i].has_key(rounded_value):
3643	#ie unless round_up_rel = round_down_rel
3644	#so the keys stay unique
3645	for j in range(len(keys)):
3646	rounded_key = list(keys[j])
3647	rounded_key[i]=rounded_value
3648	keys.append(tuple(rounded_key))
3649	return keys
3650
3651	def append(self,item):
3652	keys = self.rounded_keys(item)
3653	for key in keys:
3654	if self.__items__.has_key(key):
3655	self.__items__[key].append(item)
3656	else:
3657	self.__items__[key]=[item]
3658
3659	def __getitem__(self,key):
3660	answer = []
3661	keys = self.rounded_keys(key)
3662	for key in keys:
3663	if self.__items__.has_key(key):
3664	answer.extend(self.__items__[key])
3665	#if len(answer)==0:
3666	# raise KeyError#FIXME or return KeyError
3667	# #FIXME or just return []?
3668	else:
3669	return answer #FIXME or unique(answer)?
3670
3671	def __delete__(self,item):
3672	keys = self.rounded_keys(item)
3673	answer = False
3674	#if any of the possible keys contains
3675	#a list, return true
3676	for key in keys:
3677	if self.__items__.has_key(key):
3678	if item in self.__items__[key]:
3679	self.__items__[key].remove(item)
3680
3681	def remove(self,item):
3682	self.__delete__(item)
3683
3684	def __contains__(self,item):
3685
3686	keys = self.rounded_keys(item)
3687	answer = False
3688	#if any of the possible keys contains
3689	#a list, return true
3690	for key in keys:
3691	if self.__items__.has_key(key):
3692	if item in self.__items__[key]:
3693	answer = True
3694	return answer
3695
3696
3697	def has_item(self,item):
3698	return self.__contains__(item)
3699
3700	def round_up_rel2(self,value):
3701	t_rel=self.__t_rel__
3702	#Rounding up the value
3703	m,e = self.frexp(value)
3704	m = m/2
3705	e = e + 1
3706	#m is the mantissa, e the exponent
3707	# 0.5 < \|m\| < 1.0
3708	m = m+t_rel(10*-(self.__p_rel__))
3709	#bump m up
3710	m = round(m,self.__p_rel__)
3711	return m(2.*e)
3712
3713	def round_down_rel2(self,value):
3714	t_rel=self.__t_rel__
3715	#Rounding down the value
3716	m,e = self.frexp(value)
3717	m = m/2
3718	e = e + 1
3719	#m is the mantissa, e the exponent
3720	# 0.5 < m < 1.0
3721	m = m-t_rel(10*-(self.__p_rel__))
3722	#bump the \|m\| down, by 5% or whatever
3723	#self.p_rel dictates
3724	m = round(m,self.__p_rel__)
3725	return m(2.*e)
3726
3727	def round_flat_rel2(self,value):
3728	#redundant
3729	m,e = self.frexp(value)
3730	m = m/2
3731	e = e + 1
3732	m = round(m,self.__p_rel__)
3733	return m(2.*e)
3734
3735	def round_up_rel(self,value):
3736	t_rel=self.__t_rel__
3737	#Rounding up the value
3738	m,e = self.frexp(value)
3739	#m is the mantissa, e the exponent
3740	# 0.5 < \|m\| < 1.0
3741	m = m+t_rel(10*-(self.__p_rel__))
3742	#bump m up
3743	m = round(m,self.__p_rel__)
3744	return m(2.*e)
3745
3746	def round_down_rel(self,value):
3747	t_rel=self.__t_rel__
3748	#Rounding down the value
3749	m,e = self.frexp(value)
3750	#m is the mantissa, e the exponent
3751	# 0.5 < m < 1.0
3752	m = m-t_rel(10*-(self.__p_rel__))
3753	#bump the \|m\| down, by 5% or whatever
3754	#self.p_rel dictates
3755	m = round(m,self.__p_rel__)
3756	return m(2.*e)
3757
3758	def round_flat_rel(self,value):
3759	#redundant
3760	m,e = self.frexp(value)
3761	m = round(m,self.__p_rel__)
3762	return m(2.*e)
3763
3764	def round_up_abs(self,value):
3765	t_abs=self.__t_abs__
3766	#Rounding up the value
3767	m = value+t_abs(10*-(self.__p_abs__))
3768	#bump m up
3769	m = round(m,self.__p_abs__)
3770	return m
3771
3772	def round_down_abs(self,value):
3773	t_abs=self.__t_abs__
3774	#Rounding down the value
3775	m = value-t_abs(10*-(self.__p_abs__))
3776	#bump the \|m\| down, by 5% or whatever
3777	#self.p_rel dictates
3778	m = round(m,self.__p_abs__)
3779	return m
3780
3781	def round_flat_abs(self,value):
3782	#redundant?
3783	m = round(value,self.__p_abs__)
3784	return m
3785
3786	def keys(self):
3787	return self.__items__.keys()
3788
3789
3790	class Mapped_Discretised_Tuple_Set(Discretised_Tuple_Set):
3791	"""
3792	This is a discretised tuple set, but
3793	based on a mapping. The mapping MUST
3794	return a sequence.
3795
3796	example:
3797	def weight(animal):
3798	return [animal.weight]
3799
3800	a = Mapped_Discretised_Tuple_Set(weight)
3801	a.append[cow]
3802	a.append[fox]
3803	a.append[horse]
3804
3805	a[horse] -> [cow,horse]
3806	a[dog] -> [fox]
3807	a[elephant] -> []
3808	"""
3809	def __init__(self,mapping,p_rel = 6, t_rel=0.01):
3810	Discretised_Tuple_Set.__init__\
3811	(self,p_rel,t_rel = t_rel)
3812	self.mapping = mapping
3813
3814	def rounded_keys(self,key):
3815	mapped_key = tuple(self.mapping(key))
3816	keys = self.__rounded_keys__(mapped_key)
3817	return keys
3818
3819	class Affine_Linespace(Mapped_Discretised_Tuple_Set):
3820	"""
3821	The affine linespace creates a way to record and compare lines.
3822	Precision is a bit of a hack, but it creates a way to avoid
3823	misses caused by round offs (between two lines of different
3824	lenghts, the short one gets rounded off more).
3825	I am starting to think that a quadratic search would be faster.
3826	Nearly.
3827	"""
3828	def __init__(self,p_rel=4,t_rel=0.2):
3829	Mapped_Discretised_Tuple_Set.__init__\
3830	(self,self.affine_line,\
3831	p_rel=p_rel,t_rel=t_rel)
3832
3833	roundings = \
3834	[self.round_down_rel,self.round_up_rel,self.round_flat_rel]
3835	self.roundings = roundings
3836	#roundings = \
3837	#[self.round_down_abs,self.round_up_abs,self.round_flat_abs]
3838	#self.roundings = roundings
3839
3840	def affine_line(self,line):
3841	point_1 = line[0]
3842	point_2 = line[1]
3843	#returns the equation of a line
3844	#between two points, in the from
3845	#(a,b,-c), as in ax+by-c=0
3846	#or line dot (x,y,1) = (0,0,0)
3847
3848	#Note that it normalises the line
3849	#(a,b,-c) so Line*Line = 1.
3850	#This does not change the mathematical
3851	#properties, but it makes comparism
3852	#easier.
3853
3854	#There are probably better algorithms.
3855	x1 = point_1[0]
3856	x2 = point_2[0]
3857	y1 = point_1[1]
3858	y2 = point_2[1]
3859	dif_x = x1-x2
3860	dif_y = y1-y2
3861
3862	if dif_x == dif_y == 0:
3863	msg = 'points are the same'
3864	raise msg
3865	elif abs(dif_x)>=abs(dif_y):
3866	alpha = (-dif_y)/dif_x
3867	#a = alpha * b
3868	b = -1.
3869	c = (x1alpha+x2alpha+y1+y2)/2.
3870	a = alpha*b
3871	else:
3872	beta = dif_x/(-dif_y)
3873	#b = beta * a
3874	a = 1.
3875	c = (x1+x2+y1beta+y2beta)/2.
3876	b = beta*a
3877	mag = abs(a)+abs(b)
3878	#This does not change the mathematical
3879	#properties, but it makes comparism possible.
3880
3881	#note that the gradient is b/a, or (a/b)**-1.
3882	#so
3883
3884	#if a == 0:
3885	# sign_a = 1.
3886	#else:
3887	# sign_a = a/((a2)0.5)
3888	#if b == 0:
3889	# sign_b = 1.
3890	#else:
3891	# sign_b = b/((b2)0.5)
3892	#if c == 0:
3893	# sign_c = 1.
3894	#else:
3895	# sign_c = c/((c2)0.5)
3896	#a = a/mag*sign_a
3897	#b = b/mag*sign_b
3898	#c = c/mag*sign_c
3899	a = a/mag
3900	b = b/mag
3901	c = c/mag
3902	return a,b,c
3903
3904
3905	# FIXME (DSG-DSG)
3906	# To do-
3907	# Create a clear interface. eg
3908	# have the interface methods more at the top of this file and add comments
3909	# for the interface functions/methods, use function_name (not functionName),
3910
3911	#Currently
3912	#function_name methods assume absolute values. Geo-refs can be passed in.
3913	#
3914
3915	# instead of functionName
3916	if __name__ == "__main__":
3917	#from mesh import *
3918	# THIS CAN BE DELETED
3919	m = Mesh()
3920	dict = importUngenerateFile("ungen_test.txt")
3921	m.addVertsSegs(dict)
3922	print m3

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

Download in other formats: