Context Navigation

source: inundation/pmesh/mesh.py @ 2891

Last change on this file since 2891 was 2778, checked in by steve, 19 years ago
On our 64 bit machine, ran into problem in pmesh/mesh.py where seg[0], seg[1], triangle, and neighbor where not seen as integers. Had to wrap them in int(..)
File size: 124.0 KB

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

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

Download in other formats: