Context Navigation

source: inundation/pmesh/mesh.py @ 2243

Last change on this file since 2243 was 2243, checked in by duncan, 18 years ago
comments
File size: 118.7 KB

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

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

Download in other formats: