Context Navigation

source: inundation/pmesh/mesh.py @ 2376

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

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

Download in other formats: