Context Navigation

source: inundation/pmesh/mesh.py @ 2152

Last change on this file since 2152 was 2144, checked in by duncan, 19 years ago
different ordering found in windows and linux with dict.value(). change test so order is not tested
File size: 119.8 KB

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

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

Download in other formats: