Context Navigation

source: inundation/pmesh/mesh.py @ 2152

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

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

Download in other formats: