Context Navigation

source: trunk/anuga_core/source/anuga/pmesh/mesh.py @ 7873

Last change on this file since 7873 was 7873, checked in by hudson, 14 years ago
Added ungenerate loading functionality.
File size: 84.8 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
[4144]	7	should not be used as persistant information. It is not the 'index' of
[2525]	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
	16	import sys
	17	import math
	18	import re
	19	import os
	20	import pickle
	21
	22	import types
	23	import exceptions
	24
[7276]	25	import numpy as num
[4144]	26
[7276]	27
[4665]	28	from anuga.coordinate_transforms.geo_reference import Geo_reference, \
	29	DEFAULT_ZONE
[4182]	30	from anuga.load_mesh.loadASCII import NOMAXAREA, export_mesh_file, \
	31	import_mesh_file
[3514]	32	import anuga.alpha_shape.alpha_shape
[4165]	33	from anuga.geospatial_data.geospatial_data import Geospatial_data, \
	34	ensure_geospatial, ensure_absolute, ensure_numeric
[4144]	35	from anuga.mesh_engine.mesh_engine import generate_mesh
[7317]	36	import anuga.utilities.log as log
[2525]	37
[7873]	38	from anuga.file.ungenerate import load_ungenerate
[7317]	39
[3193]	40	try:
	41	import kinds
	42	except ImportError:
[4665]	43	# Hand-built mockup of the things we need from the kinds package, since it
[7276]	44	# was recently removed from the standard numeric distro. Some users may
[3193]	45	# not have it by default.
	46	class _bunch:
	47	pass
	48
	49	class _kinds(_bunch):
	50	default_float_kind = _bunch()
	51	default_float_kind.MIN = 2.2250738585072014e-308 #smallest +ve number
	52	default_float_kind.MAX = 1.7976931348623157e+308
	53
	54	kinds = _kinds()
	55
[2525]	56
	57	# 1st and third values must be the same
	58	# FIXME: maybe make this a switch that the user can change? - DSG
	59	initialconversions = ['', 'exterior', '']
[5188]	60	SEG_COLOUR = 'blue'
[2525]	61
	62
	63	class MeshObject:
	64	"""
	65	An abstract superclass for the basic mesh objects, eg vertex, segment,
	66	triangle.
	67	"""
	68	def __init__(self):
	69	pass
	70
	71	class Point(MeshObject):
	72	"""
	73	Define a point in a 2D space.
	74	"""
	75	def __init__(self,X,Y):
	76	__slots__ = ['x','y']
	77	self.x=X
	78	self.y=Y
	79
	80	def DistanceToPoint(self, OtherPoint):
	81	"""
	82	Returns the distance from this point to another
	83	"""
[4665]	84	SumOfSquares = ((self.x - OtherPoint.x)**2) + \
	85	((self.y - OtherPoint.y)**2)
[2525]	86	return math.sqrt(SumOfSquares)
	87
	88	def IsInsideCircle(self, Center, Radius):
	89	"""
	90	Return 1 if this point is inside the circle,
	91	0 otherwise
	92	"""
	93
	94	if (self.DistanceToPoint(Center)<Radius):
	95	return 1
	96	else:
	97	return 0
	98
	99	def __repr__(self):
	100	return "(%f,%f)" % (self.x,self.y)
	101
	102	def cmp_xy(self, point):
	103	if self.x < point.x:
	104	return -1
	105	elif self.x > point.x:
	106	return 1
	107	else:
	108	if self.y < point.y:
	109	return -1
	110	elif self.y > point.y:
	111	return 1
	112	else:
	113	return 0
	114
	115	def same_x_y(self, point):
	116	if self.x == point.x and self.y == point.y:
	117	return True
	118	else:
	119	return False
	120
	121
	122
	123	class Vertex(Point):
[7276]	124	"""A point on the mesh.
[2525]	125	Object attributes based on the Triangle program
	126	"""
[7276]	127
	128	VERTEXSQUARESIDELENGTH = 6
	129
	130	def __init__(self, X, Y, attributes=None):
[2525]	131	__slots__ = ['x','y','attributes']
[7276]	132
	133	# we don't care what we get, as long as we can get float value
	134	self.x = float(X)
	135	self.y = float(Y)
	136
	137	self.attributes = []
	138	if not attributes is None:
	139	self.attributes = attributes
[2525]	140
[7276]	141	def setAttributes(self,attributes):
	142	"""attributes is a list. """
[2525]	143
	144	self.attributes = attributes
	145
[4665]	146	def draw(self, canvas, tags, colour = 'black',scale = 1, xoffset = 0,
	147	yoffset =0, ):
[2525]	148	x = scale*(self.x + xoffset)
	149	y = -1scale(self.y + yoffset) # - since for a canvas - is up
	150	cornerOffset= self.VERTEXSQUARESIDELENGTH/2
	151
	152	# A hack to see the vert tags
	153	# note: there will be many tags, since tags will not be removed
	154	#when zooming
	155	#canvas.create_text(x+ 2*cornerOffset,
	156	# y+ 2*cornerOffset,
	157	# text=tags)
	158
[5188]	159	# This gives points info. It is a mess though, since numbers are
	160	# not deleted when zaooming in.
	161	#canvas.create_text(x+ 2*cornerOffset,
	162	# y+ 2*cornerOffset,
	163	# text=str(x)+','+str(y))
	164
[2525]	165	return canvas.create_rectangle(x-cornerOffset,
	166	y-cornerOffset,
	167	x+cornerOffset,
	168	y+cornerOffset,
	169	tags = tags,
	170	outline=colour,
	171	fill = 'white')
	172
	173	#return tags
	174
	175	def __repr__(self):
	176	return "[(%f,%f),%r]" % (self.x,self.y,self.attributes)
[7276]	177
	178
[2525]	179	class Hole(Point):
[7276]	180	"""A region of the mesh were no triangles are generated.
[2525]	181	Defined by a point in the hole enclosed by segments.
	182	"""
	183
	184	HOLECORNERLENGTH = 6
	185
[4665]	186	def draw(self, canvas, tags, colour = 'purple',scale = 1, xoffset = 0,
	187	yoffset =0, ):
[2525]	188	x = scale*(self.x + xoffset)
	189	y = -1scale(self.y + yoffset) # - since for a canvas - is up
	190	cornerOffset= self.HOLECORNERLENGTH/2
	191	return canvas.create_oval(x-cornerOffset,
	192	y-cornerOffset,
	193	x+cornerOffset,
	194	y+cornerOffset,
	195	tags = tags,
	196	outline=colour,
	197	fill = 'white')
	198
	199	class Region(Point):
[7276]	200	"""A region of the mesh.
	201	Defined by a point in the region enclosed by segments. Used to tag areas.
[2525]	202	"""
[7276]	203
[2525]	204	CROSSLENGTH = 6
	205	TAG = 0
	206	MAXAREA = 1
	207
[7276]	208	def __init__(self, X, Y, tag=None, maxArea=None):
	209	"""Precondition: tag is a string and maxArea is a real"""
	210
[2525]	211	self.x=X
	212	self.y=Y
[7276]	213	self.attributes = [] # index 0 is the tag string
	214	# optional index 1 is the max triangle area
	215	# NOTE the size of this attribute is assumed
	216	# to be 1 or 2 in regionstrings2int
[2525]	217	if tag is None:
	218	self.attributes.append("")
	219	else:
	220	self.attributes.append(tag) #this is a string
	221
	222	if maxArea is not None:
	223	self.setMaxArea(maxArea) # maxArea is a number
	224
	225	def getTag(self,):
	226	return self.attributes[self.TAG]
	227
	228	def setTag(self,tag):
	229	self.attributes[self.TAG] = tag
	230
	231	def getMaxArea(self):
	232	""" Returns the Max Area of a Triangle or
	233	None, if the Max Area has not been set.
	234	"""
	235	if self.isMaxArea():
	236	return self.attributes[1]
	237	else:
	238	return None
	239
	240	def setMaxArea(self,MaxArea):
[3193]	241	if MaxArea is not None:
	242	if self.isMaxArea():
	243	self.attributes[self.MAXAREA] = float(MaxArea)
	244	else:
	245	self.attributes.append( float(MaxArea) )
[2525]	246
	247	def deleteMaxArea(self):
	248	if self.isMaxArea():
	249	self.attributes.pop(self.MAXAREA)
	250
	251	def isMaxArea(self):
	252	return len(self.attributes)> 1
	253
[4665]	254	def draw(self, canvas, tags, scale=1, xoffset = 0, yoffset =0,
	255	colour = "black"):
[2525]	256	"""
	257	Draw a black cross, returning the objectID
	258	"""
	259	x = scale*(self.x + xoffset)
	260	y = -1scale(self.y + yoffset)
	261	cornerOffset= self.CROSSLENGTH/2
	262	return canvas.create_polygon(x,
	263	y-cornerOffset,
	264	x,
	265	y,
	266	x+cornerOffset,
	267	y,
	268	x,
	269	y,
	270	x,
	271	y+cornerOffset,
	272	x,
	273	y,
	274	x-cornerOffset,
	275	y,
	276	x,
	277	y,
	278	tags = tags,
	279	outline = colour,fill = '')
	280
	281	def __repr__(self):
	282	if self.isMaxArea():
	283	area = self.getMaxArea()
	284	return "(%f,%f,%s,%f)" % (self.x,self.y,
	285	self.getTag(), self.getMaxArea())
	286	else:
	287	return "(%f,%f,%s)" % (self.x,self.y,
	288	self.getTag())
	289
	290	class Segment(MeshObject):
	291	"""
	292	Segments are edges whose presence in the triangulation is enforced.
	293
	294	"""
	295	def __init__(self, vertex1, vertex2, tag = None ):
	296	"""
	297	Each segment is specified by listing the vertices of its endpoints
	298	The vertices are Vertex objects
	299	"""
	300	assert(vertex1 != vertex2)
	301	self.vertices = [vertex1,vertex2 ]
	302
	303	if tag is None:
	304	self.tag = self.__class__.default
	305	else:
	306	self.tag = tag #this is a string
	307
	308	def __repr__(self):
	309	return "[%s,%s]" % (self.vertices,self.tag)
	310
	311
[5188]	312	def draw(self, canvas, tags,scale=1, xoffset=0,
	313	yoffset=0, colour=SEG_COLOUR):
[2525]	314	x=[]
	315	y=[]
	316	for end in self.vertices:
	317	#end.draw(canvas,scale, xoffset, yoffset ) # draw the vertices
	318	x.append(scale*(end.x + xoffset))
	319	y.append(-1scale(end.y + yoffset)) # - since for a canvas - is up
	320
	321	return canvas.create_line(x[0], y[0], x[1], y[1],
	322	tags = tags,fill=colour)
	323	def set_tag(self,tag):
	324	self.tag = tag
	325
	326	# Class methods
	327	def set_default_tag(cls, default):
	328	cls.default = default
	329
	330	def get_default_tag(cls):
	331	return cls.default
	332
	333	set_default_tag = classmethod(set_default_tag)
	334	get_default_tag = classmethod(get_default_tag)
	335
	336	Segment.set_default_tag("")
	337
[4144]	338	class Rigid_triangulation:
	339	"""
	340	This is a triangulation that can't have triangles added or taken away.
	341
	342	It just represents the triangulation, not the mesh outline needed to
	343	build the triangulation.
	344
	345	This is the guts of the data structure;
	346	generated vertex list: [(x1,y1),(x2,y2),...] (Tuples of doubles)
	347	generated segment list: [(point1,point2),(p3,p4),...]
	348	(Tuples of integers)
	349	generated segment tag list: [tag,tag,...] list of strings
	350
	351	generated triangle list: [(p1,p2,p3), (p4,p5,p6),....] tuple of points
	352
	353	generated triangle attribute list: [s1,s2,...] list of strings
	354
	355	generated triangle neighbor list: [(t1,t2,t3), (t4,t5,t6),....]
	356	tuple of triangles
	357
	358	Should I do the basic structure like triangle, general
	359	mesh or something else? How about like triangle, since
	360	that's where the info is from, and it'll fit easier into
	361	this file..
	362
	363	Removing the loners is difficult, since all the vert's
	364	after it must be removed.
	365
	366	This happens in set_triangulation.
	367
	368	"""
	369
	370	def __init__(self,
	371	triangles,
	372	segments,
	373	vertices,
	374	triangle_tags,
[4901]	375	triangle_neighbors,
[4144]	376	segment_tags,
[4902]	377	vertex_attributes,
	378	vertex_attribute_titles=None
[4144]	379	):
	380
[4901]	381	self.triangles = ensure_numeric(triangles)
	382	self.triangle_neighbors = ensure_numeric(triangle_neighbors)
	383	self.triangle_tags = triangle_tags # list of strings
	384	self.segments = ensure_numeric(segments)
[4144]	385
[4901]	386	self.segment_tags = segment_tags # list of strings
	387	self.vertices = ensure_numeric(vertices)
[4902]	388	# This is needed for __cmp__
	389	if vertex_attributes is None:
	390	self.vertex_attributes = []
	391	else:
	392	self.vertex_attributes = ensure_numeric(vertex_attributes)
	393	if vertex_attribute_titles is None:
	394	self.vertex_attribute_titles = []
	395	else:
	396	self.vertex_attribute_titles = vertex_attribute_titles
[4901]	397
	398	def draw_triangulation(self, canvas, scale=1, xoffset=0, yoffset=0,
	399	colour="green"):
	400	"""
	401	Draw a triangle, returning the objectID
	402	"""
[4144]	403
[4901]	404	# FIXME(DSG-DSG) This could be a data structure that is
	405	# remembered and doesn't have any duplicates.
	406	# but I wouldn't be able to use create_polygon.
	407	# regard it as drawing a heap of segments
[4144]	408
[4901]	409	for tri in self.triangles:
	410	vertices = []
	411	for v_index in range(3):
	412	vertices.append(self.vertices[tri[v_index]])
	413
	414	objectID = canvas.create_polygon(
	415	scale*(vertices[1][0] + xoffset),
	416	scale-1(vertices[1][1] + yoffset),
	417	scale*(vertices[0][0] + xoffset),
	418	scale-1(vertices[0][1] + yoffset),
	419	scale*(vertices[2][0] + xoffset),
	420	scale-1(vertices[2][1] + yoffset),
	421	outline = colour,fill = ''
	422	)
[4905]	423
	424	def calc_mesh_area(self):
	425	area = 0
	426	for tri in self.triangles:
	427	vertices = []
	428	for v_index in range(3):
	429	vertices.append(self.vertices[tri[v_index]])
	430	ax = vertices[0][0]
	431	ay = vertices[0][1]
	432
	433	bx = vertices[1][0]
	434	by = vertices[1][1]
	435
	436	cx = vertices[2][0]
	437	cy = vertices[2][1]
	438
	439	area += abs((bxay-axby)+(cxby-bxcy)+(axcy-cxay))/2
	440	return area
[4901]	441
[2525]	442	class Mesh:
	443	"""
	444	Representation of a 2D triangular mesh.
	445	User attributes describe the mesh region/segments/vertices/attributes
	446
[4665]	447	mesh attributes describe the mesh that is produced eg triangles and
	448	vertices.
[2525]	449	All point information is relative to the geo_reference passed in
	450
	451
	452	"""
	453
	454	def __repr__(self):
	455	return """
	456	mesh Triangles: %s
	457	mesh Attribute Titles: %s
	458	mesh Segments: %s
	459	mesh Vertices: %s
	460	user Segments: %s
	461	user Vertices: %s
	462	holes: %s
	463	regions: %s""" % (self.meshTriangles,
	464	self.attributeTitles,
	465	self.meshSegments,
	466	self.meshVertices,
	467	self.getUserSegments(),
	468	self.userVertices,
	469	self.holes,
	470	self.regions)
	471
	472	def __init__(self,
	473	userSegments=None,
	474	userVertices=None,
	475	holes=None,
	476	regions=None,
	477	geo_reference=None):
	478	self.meshTriangles=[]
	479	self.attributeTitles=[]
	480	self.meshSegments=[]
	481	self.meshVertices=[]
	482
[4901]	483	# Rigid
	484	self.tri_mesh=None
	485
[2525]	486
	487	self.visualise_graph = True
	488
	489	if userSegments is None:
	490	self.userSegments=[]
	491	else:
	492	self.userSegments=userSegments
	493	self.alphaUserSegments=[]
	494
	495	if userVertices is None:
	496	self.userVertices=[]
	497	else:
	498	self.userVertices=userVertices
	499
	500	if holes is None:
	501	self.holes=[]
	502	else:
	503	self.holes=holes
	504
	505	if regions is None:
	506	self.regions=[]
	507	else:
	508	self.regions=regions
	509
	510	if geo_reference is None:
	511	self.geo_reference = Geo_reference(DEFAULT_ZONE,0,0)
	512	else:
	513	self.geo_reference = geo_reference
[3435]	514
	515	self.shape = None
[2525]	516
	517	def __cmp__(self,other):
	518
	519	# A dic for the initial m
	520	dic = self.Mesh2triangList()
	521	dic_mesh = self.Mesh2MeshList()
	522	for element in dic_mesh.keys():
	523	dic[element] = dic_mesh[element]
	524	for element in dic.keys():
	525	dic[element].sort()
	526
	527	# A dic for the exported/imported m
	528	dic_other = other.Mesh2triangList()
	529	dic_mesh = other.Mesh2MeshList()
	530	for element in dic_mesh.keys():
	531	dic_other[element] = dic_mesh[element]
	532	for element in dic.keys():
	533	dic_other[element].sort()
	534
	535	#print "dsg************************8"
	536	#print "dic ",dic
	537	#print "*******8"
	538	#print "mesh",dic_other
	539	#print "dic.__cmp__(dic_o)",dic.__cmp__(dic_other)
	540	#print "dsg************************8"
	541
	542	return (dic.__cmp__(dic_other))
	543
	544	def addUserPoint(self, pointType, x,y):
	545	if pointType == Vertex:
	546	point = self.addUserVertex(x,y)
	547	if pointType == Hole:
	548	point = self._addHole(x,y)
	549	if pointType == Region:
	550	point = self._addRegion(x,y)
	551	return point
	552
	553	def addUserVertex(self, x,y):
	554	v=Vertex(x, y)
	555	self.userVertices.append(v)
	556	return v
	557
	558	def _addHole(self, x,y):
	559	h=Hole(x, y)
	560	self.holes.append(h)
	561	return h
	562
	563	def add_hole(self, x,y, geo_reference=None):
	564	"""
	565	adds a point, which represents a hole.
	566
	567	The point data can have it's own geo_refernece.
	568	If geo_reference is None the data is asumed to be absolute
	569	"""
	570	[[x,y]] = self.geo_reference.change_points_geo_ref([x,y],
	571	points_geo_ref=geo_reference)
	572	return self._addHole(x, y)
	573
	574	def _addRegion(self, x,y):
	575	h=Region(x, y)
	576	self.regions.append(h)
	577	return h
	578
[5207]	579	def add_region(self, x,y, geo_reference=None, tag=None):
[2525]	580	"""
	581	adds a point, which represents a region.
	582
	583	The point data can have it's own geo_refernece.
	584	If geo_reference is None the data is asumed to be absolute
	585	"""
[3180]	586	#FIXME: have the user set the tag and resolution here,
	587	# but still return the instance, just in case.
[2525]	588	[[x,y]] = self.geo_reference.change_points_geo_ref([x,y],
	589	points_geo_ref=geo_reference)
[5207]	590	region = self._addRegion(x, y)
	591	if tag is not None:
	592	region.setTag(tag)
	593	return region
[2525]	594
[4665]	595	def build_grid(self, vert_rows, vert_columns):
	596	"""
	597	Build a grid with vert_rows number of vertex rows and
	598	vert_columns number if vertex columns
	599
	600	Grid spacing of 1, the origin is the lower left hand corner.
	601
	602	FIXME(DSG-DSG) no test.
	603	"""
	604
	605	for i in range(vert_rows):
	606	for j in range(vert_columns):
	607	self.addUserVertex(j,i)
	608	self.auto_segment()
	609	self.generateMesh(mode = "Q", minAngle=20.0)
	610
[2525]	611
[3180]	612	def add_vertices(self, point_data):
	613	"""
	614	Add user vertices.
[3177]	615
[3180]	616	The point_data can be a list of (x,y) values, a numeric
	617	array or a geospatial_data instance.
	618	"""
	619	point_data = ensure_geospatial(point_data)
[3351]	620	#print "point_data",point_data
[3177]	621	# get points relative to the mesh geo_ref
[3180]	622	points = point_data.get_data_points(geo_reference=self.geo_reference)
[2525]	623
[3177]	624	for point in points:
	625	v=Vertex(point[0], point[1])
	626	self.userVertices.append(v)
	627
[4888]	628	def add_hole_from_polygon(self, polygon, segment_tags=None,
	629	geo_reference=None):
[2525]	630	"""
	631	Add a polygon with tags to the current mesh, as a region.
	632	The maxArea of the region can be specified.
	633
	634	If a geo_reference of the polygon points is given, this is used.
	635	If not;
	636	The x,y info is assumed to be Easting and Northing, absolute,
	637	for the meshes zone.
	638
	639	polygon a list of points, in meters that describe the polygon
	640	(e.g. [[x1,y1],[x2,y2],...]
	641	tags (e.g.{'wall':[0,1,3],'ocean':[2]})
	642
	643	This returns the region instance, so if the user whats to modify
	644	it they can.
	645	"""
[4888]	646	return self._add_area_from_polygon(polygon,
	647	segment_tags=segment_tags,
[2525]	648	hole=True,
	649	geo_reference=geo_reference
	650	)
	651
	652
[4888]	653	def add_region_from_polygon(self, polygon, segment_tags=None,
	654	max_triangle_area=None, geo_reference=None,
	655	region_tag=None):
[2525]	656	"""
	657	Add a polygon with tags to the current mesh, as a region.
	658	The maxArea of the region can be specified.
	659
	660	If a geo_reference of the polygon points is given, this is used.
	661	If not;
	662	The x,y info is assumed to be Easting and Northing, absolute,
	663	for the meshes zone.
	664
	665	polygon a list of points, in meters that describe the polygon
	666	(e.g. [[x1,y1],[x2,y2],...]
[4888]	667	segment_tags (e.g.{'wall':[0,1,3],'ocean':[2]}) add tags to the
	668	segements of the region
	669	region_tags - add a tag to all of the triangles in the region.
[2525]	670
	671	This returns the region instance (if a max_triangle_area is given),
	672	so if the user whats to modify it they can.
	673	"""
[5207]	674	if max_triangle_area is None and region_tag is None:
[2525]	675	create_region = False
	676	else:
	677	create_region = True
	678
[4888]	679	region = self._add_area_from_polygon(polygon,
	680	segment_tags=segment_tags,
[2525]	681	geo_reference=geo_reference,
	682	region=create_region)
	683	if max_triangle_area is not None:
	684	region.setMaxArea(max_triangle_area)
[4888]	685	if region_tag is not None:
	686	region.setTag(region_tag)
	687
	688
[2525]	689	return region
	690
	691
[4888]	692	def _add_area_from_polygon(self, polygon, segment_tags=None,
[2525]	693	geo_reference=None,
	694	hole=False,
	695	region=False):
	696	"""
	697	Add a polygon with tags to the current mesh, as a region.
	698	The maxArea of the region can be specified.
	699
	700	If a geo_reference of the polygon points is given, this is used.
	701	If not;
	702	The x,y info is assumed to be Easting and Northing, absolute,
	703	for the meshes zone.
	704
	705	polygon a list of points, in meters that describe the polygon
	706	(e.g. [[x1,y1],[x2,y2],...]
[4888]	707	segment_tags (e.g.{'wall':[0,1,3],'ocean':[2]})add tags to the
	708	segements of the region.
[2525]	709
	710	This returns the region instance, so if the user whats to modify
	711	it they can.
	712
	713	"""
[4830]	714	# Only import this if necessary.
	715	# Trying to get pmesh working in an uncompiled environment
[7711]	716	from anuga.geometry.polygon import point_in_polygon
[4830]	717
[2525]	718	#get absolute values
	719	if geo_reference is not None:
	720	polygon = geo_reference.get_absolute(polygon)
	721	# polygon is now absolute
	722	#print "polygon should be absolute",polygon
	723
	724	#create points, segs and tags
	725	region_dict = {}
	726	region_dict['points'] = polygon
	727
	728	#Create segments
	729	#E.g. [[0,1], [1,2], [2,3], [3,0]]
	730	#from polygon
	731	#[0,1,2,3]
	732	segments = []
	733	N = len(polygon)
	734	for i in range(N):
	735	lo = i
	736	hi = (lo + 1) % N
	737	segments.append( [lo, hi] )
	738	region_dict['segments'] = segments
[4888]	739	region_dict['segment_tags'] = self._tag_dict2list(segment_tags, N)
[3541]	740
[2525]	741
	742	self.addVertsSegs(region_dict) #this is passing absolute values
	743
	744	if region is True:
	745	#get inner point - absolute values
	746	inner_point = point_in_polygon(polygon)
	747	inner = self.add_region(inner_point[0], inner_point[1],
	748	geo_reference=None)
	749	elif hole is True:
	750	#get inner point - absolute values
	751	inner_point = point_in_polygon(polygon)
	752	inner = self.add_hole(inner_point[0], inner_point[1],
	753	geo_reference=None)
	754	else:
	755	inner = None
	756
	757	return inner
	758
[3541]	759	def _tag_dict2list(self, tags, number_of_segs):
	760	"""
	761	Convert a tag dictionary from this sort of format;
	762	#{'wall':[0,3],'ocean':[2]}
	763
	764	To a list format
	765	# ['wall', '', 'ocean', 'wall']
	766
	767	Note: '' is a default value of nothing
	768	"""
	769	# FIXME (DSG-DSG): Using '' as a default isn't good.
	770	#Try None.
	771	# Due to this default this method is too connected to
	772	# _add_area_from_polygon
	773	segment_tags = ['']*number_of_segs
	774	if tags is not None:
	775	for key in tags:
	776	indices = tags[key]
	777	for i in indices:
	778	segment_tags[i] = key
	779	return segment_tags
	780
[2525]	781	def add_circle(self, center, radius, segment_count=100,
	782	center_geo_reference=None, tag = None,
	783	region=False, hole=False):
	784	"""
	785	center is a point, in absulute or relative to center_geo_ref
	786	radius is the radius of the circle
	787	segment_count is the number of segments used to represent the circle
	788	tag is a string name, given to the segments.
	789	If region is True a region object is returned.
	790	If hole is True a hole object is returned.
	791	(Don't have them both True.)
	792
	793
	794	"""
	795	# convert center and radius to a polygon
	796	cuts = []
	797	factor = 2* math.pi/segment_count
	798	for cut in range(segment_count):
	799	cuts.append(cut*factor)
	800
	801	polygon = []
	802	for cut in cuts:
	803
	804	x = center[0] + radius * math.cos(cut)
	805	y = center[1] + radius * math.sin(cut)
	806	polygon.append([x,y])
	807	# build the tags
	808	tags = {tag:range(segment_count)}
	809
[4888]	810	return self._add_area_from_polygon(polygon, segment_tags=tags,
[2525]	811	region=region, hole=hole,
	812	geo_reference=center_geo_reference)
	813
	814	def auto_set_geo_reference(self):
	815	"""
	816	Automatically set the georeference, based in the minimum x and y
	817	user vertex values.
	818
	819	Not to be used with the graphical interface
	820
	821	Not implemented.
	822	Don't implement now. using the new georeferenced points class
	823	will change this?
	824	"""
	825	#to do
	826	# find the lower left hand corner
	827	[xmin, ymin, xmax, ymax] = self.boxsize()
	828
	829	# set all points to that lower left hand corner.
	830	# use change_geo_reference
	831	new_geo = Geo_reference(self.geo_reference.get_zone(), xmin, ymin)
	832	self.change_geo_reference(new_geo)
	833
	834	def change_geo_reference(self, new_geo_reference):
	835	"""
	836	Change from one geo_reference to another.
	837	Not to be used with the graphical interface
	838	"""
	839	# FIXME
	840	# change georeference of;
	841	#self.userVertices = self.geo_reference.change_points_geo_ref( \
	842	#self.userVertices)
	843	#self.holes = self.geo_reference.change_points_geo_ref(self.holes)
	844	#self.regions = self.geo_reference.change_points_geo_ref(self.regions)
	845	# The above will not work.
	846	# since userVertices (etc) is a list of point objects,
	847	#not a list of lists.
	848	# add a method to the points class to fix this up.
[3541]	849
	850	def add_segment(self, v1, v2, tag):
	851	"""
	852	Don't do this function.
	853	what will the v's be objects? How is the user suppost to get them?
	854
	855	Indexes? If add vertstosegs or add_region_from_polygon is called
	856	more than once then the actual index is not obvious. Making this
	857	function confusing.
	858	"""
	859	pass
	860
	861
	862	def add_points_and_segments(self, points,
[5673]	863	segments=None, segment_tags = None):
[3541]	864	"""
	865	Add an outline of the mesh.
[5673]	866	Points is a list of points a standard representation of points.
[3541]	867	Segments is a list of tuples of integers. Each tuple defines the
[5673]	868	start and end of the segment by it's point index.
[3541]	869	segment_tags is an optional dictionary which is used to add tags to
	870	the segments. The key is the tag name, value is the list of segment
	871	indexes the tag will apply to.
	872	eg. {'wall':[0,3],'ocean':[2]}
	873
	874	"""
[5931]	875	# make sure the points are absolute
	876	# Since addVertsSegs will deal with georeferencing.
[3541]	877	points = ensure_absolute(points)
[5673]	878
	879	if segments is None:
	880	segments = []
	881	for i in range(len(points)-1):
	882	segments.append([i, i+1])
	883
[3541]	884	#create points, segs and tags
	885	region_dict = {}
	886	region_dict['points'] = points
	887	region_dict['segments'] = segments
	888	region_dict['segment_tags'] = self._tag_dict2list(segment_tags,
	889	len(segments))
	890	self.addVertsSegs(region_dict)
	891
[2525]	892	def addVertsSegs(self, outlineDict):
	893	"""
	894	Add out-line (user Mesh) attributes given a dictionary of the lists
	895	points: [(x1,y1),(x2,y2),...] (Tuples of doubles)
	896	segments: [(point1,point2),(p3,p4),...] (Tuples of integers)
	897	segment_tags: [S1Tag, S2Tag, ...] (list of strings)
	898
	899	Assume the values are in Eastings and Northings, with no reference
	900	point. eg absolute
	901	"""
	902	if not outlineDict.has_key('segment_tags'):
	903	outlineDict['segment_tags'] = []
	904	for i in range(len(outlineDict['segments'])):
	905	outlineDict['segment_tags'].append('')
	906	#print "outlineDict['segment_tags']",outlineDict['segment_tags']
	907	#print "outlineDict['points']",outlineDict['points']
	908	#print "outlineDict['segments']",outlineDict['segments']
	909
[3541]	910	i_offset = len(self.userVertices)
	911	#print "self.userVertices",self.userVertices
	912	#print "index_offset",index_offset
[2525]	913	for point in outlineDict['points']:
	914	v=Vertex(point[0]-self.geo_reference.xllcorner,
	915	point[1]-self.geo_reference.yllcorner)
	916	self.userVertices.append(v)
	917
	918	for seg,seg_tag in map(None,outlineDict['segments'],
	919	outlineDict['segment_tags']):
[3541]	920	segObject = Segment(self.userVertices[int(seg[0])+i_offset],
	921	self.userVertices[int(seg[1])+i_offset] )
[2525]	922	if not seg_tag == '':
	923	segObject.set_tag(seg_tag)
	924	self.userSegments.append(segObject)
[3541]	925
[2525]	926
	927	def get_triangle_count(self):
	928	return len(self.getTriangulation())
	929
	930	def getUserVertices(self):
[3177]	931	"""
	932	Note: The x,y values will be relative to the mesh geo_reference
	933	This returns a list of vertex objects
	934	"""
[2525]	935	return self.userVertices
[3177]	936
	937	def get_user_vertices(self, absolute=True):
	938	"""
	939	This returns a list of (x,y) values
	940	(maybe it should return a geospatical object?
	941	It makes it a bit confusing though.)
	942	"""
	943	pointlist=[]
	944	for vertex in self.userVertices:
	945	pointlist.append([vertex.x,vertex.y])
	946	spat = Geospatial_data(pointlist, geo_reference=self.geo_reference)
	947	return spat.get_data_points(absolute=absolute)
[2525]	948
	949	def getUserSegments(self):
	950	allSegments = self.userSegments + self.alphaUserSegments
	951	#print "self.userSegments",self.userSegments
	952	#print "self.alphaUserSegments",self.alphaUserSegments
	953	#print "allSegments",allSegments
	954	return allSegments
	955
	956	def deleteUserSegments(self,seg):
	957	if self.userSegments.count(seg) == 0:
	958	self.alphaUserSegments.remove(seg)
	959	pass
	960	else:
	961	self.userSegments.remove(seg)
	962
	963	def clearUserSegments(self):
	964	self.userSegments = []
	965	self.alphaUserSegments = []
[4902]	966
	967	#FIXME see where this is used. return an array instead
[2525]	968	def getTriangulation(self):
[4902]	969	#return self.meshTriangles
	970	return self.tri_mesh.triangles.tolist()
[2525]	971
	972	def getMeshVertices(self):
[4902]	973	#return self.meshVertices
	974	return self.tri_mesh.vertices
[2525]	975
[4902]	976	def getMeshVerticeAttributes(self):
	977	#return self.meshVertices
	978	return self.tri_mesh.vertex_attributes
	979
[2525]	980	def getMeshSegments(self):
[4902]	981	#return self.meshSegments
	982	return self.tri_mesh.segments
[2525]	983
[4902]	984	def getMeshSegmentTags(self):
	985	#return self.meshSegments
	986	return self.tri_mesh.segment_tags
	987
[2525]	988	def getHoles(self):
	989	return self.holes
	990
	991	def getRegions(self):
	992	return self.regions
	993
	994	def isTriangulation(self):
	995	if self.meshVertices == []:
	996	return False
	997	else:
	998	return True
	999
	1000	def addUserSegment(self, v1,v2):
	1001	"""
	1002	PRECON: A segment between the two vertices is not already present.
	1003	Check by calling isUserSegmentNew before calling this function.
	1004
	1005	"""
	1006	s=Segment( v1,v2)
	1007	self.userSegments.append(s)
	1008	return s
	1009
	1010	def generate_mesh(self,
[3193]	1011	maximum_triangle_area="",
[2525]	1012	minimum_triangle_angle=28.0,
[3193]	1013	verbose=True):
[2525]	1014	if verbose is True:
	1015	silent = ''
	1016	else:
	1017	silent = 'Q'
	1018	self.generateMesh(mode = silent +"pzq"+str(minimum_triangle_angle)
	1019	+"a"+str(maximum_triangle_area)
	1020	+"a")
	1021	#The last a is so areas for regions will be used
	1022
	1023	def generateMesh(self, mode = None, maxArea = None, minAngle= None,
	1024	isRegionalMaxAreas = True):
	1025	"""
	1026	Based on the current user vaules, holes and regions
	1027	generate a new mesh
	1028	mode is a string that sets conditions on the mesh generations
	1029	see triangle_instructions.txt for a definition of the commands
[4831]	1030
	1031	PreCondition: maxArea is a double between 1e-20 and 1e30 or is a
	1032	string.
[2525]	1033	"""
	1034	#print "mode ",mode
	1035	if mode == None:
	1036	self.mode = ""
	1037	else:
	1038	self.mode = mode
	1039
	1040	if self.mode.find('p') < 0:
	1041	self.mode += 'p' #p - read a planar straight line graph.
	1042	#there must be segments to use this switch
	1043	# TODO throw an aception if there aren't seg's
	1044	# it's more comlex than this. eg holes
	1045	if self.mode.find('z') < 0:
	1046	self.mode += 'z' # z - Number all items starting from zero
	1047	# (rather than one)
	1048	if self.mode.find('n'):
	1049	self.mode += 'n' # n - output a list of neighboring triangles
	1050	if self.mode.find('A') < 0:
	1051	self.mode += 'A' # A - output region attribute list for triangles
	1052
	1053	if not self.mode.find('V') and not self.mode.find('Q'):
	1054	self.mode += 'V' # V - output info about what Triangle is doing
	1055
	1056	if self.mode.find('q') < 0 and minAngle is not None:
	1057	# print "********8minAngle****** ",minAngle
	1058	min_angle = 'q' + str(minAngle)
	1059	self.mode += min_angle # z - Number all items starting from zero
	1060	# (rather than one)
	1061	if maxArea != None:
	1062	self.mode += 'a' + str(maxArea)
[4831]	1063	try:
	1064	self.mode += 'a' + '%20.20f' %maxArea
	1065	except TypeError:
	1066	self.mode += 'a' + str(maxArea)
	1067	#print "self.mode", self.mode
[3193]	1068	#FIXME (DSG-DSG) This isn't explained.
[2525]	1069	if isRegionalMaxAreas:
	1070	self.mode += 'a'
[3193]	1071	#print "mesh#generateMesh# self.mode",self.mode
[2525]	1072	meshDict = self.Mesh2triangList()
[4144]	1073
	1074	#FIXME (DSG-DSG) move below section into generate_mesh.py
	1075	# & 4 functions eg segment_strings2ints
[4898]	1076	# Actually, because of region_list.append((1.0,2.0,""))
	1077	# don't move it, without careful thought
[2525]	1078	#print "*************************!@!@ This is going to triangle !@!@"
	1079	#print meshDict
	1080	#print "************************!@!@ This is going to triangle !@!@"
	1081
	1082	#print "meshDict['segmenttaglist']", meshDict['segmenttaglist']
	1083	[meshDict['segmenttaglist'],
	1084	segconverter] = segment_strings2ints(meshDict['segmenttaglist'],
	1085	initialconversions)
	1086	#print "regionlist",meshDict['regionlist']
	1087	[meshDict['regionlist'],
	1088	regionconverter] = region_strings2ints(meshDict['regionlist'])
	1089	#print "%%%%%%%%%%%%%%%%%%%%%%%%%%%regionlist",meshDict['regionlist']
	1090	#print "meshDict['segmenttaglist']", meshDict['segmenttaglist'
[4741]	1091	#print "self.mode", self.mode
[4144]	1092	generatedMesh = generate_mesh(
[2525]	1093	meshDict['pointlist'],
	1094	meshDict['segmentlist'],
	1095	meshDict['holelist'],
	1096	meshDict['regionlist'],
	1097	meshDict['pointattributelist'],
	1098	meshDict['segmenttaglist'],
[4144]	1099	self.mode,
	1100	meshDict['pointlist'])
[2525]	1101	#print "%%%%%%%%%%%%%%%%%%%%%%%%%%%generated",generatedMesh
[4668]	1102	generatedMesh['qaa'] = 1
[4899]	1103	if generatedMesh['generatedsegmentmarkerlist'] is not None:
	1104	generatedMesh['generatedsegmentmarkerlist'] = \
	1105	segment_ints2strings(generatedMesh['generatedsegmentmarkerlist'],
	1106	segconverter)
[2525]	1107	#print "processed gen",generatedMesh['generatedsegmentmarkerlist']
[4899]	1108	#print "pmesh mesh generatedMesh['generatedtriangleattributelist']", generatedMesh['generatedtriangleattributelist']
	1109	if generatedMesh['generatedtriangleattributelist'] is not None:
	1110	generatedMesh['generatedtriangleattributelist'] = \
	1111	region_ints2strings(generatedMesh['generatedtriangleattributelist'],
[2525]	1112	regionconverter)
	1113
[4899]	1114	#print "pmesh mesh generatedMesh['generatedtriangleattributelist']", generatedMesh['generatedtriangleattributelist']
[4144]	1115	#FIXME (DSG-DSG) move above section into generate_mesh.py
[4899]	1116
[4902]	1117	if generatedMesh['generatedpointattributelist'] is None or \
	1118	generatedMesh['generatedpointattributelist'].shape[1] ==0:
[2525]	1119	self.attributeTitles = []
	1120	generatedMesh['generatedpointattributetitlelist']= \
	1121	self.attributeTitles
	1122	#print "################ FROM TRIANGLE"
	1123	#print "generatedMesh",generatedMesh
	1124	#print "##################"
	1125	self.setTriangulation(generatedMesh)
	1126
	1127	def clearTriangulation(self):
	1128
	1129	#Clear the current generated mesh values
	1130	self.meshTriangles=[]
	1131	self.meshSegments=[]
	1132	self.meshVertices=[]
	1133
	1134	def removeDuplicatedUserVertices(self):
	1135	"""Pre-condition: There are no user segments
	1136	This function will keep the first duplicate
	1137	"""
	1138	assert self.getUserSegments() == []
	1139	self.userVertices, counter = self.removeDuplicatedVertices(
	1140	self.userVertices)
	1141	return counter
	1142
	1143	def removeDuplicatedVertices(self, Vertices):
	1144	"""
	1145	This function will keep the first duplicate, remove all others
	1146	Precondition: Each vertex has a dupindex, which is the list
	1147	index.
	1148
	1149	Note: this removes vertices that have the same x,y values,
	1150	not duplicate instances in the Vertices list.
	1151	"""
	1152	remove = []
	1153	index = 0
	1154	for v in Vertices:
	1155	v.dupindex = index
	1156	index += 1
	1157	t = list(Vertices)
	1158	t.sort(Point.cmp_xy)
	1159
	1160	length = len(t)
	1161	behind = 0
	1162	ahead = 1
	1163	counter = 0
	1164	while ahead < length:
	1165	b = t[behind]
	1166	ah = t[ahead]
	1167	if (b.y == ah.y and b.x == ah.x):
	1168	remove.append(ah.dupindex)
	1169	behind += 1
	1170	ahead += 1
	1171
	1172	# remove the duplicate vertices
	1173	remove.sort()
	1174	remove.reverse()
	1175	for i in remove:
	1176	Vertices.pop(i)
	1177	pass
	1178
	1179	#Remove the attribute that this function added
	1180	for v in Vertices:
	1181	del v.dupindex
	1182	return Vertices,counter
[4663]	1183
	1184	# FIXME (DSG-DSG) Move this to geospatial
[2525]	1185	def thinoutVertices(self, delta):
	1186	"""Pre-condition: There are no user segments
	1187	This function will keep the first duplicate
	1188	"""
	1189	assert self.getUserSegments() == []
	1190	#t = self.userVertices
	1191	#self.userVertices =[]
	1192	boxedVertices = {}
	1193	thinnedUserVertices =[]
	1194	delta = round(delta,1)
	1195
	1196	for v in self.userVertices :
	1197	# tag is the center of the boxes
	1198	tag = (round(v.x/delta,0)delta,round(v.y/delta,0)delta)
	1199	#this creates a dict of lists of faces, indexed by tag
	1200	boxedVertices.setdefault(tag,[]).append(v)
	1201
	1202	for [tag,verts] in boxedVertices.items():
	1203	min = delta
	1204	bestVert = None
	1205	tagVert = Vertex(tag[0],tag[1])
	1206	for v in verts:
	1207	dist = v.DistanceToPoint(tagVert)
	1208	if (dist<min):
	1209	min = dist
	1210	bestVert = v
	1211	thinnedUserVertices.append(bestVert)
	1212	self.userVertices =thinnedUserVertices
	1213
	1214
	1215	def isUserSegmentNew(self, v1,v2):
	1216	identicalSegs= [x for x in self.getUserSegments() \
	1217	if (x.vertices[0] == v1 and x.vertices[1] == v2)
	1218	or (x.vertices[0] == v2 and x.vertices[1] == v1) ]
	1219
	1220	return len(identicalSegs) == 0
	1221
	1222
	1223	def deleteSegsOfVertex(self, delVertex):
	1224	"""
	1225	Delete this vertex and any segments that connect to it.
	1226	"""
	1227	#Find segments that connect to delVertex
	1228	deletedSegments = []
	1229	for seg in self.getUserSegments():
	1230	if (delVertex in seg.vertices):
	1231	deletedSegments.append(seg)
	1232	# Delete segments that connect to delVertex
	1233	for seg in deletedSegments:
	1234	self.deleteUserSegments(seg)
	1235	return deletedSegments
	1236
	1237
	1238	def deleteMeshObject(self, MeshObject):
	1239	"""
	1240	Returns a list of all objects that were removed
	1241	"""
	1242	deletedObs = []
	1243	if isinstance(MeshObject, Vertex ):
	1244	deletedObs = self.deleteSegsOfVertex(MeshObject)
	1245	deletedObs.append(MeshObject)
	1246	self.userVertices.remove(MeshObject)
	1247	elif isinstance(MeshObject, Segment):
	1248	deletedObs.append(MeshObject)
	1249	self.deleteUserSegments(MeshObject)
	1250	elif isinstance(MeshObject, Hole):
	1251	deletedObs.append(MeshObject)
	1252	self.holes.remove(MeshObject)
	1253	elif isinstance(MeshObject, Region):
	1254	deletedObs.append(MeshObject)
	1255	self.regions.remove(MeshObject)
	1256	return deletedObs
	1257
	1258	def Mesh2triangList(self, userVertices=None,
	1259	userSegments=None,
	1260	holes=None,
	1261	regions=None):
	1262	"""
	1263	Convert the Mesh to a dictionary of the lists needed for the
	1264	triang module
	1265	points list: [(x1,y1),(x2,y2),...] (Tuples of doubles)
	1266	pointattributelist: [(a11,a12,...),(a21,a22),...] (Tuples of doubles)
	1267	segment list: [(point1,point2),(p3,p4),...] (Tuples of integers)
	1268	hole list: [(x1,y1),...](Tuples of doubles, one inside each hole)
	1269	regionlist: [ (x1,y1,tag, max area),...] (Tuple of 3-4 doubles)
	1270
	1271	Note, this adds an index attribute to the user Vertex objects.
	1272
	1273	Used to produce output to triangle
	1274	"""
	1275	if userVertices is None:
	1276	userVertices = self.getUserVertices()
	1277	if userSegments is None:
	1278	userSegments = self.getUserSegments()
	1279	if holes is None:
	1280	holes = self.getHoles()
	1281	if regions is None:
	1282	regions = self.getRegions()
	1283
	1284	meshDict = {}
	1285
	1286	pointlist=[]
	1287	pointattributelist=[]
	1288	index = 0
	1289	for vertex in userVertices:
	1290	vertex.index = index
	1291	pointlist.append((vertex.x,vertex.y))
	1292	pointattributelist.append((vertex.attributes))
	1293
	1294	index += 1
	1295	meshDict['pointlist'] = pointlist
	1296	meshDict['pointattributelist'] = pointattributelist
	1297
	1298	segmentlist=[]
	1299	segmenttaglist=[]
	1300	for seg in userSegments:
	1301	segmentlist.append((seg.vertices[0].index,seg.vertices[1].index))
	1302	segmenttaglist.append(seg.tag)
	1303	meshDict['segmentlist'] =segmentlist
	1304	meshDict['segmenttaglist'] =segmenttaglist
	1305
	1306	holelist=[]
	1307	for hole in holes:
	1308	holelist.append((hole.x,hole.y))
	1309	meshDict['holelist'] = holelist
	1310
	1311	regionlist=[]
	1312	for region in regions:
	1313	if (region.getMaxArea() != None):
	1314	regionlist.append((region.x,region.y,region.getTag(),
	1315	region.getMaxArea()))
	1316	else:
	1317	regionlist.append((region.x,region.y,region.getTag()))
	1318	meshDict['regionlist'] = regionlist
	1319	#print "(("
	1320	#print meshDict
	1321	#print meshDict['regionlist']
	1322	#print "(("
	1323	return meshDict
	1324
	1325	def Mesh2MeshList(self):
	1326	"""
	1327	Convert the Mesh to a dictionary of lists describing the
	1328	triangulation variables;
[4902]	1329
	1330	This is used by __cmp__
[2525]	1331	generated point list: [(x1,y1),(x2,y2),...] (Tuples of doubles)
	1332	generated point attribute list: [(a11,a12,...),(a21,a22),...]
	1333	(Tuples of doubles)
	1334	generated point attribute title list:[A1Title, A2Title ...]
	1335	(list of strings)
	1336	generated segment list: [(point1,point2),(p3,p4),...]
	1337	(Tuples of integers)
	1338	generated segment tag list: [tag,tag,...] list of strings
	1339
	1340	generated triangle list: [(p1,p2,p3), (p4,p5,p6),....] tuple of points
	1341
	1342	generated triangle attribute list: [s1,s2,...] list of strings
	1343
	1344	generated triangle neighbor list: [(t1,t2,t3), (t4,t5,t6),....]
	1345	tuple of triangles
	1346
	1347	Used to produce .tsh file
	1348	"""
	1349
[4902]	1350	meshDict = {}
	1351	#print "old meshDict['generatedpointattributetitlelist']",meshDict['generatedpointattributetitlelist']
	1352	#print "self.tri_mesh", self.tri_mesh
	1353	if self.tri_mesh is not None:
	1354	#print "self.tri_mesh.triangles.tolist()", self.tri_mesh.triangles.tolist()
	1355	meshDict['generatedtrianglelist'] = self.tri_mesh.triangles.tolist()
	1356	meshDict['generatedtriangleattributelist'] = self.tri_mesh.triangle_tags
	1357	meshDict['generatedtriangleneighborlist'] = self.tri_mesh.triangle_neighbors.tolist()
	1358	meshDict['generatedpointlist'] = self.tri_mesh.vertices.tolist()
	1359	if self.tri_mesh.vertex_attributes == []:
	1360	meshDict['generatedpointattributelist'] = []
	1361	#meshDict['generatedpointattributelist'] = self.tri_mesh.vertex_attributes
	1362	meshDict['generatedpointattributetitlelist'] = \
	1363	self.tri_mesh.vertex_attribute_titles
	1364	meshDict['generatedsegmentlist'] = self.tri_mesh.segments.tolist()
	1365	meshDict['generatedsegmenttaglist'] = self.tri_mesh.segment_tags
	1366	else:
	1367	meshDict['generatedtrianglelist'] = []
	1368	meshDict['generatedtriangleattributelist'] = []
	1369	meshDict['generatedtriangleneighborlist'] = []
	1370	meshDict['generatedpointlist'] = []
	1371	meshDict['generatedpointattributelist'] = []
	1372	meshDict['generatedpointattributetitlelist'] = []
	1373	meshDict['generatedsegmentlist'] = []
	1374	meshDict['generatedsegmenttaglist'] = []
	1375
	1376	#print "new meshDict['generatedpointattributetitlelist']",meshDict['generatedpointattributetitlelist']
[2525]	1377	#print "mesh.Mesh2MeshList))"
	1378	#print meshDict
	1379	#print "mesh.Mesh2MeshList))"
	1380
	1381	return meshDict
	1382
	1383
	1384	def Mesh2MeshDic(self):
	1385	"""
	1386	Convert the user and generated info of a mesh to a dictionary
	1387	structure
	1388	"""
	1389	dic = self.Mesh2triangList()
	1390	dic_mesh = self.Mesh2MeshList()
	1391	for element in dic_mesh.keys():
	1392	dic[element] = dic_mesh[element]
	1393	return dic
	1394
	1395	def setTriangulation(self, genDict):
	1396	"""
	1397	Set the mesh attributes given a dictionary of the lists
	1398	returned from the triang module
	1399	generated point list: [(x1,y1),(x2,y2),...] (Tuples of doubles)
	1400	generated point attribute list:[(P1att1,P1attt2, ...),
[4902]	1401	(P2att1,P2attt2,...),...]- not implemented
[2525]	1402	generated point attribute title list:[A1Title, A2Title ...]
[4902]	1403	(list of strings) - not implemented
[2525]	1404	generated segment list: [(point1,point2),(p3,p4),...]
	1405	(Tuples of integers)
[4902]	1406	generated segment marker list: [S1Tag, S2Tag, ...] (list of strings)
[2525]	1407	triangle list: [(point1,point2, point3),(p5,p4, p1),...]
	1408	(Tuples of integers)
	1409	triangle neighbor list: [(triangle1,triangle2, triangle3),
	1410	(t5,t4, t1),...] (Tuples of integers)
	1411	-1 means there's no triangle neighbor
[4902]	1412	triangle attribute list: [(T1att), (T2att), ...](list of strings)
[4144]	1413
[2525]	1414	"""
[4901]	1415	# Setting up the rigid triangulation
	1416	self.tri_mesh = Rigid_triangulation(
	1417	genDict['generatedtrianglelist']
	1418	,genDict['generatedsegmentlist']
	1419	,genDict['generatedpointlist']
	1420	,genDict['generatedtriangleattributelist']
	1421	,genDict['generatedtriangleneighborlist']
	1422	,genDict['generatedsegmentmarkerlist']
[4902]	1423	,genDict['generatedpointattributelist']
	1424	,genDict['generatedpointattributetitlelist']
[4901]	1425	)
[4668]	1426
[2525]	1427	def setMesh(self, genDict):
	1428	"""
	1429	Set the user Mesh attributes given a dictionary of the lists
	1430	point list: [(x1,y1),(x2,y2),...] (Tuples of doubles)
	1431	point attribute list:[(P1att1,P1attt2, ...),(P2att1,P2attt2,...),...]
	1432	segment list: [(point1,point2),(p3,p4),...] (Tuples of integers)
	1433	segment tag list: [S1Tag, S2Tag, ...] (list of ints)
	1434	region list: [(x1,y1),(x2,y2),...] (Tuples of doubles)
	1435	region attribute list: ["","reservoir",""] list of strings
	1436	region max area list:[real, None, Real,...] list of None and reals
	1437
	1438	mesh is an instance of a mesh object
	1439	"""
	1440	#Clear the current user mesh values
	1441	self.clearUserSegments()
	1442	self.userVertices=[]
	1443	self.Holes=[]
	1444	self.Regions=[]
	1445
	1446	#print "mesh.setMesh@#@#@#"
	1447	#print genDict
	1448	#print "@#@#@#"
	1449
	1450	#index = 0
	1451	for point in genDict['pointlist']:
	1452	v=Vertex(point[0], point[1])
	1453	#v.index = index
	1454	#index +=1
	1455	self.userVertices.append(v)
	1456
	1457	#index = 0
	1458	for seg,tag in map(None,genDict['segmentlist'],
	1459	genDict['segmenttaglist']):
[2778]	1460	segObject = Segment( self.userVertices[int(seg[0])],
	1461	self.userVertices[int(seg[1])], tag = tag )
[2525]	1462	#segObject.index = index
	1463	#index +=1
	1464	self.userSegments.append(segObject)
	1465
	1466	# Remove the loading of attribute info.
	1467	# Have attribute info added using least_squares in pyvolution
	1468	# index = 0
	1469	# for att in genDict['pointattributelist']:
	1470	# if att == None:
	1471	# self.userVertices[index].setAttributes([])
	1472	# else:
	1473	# self.userVertices[index].setAttributes(att)
	1474	# index += 1
	1475
	1476	#index = 0
	1477	for point in genDict['holelist']:
	1478	h=Hole(point[0], point[1])
	1479	#h.index = index
	1480	#index +=1
	1481	self.holes.append(h)
	1482
	1483	#index = 0
	1484	for reg,att,maxArea in map(None,
	1485	genDict['regionlist'],
	1486	genDict['regionattributelist'],
	1487	genDict['regionmaxarealist']):
	1488	Object = Region( reg[0],
	1489	reg[1],
	1490	tag = att,
	1491	maxArea = maxArea)
	1492	#Object.index = index
	1493	#index +=1
	1494	self.regions.append(Object)
	1495
[3437]	1496	def Testauto_segment(self):
[2525]	1497	newsegs = []
	1498	s1 = Segment(self.userVertices[0],
	1499	self.userVertices[1])
	1500	s2 = Segment(self.userVertices[0],
	1501	self.userVertices[2])
	1502	s3 = Segment(self.userVertices[2],
	1503	self.userVertices[1])
	1504	if self.isUserSegmentNew(s1.vertices[0],s1.vertices[1]):
	1505	newsegs.append(s1)
	1506	if self.isUserSegmentNew(s2.vertices[0],s2.vertices[1]):
	1507	newsegs.append(s2)
	1508	if self.isUserSegmentNew(s3.vertices[0],s3.vertices[1]):
	1509	newsegs.append(s3)
	1510	#DSG!!!
	1511	self.userSegments.extend(newsegs)
	1512	return newsegs
	1513
	1514
	1515	def savePickle(self, currentName):
	1516	fd = open(currentName, 'w')
	1517	pickle.dump(self,fd)
	1518	fd.close()
	1519
[3437]	1520	def auto_segmentFilter(self,raw_boundary=True,
[2525]	1521	remove_holes=False,
	1522	smooth_indents=False,
	1523	expand_pinch=False):
	1524	"""
[3435]	1525	Change the filters applied on the alpha shape boundary
[2525]	1526	"""
[3435]	1527	if self.shape is None:
	1528	return [],[],0.0
[2525]	1529	return self._boundary2mesh(raw_boundary=raw_boundary,
	1530	remove_holes=remove_holes,
	1531	smooth_indents=smooth_indents,
	1532	expand_pinch=expand_pinch)
	1533
	1534
	1535
[3437]	1536	def auto_segment(self, alpha = None,
[2525]	1537	raw_boundary=True,
	1538	remove_holes=False,
	1539	smooth_indents=False,
	1540	expand_pinch=False):
	1541	"""
	1542	Precon: There must be 3 or more vertices in the userVertices structure
[5188]	1543
	1544	This returns alpha_segs_no_user_segs, segs2delete, optimum_alpha
[2525]	1545	"""
	1546	self._createBoundary(alpha=alpha)
	1547	return self._boundary2mesh(raw_boundary=raw_boundary,
	1548	remove_holes=remove_holes,
	1549	smooth_indents=smooth_indents,
	1550	expand_pinch=expand_pinch)
	1551
	1552	def _createBoundary(self,alpha=None):
	1553	"""
	1554	"""
	1555	points=[]
	1556	for vertex in self.getUserVertices():
	1557	points.append((vertex.x,vertex.y))
[4665]	1558	self.shape = anuga.alpha_shape.alpha_shape.Alpha_Shape(points,
	1559	alpha=alpha)
[2525]	1560
	1561
	1562	def _boundary2mesh(self, raw_boundary=True,
	1563	remove_holes=False,
	1564	smooth_indents=False,
	1565	expand_pinch=False):
	1566	"""
	1567	Precon there must be a shape object.
	1568	"""
	1569	self.shape.set_boundary_type(raw_boundary=raw_boundary,
	1570	remove_holes=remove_holes,
	1571	smooth_indents=smooth_indents,
	1572	expand_pinch=expand_pinch)
	1573	boundary_segs = self.shape.get_boundary()
	1574	#print "boundary_segs",boundary_segs
	1575	segs2delete = self.alphaUserSegments
	1576	#FIXME(DSG-DSG) this algorithm needs comments
	1577	new_segs = {}
	1578	#alpha_segs = []
	1579	#user_segs = []
	1580	for seg in boundary_segs:
	1581	v1 = self.userVertices[int(seg[0])]
	1582	v2 = self.userVertices[int(seg[1])]
	1583	boundary_seg = Segment(v1, v2)
	1584	new_segs[(v1,v2)] = boundary_seg
	1585
	1586	for user_seg in self.userSegments:
	1587	if new_segs.has_key((user_seg.vertices[0],
	1588	user_seg.vertices[1])):
	1589	del new_segs[user_seg.vertices[0],
	1590	user_seg.vertices[1]]
	1591	elif new_segs.has_key((user_seg.vertices[1],
	1592	user_seg.vertices[0])):
	1593	del new_segs[user_seg.vertices[1],
	1594	user_seg.vertices[0]]
	1595
	1596	optimum_alpha = self.shape.get_alpha()
	1597	alpha_segs_no_user_segs = new_segs.values()
	1598	self.alphaUserSegments = alpha_segs_no_user_segs
	1599	return alpha_segs_no_user_segs, segs2delete, optimum_alpha
	1600
	1601
	1602	def representedAlphaUserSegment(self, v1,v2):
	1603	identicalSegs= [x for x in self.alphaUserSegments \
	1604	if (x.vertices[0] == v1 and x.vertices[1] == v2)
	1605	or (x.vertices[0] == v2 and x.vertices[1] == v1) ]
	1606
	1607	if identicalSegs == []:
	1608	return None
	1609	else:
	1610	# Only return the first one.
	1611	return identicalSegs[0]
	1612
	1613	def representedUserSegment(self, v1,v2):
	1614	identicalSegs= [x for x in self.userSegments \
	1615	if (x.vertices[0] == v1 and x.vertices[1] == v2)
	1616	or (x.vertices[0] == v2 and x.vertices[1] == v1) ]
	1617
	1618	if identicalSegs == []:
	1619	return None
	1620	else:
	1621	# Only return the first one.
	1622	return identicalSegs[0]
	1623
	1624	def joinVertices(self):
	1625	"""
	1626	Return list of segments connecting all userVertices
	1627	in the order they were given
	1628
	1629	Precon: There must be 3 or more vertices in the userVertices structure
	1630	"""
	1631
	1632	newsegs = []
	1633
	1634	v1 = self.userVertices[0]
	1635	for v2 in self.userVertices[1:]:
	1636	if self.isUserSegmentNew(v1,v2):
	1637	newseg = Segment(v1, v2)
	1638	newsegs.append(newseg)
	1639	v1 = v2
	1640
	1641	#Connect last point to the first
	1642	v2 = self.userVertices[0]
	1643	if self.isUserSegmentNew(v1,v2):
	1644	newseg = Segment(v1, v2)
	1645	newsegs.append(newseg)
	1646
	1647
	1648	#Update list of user segments
	1649	#DSG!!!
	1650	self.userSegments.extend(newsegs)
	1651	return newsegs
	1652
	1653	def normaliseMesh(self,scale, offset, height_scale):
	1654	[xmin, ymin, xmax, ymax] = self.boxsize()
	1655	[attmin0, attmax0] = self.maxMinVertAtt(0)
	1656	#print "[attmin0, attmax0]" ,[attmin0, attmax0]
	1657	[attmin1, attmax1] = self.maxMinVertAtt(1)
	1658	#print [xmin, ymin, xmax, ymax]
	1659	xrange = xmax - xmin
	1660	yrange = ymax - ymin
	1661	if xrange > yrange:
	1662	min,max = xmin, xmax
	1663	else:
	1664	min,max = ymin, ymax
	1665
	1666	for obj in self.getUserVertices():
	1667	obj.x = (obj.x - xmin)/(max- min)*scale + offset
	1668	obj.y = (obj.y - ymin)/(max- min)*scale + offset
	1669	if len(obj.attributes) > 0 and attmin0 != attmax0:
	1670	obj.attributes[0] = (obj.attributes[0]-attmin0)/ \
	1671	(attmax0-attmin0)*height_scale
	1672	if len(obj.attributes) > 1 and attmin1 != attmax1:
	1673	obj.attributes[1] = (obj.attributes[1]-attmin1)/ \
	1674	(attmax1-attmin1)*height_scale
	1675
	1676	for obj in self.getMeshVertices():
	1677	obj.x = (obj.x - xmin)/(max- min)*scale + offset
	1678	obj.y = (obj.y - ymin)/(max- min)*scale + offset
	1679	if len(obj.attributes) > 0 and attmin0 != attmax0:
	1680	obj.attributes[0] = (obj.attributes[0]-attmin0)/ \
	1681	(attmax0-attmin0)*height_scale
	1682	if len(obj.attributes) > 1 and attmin1 != attmax1:
	1683	obj.attributes[1] = (obj.attributes[1]-attmin1)/ \
	1684	(attmax1-attmin1)*height_scale
	1685
	1686	for obj in self.getHoles():
	1687	obj.x = (obj.x - xmin)/(max- min)*scale + offset
	1688	obj.y = (obj.y - ymin)/(max- min)*scale + offset
	1689	for obj in self.getRegions():
	1690	obj.x = (obj.x - xmin)/(max- min)*scale + offset
	1691	obj.y = (obj.y - ymin)/(max- min)*scale + offset
	1692	[xmin, ymin, xmax, ymax] = self.boxsize()
	1693	#print [xmin, ymin, xmax, ymax]
	1694
	1695	def boxsizeVerts(self):
	1696	"""
	1697	Returns a list of verts denoting a box or triangle that contains
	1698	verts on the xmin, ymin, xmax and ymax axis.
	1699	Structure: list of verts
	1700	"""
[3193]	1701
	1702	large = kinds.default_float_kind.MAX
[2525]	1703	xmin= large
	1704	xmax=-large
	1705	ymin= large
	1706	ymax=-large
	1707	for vertex in self.userVertices:
	1708	if vertex.x < xmin:
	1709	xmin = vertex.x
	1710	xminVert = vertex
	1711	if vertex.x > xmax:
	1712	xmax = vertex.x
	1713	xmaxVert = vertex
	1714
	1715	if vertex.y < ymin:
	1716	ymin = vertex.y
	1717	yminVert = vertex
	1718	if vertex.y > ymax:
	1719	ymax = vertex.y
	1720	ymaxVert = vertex
	1721	verts, count = self.removeDuplicatedVertices([xminVert,
	1722	xmaxVert,
	1723	yminVert,
	1724	ymaxVert])
	1725
	1726	return verts
	1727
	1728	def boxsize(self):
	1729	"""
	1730	Returns a list denoting a box that contains the entire
	1731	structure of vertices
	1732	Structure: [xmin, ymin, xmax, ymax]
	1733	"""
[3193]	1734
	1735	large = kinds.default_float_kind.MAX
[2525]	1736	xmin= large
	1737	xmax=-large
	1738	ymin= large
	1739	ymax=-large
	1740	for vertex in self.userVertices:
	1741	if vertex.x < xmin:
	1742	xmin = vertex.x
	1743	if vertex.x > xmax:
	1744	xmax = vertex.x
	1745
	1746	if vertex.y < ymin:
	1747	ymin = vertex.y
	1748	if vertex.y > ymax:
	1749	ymax = vertex.y
	1750	return [xmin, ymin, xmax, ymax]
	1751
	1752	def maxMinVertAtt(self, iatt):
	1753	"""
	1754	Returns a list denoting a box that contains the entire structure
	1755	of vertices
	1756	Structure: [xmin, ymin, xmax, ymax]
	1757	"""
[3193]	1758
	1759	large = kinds.default_float_kind.MAX
[2525]	1760	min= large
	1761	max=-large
	1762	for vertex in self.userVertices:
	1763	if len(vertex.attributes) > iatt:
	1764	if vertex.attributes[iatt] < min:
	1765	min = vertex.attributes[iatt]
	1766	if vertex.attributes[iatt] > max:
	1767	max = vertex.attributes[iatt]
	1768	for vertex in self.meshVertices:
	1769	if len(vertex.attributes) > iatt:
	1770	if vertex.attributes[iatt] < min:
	1771	min = vertex.attributes[iatt]
	1772	if vertex.attributes[iatt] > max:
	1773	max = vertex.attributes[iatt]
	1774	return [min, max]
	1775
	1776	def scaleoffset(self, WIDTH, HEIGHT):
	1777	"""
	1778	Returns a list denoting the scale and offset terms that need to be
	1779	applied when converting mesh co-ordinates onto grid co-ordinates
	1780	Structure: [scale, xoffset, yoffset]
	1781	"""
	1782	OFFSET = 0.05*min([WIDTH, HEIGHT])
	1783	[xmin, ymin, xmax, ymax] = self.boxsize()
	1784	SCALE = min([0.9WIDTH, 0.9HEIGHT])/max([xmax-xmin, ymax-ymin])
	1785
	1786	if SCALE*xmin < OFFSET:
	1787	xoffset = abs(SCALE*xmin) + OFFSET
	1788	if SCALE*xmax > WIDTH - OFFSET:
	1789	xoffset= -(SCALE*xmax - WIDTH + OFFSET)
	1790	if SCALE*ymin < OFFSET:
	1791	b = abs(SCALE*ymin)+OFFSET
	1792	if SCALE*ymax > HEIGHT-OFFSET:
	1793	b = -(SCALE*ymax - HEIGHT + OFFSET)
	1794	yoffset = HEIGHT - b
	1795	return [SCALE, xoffset, yoffset]
	1796
	1797
	1798	def exportASCIIobj(self,ofile):
	1799	"""
	1800	export a file, ofile, with the format
	1801	lines: v <x> <y> <first attribute>
	1802	f <vertex #> <vertex #> <vertex #> (of the triangles)
	1803	"""
	1804	fd = open(ofile,'w')
	1805	self.writeASCIIobj(fd)
	1806	fd.close()
	1807
	1808
	1809	def writeASCIIobj(self,fd):
	1810	fd.write(" # Triangulation as an obj file\n")
	1811	numVert = str(len(self.meshVertices))
	1812
	1813	index1 = 1
	1814	for vert in self.meshVertices:
	1815	vert.index1 = index1
	1816	index1 += 1
	1817
	1818	fd.write("v "
	1819	+ str(vert.x) + " "
	1820	+ str(vert.y) + " "
	1821	+ str(vert.attributes[0]) + "\n")
	1822
	1823	for tri in self.meshTriangles:
	1824	fd.write("f "
	1825	+ str(tri.vertices[0].index1) + " "
	1826	+ str(tri.vertices[1].index1) + " "
	1827	+ str(tri.vertices[2].index1) + "\n")
	1828
	1829	def exportASCIIsegmentoutlinefile(self,ofile):
	1830	"""Write the boundary user mesh info, eg
	1831	vertices that are connected to segments,
	1832	segments
	1833	"""
	1834
	1835	verts = {}
	1836	for seg in self.getUserSegments():
	1837	verts[seg.vertices[0]] = seg.vertices[0]
	1838	verts[seg.vertices[1]] = seg.vertices[1]
	1839	meshDict = self.Mesh2IOOutlineDict(userVertices=verts.values())
[4182]	1840	export_mesh_file(ofile,meshDict)
[2525]	1841
	1842	# exportASCIImeshfile - this function is used
	1843	def export_mesh_file(self,ofile):
	1844	"""
	1845	export a file, ofile, with the format
	1846	"""
	1847
	1848	dict = self.Mesh2IODict()
[4182]	1849	export_mesh_file(ofile,dict)
[2525]	1850
[3013]	1851	# FIXME(DSG-DSG):Break this into two functions.
	1852	#One for the outline points.
	1853	#One for the mesh points.
	1854	# Note: this function is not in the gui
[2525]	1855	def exportPointsFile(self,ofile):
	1856	"""
[4663]	1857	export a points file, ofile.
[2525]	1858
	1859	"""
	1860
	1861	mesh_dict = self.Mesh2IODict()
[4165]	1862	#point_dict = {}
	1863	#point_dict['attributelist'] = {} #this will need to be expanded..
[2525]	1864	# if attributes are brought back in.
[4165]	1865	#point_dict['geo_reference'] = self.geo_reference
[2525]	1866	if mesh_dict['vertices'] == []:
[4165]	1867	#point_dict['pointlist'] = mesh_dict['points']
	1868	geo = Geospatial_data(mesh_dict['points'],
	1869	geo_reference=self.geo_reference)
[2525]	1870	else:
[4165]	1871	#point_dict['pointlist'] = mesh_dict['vertices']
	1872	geo = Geospatial_data(mesh_dict['vertices'],
	1873	geo_reference=self.geo_reference)
[2525]	1874
[4165]	1875	geo.export_points_file(ofile, absolute=True)
	1876
[2525]	1877
[3100]	1878
[5207]	1879	def import_ungenerate_file(self,ofile, tag=None, region_tag=None):
[3100]	1880	"""
	1881	Imports an ungenerate file, from arcGIS into mesh.
[5209]	1882	This file describes many polygons.
[3100]	1883
	1884	ofile is the name of the ungenerated file.
[5207]	1885	Tag is a string name to be taggged on each segment.
	1886
[5209]	1887	region_tag is the tag applied to the polygon regions.
[5207]	1888	if it is a string the one value will be assigned to all regions
	1889	if it is a list the first value in the list will be applied to the first polygon etc.
	1890	WARNING: size of list and number of polygons isn't checked
[3100]	1891
	1892	WARNING values are assumed to be absolute.
	1893	geo-refs are not taken into account..
	1894	"""
	1895
[7873]	1896	dict = load_ungenerate(ofile)
[3100]	1897	default_tag = Segment.get_default_tag()
	1898	if tag is not None:
	1899	Segment.set_default_tag(str(tag))
[5207]	1900
	1901	if region_tag is None:
	1902	self.addVertsSegs(dict)
	1903	else:
	1904	if not isinstance(region_tag, list):
	1905	region_tag = [region_tag]*len(dict['polygons'])
	1906	for a_tag,polygon in map(None, region_tag, dict['polygons']):
	1907	segment_tags = {tag:range(len(polygon))}
	1908	self.add_region_from_polygon(polygon,segment_tags=segment_tags,
	1909	region_tag=a_tag)
	1910
	1911
[3100]	1912	Segment.set_default_tag(default_tag)
	1913
[5207]	1914	# change the tag back to it's default
[3100]	1915
[2525]	1916
	1917	########### IO CONVERTERS ##################
	1918	"""
	1919	The dict fromat for IO with .tsh files is;
	1920	(the triangulation)
	1921	vertices: [[x1,y1],[x2,y2],...] (lists of doubles)
	1922	vertex_attributes: [[a11,a12,...],[a21,a22],...] (lists of doubles)
	1923	vertex_attribute_titles:[A1Title, A2Title ...] (A list of strings)
	1924	segments: [[v1,v2],[v3,v4],...] (lists of integers)
	1925	segment_tags : [tag,tag,...] list of strings
	1926	triangles : [(v1,v2,v3), (v4,v5,v6),....] lists of points
	1927	triangle tags: [s1,s2,...] A list of strings
	1928	triangle neighbors: [[t1,t2,t3], [t4,t5,t6],..] lists of triangles
	1929
	1930	(the outline)
	1931	points: [[x1,y1],[x2,y2],...] (lists of doubles)
	1932	point_attributes: [[a11,a12,...],[a21,a22],...] (lists of doubles)
	1933	outline_segments: [[point1,point2],[p3,p4],...] (lists of integers)
	1934	outline_segment_tags : [tag1,tag2,...] list of strings
	1935	holes : [[x1,y1],...](List of doubles, one inside each hole region)
	1936	regions : [ [x1,y1],...] (List of 4 doubles)
	1937	region_tags : [tag1,tag2,...] (list of strings)
	1938	region_max_areas: [ma1,ma2,...] (A list of doubles)
	1939	{Convension: A -ve max area means no max area}
	1940
	1941	"""
	1942
	1943
	1944
	1945	def Mesh2IODict(self):
	1946	"""
	1947	Convert the triangulation and outline info of a mesh to a dictionary
	1948	structure
	1949	"""
	1950	dict = self.Mesh2IOTriangulationDict()
	1951	dict_mesh = self.Mesh2IOOutlineDict()
	1952	for element in dict_mesh.keys():
	1953	dict[element] = dict_mesh[element]
	1954
	1955	# add the geo reference
	1956	dict['geo_reference'] = self.geo_reference
	1957	return dict
	1958
	1959	def Mesh2IOTriangulationDict(self):
	1960	"""
	1961	Convert the Mesh to a dictionary of lists describing the
	1962	triangulation variables;
	1963
	1964	Used to produce .tsh file
	1965	"""
[4902]	1966	meshDict = {}
[4901]	1967	if self.tri_mesh is not None:
	1968	meshDict['triangles'] = self.tri_mesh.triangles
	1969	meshDict['triangle_tags'] = self.tri_mesh.triangle_tags
	1970	#print "mesh meshDict['triangle_tags']", meshDict['triangle_tags']
	1971	meshDict['triangle_neighbors'] = self.tri_mesh.triangle_neighbors
[4902]	1972	meshDict['vertices'] = self.tri_mesh.vertices
	1973	meshDict['vertex_attributes'] = self.tri_mesh.vertex_attributes
	1974	meshDict['vertex_attribute_titles'] = \
	1975	self.tri_mesh.vertex_attribute_titles
	1976	meshDict['segments'] = self.tri_mesh.segments
	1977	meshDict['segment_tags'] = self.tri_mesh.segment_tags
[4901]	1978	else:
	1979	meshDict['triangles'] = []
	1980	meshDict['triangle_tags'] = []
	1981	meshDict['triangle_neighbors'] = []
[4902]	1982	meshDict['vertices'] = []
	1983	meshDict['vertex_attributes'] = []
	1984	meshDict['vertex_attribute_titles'] = []
	1985	meshDict['segments'] = []
	1986	meshDict['segment_tags'] = []
[2525]	1987	#print "mesh.Mesh2IOTriangulationDict))"
	1988	#print meshDict
	1989	#print "mesh.Mesh2IOTriangulationDict))"
	1990
	1991	return meshDict
	1992
	1993
	1994	def Mesh2IOOutlineDict(self, userVertices=None,
	1995	userSegments=None,
	1996	holes=None,
	1997	regions=None):
	1998	"""
	1999	Convert the mesh outline to a dictionary of the lists needed for the
	2000	triang module;
	2001
	2002	Note, this adds an index attribute to the user Vertex objects.
	2003
	2004	Used to produce .tsh file and output to triangle
	2005	"""
[4182]	2006
[2525]	2007	if userVertices is None:
	2008	userVertices = self.getUserVertices()
	2009	if userSegments is None:
	2010	userSegments = self.getUserSegments()
	2011	if holes is None:
	2012	holes = self.getHoles()
	2013	if regions is None:
	2014	regions = self.getRegions()
	2015
	2016	meshDict = {}
	2017	#print "userVertices",userVertices
	2018	#print "userSegments",userSegments
	2019	pointlist=[]
	2020	pointattributelist=[]
	2021	index = 0
	2022	for vertex in userVertices:
	2023	vertex.index = index
	2024	pointlist.append([vertex.x,vertex.y])
	2025	pointattributelist.append(vertex.attributes)
	2026
	2027	index += 1
	2028	meshDict['points'] = pointlist
	2029	meshDict['point_attributes'] = pointattributelist
	2030
	2031	segmentlist=[]
	2032	segmenttaglist=[]
	2033	for seg in userSegments:
	2034	segmentlist.append([seg.vertices[0].index,seg.vertices[1].index])
	2035	segmenttaglist.append(seg.tag)
	2036	meshDict['outline_segments'] =segmentlist
	2037	meshDict['outline_segment_tags'] =segmenttaglist
	2038
	2039	holelist=[]
	2040	for hole in holes:
	2041	holelist.append([hole.x,hole.y])
	2042	meshDict['holes'] = holelist
	2043
	2044	regionlist=[]
	2045	regiontaglist = []
	2046	regionmaxarealist = []
	2047	for region in regions:
	2048	regionlist.append([region.x,region.y])
	2049	regiontaglist.append(region.getTag())
	2050
	2051	if (region.getMaxArea() != None):
	2052	regionmaxarealist.append(region.getMaxArea())
	2053	else:
[4182]	2054	regionmaxarealist.append(NOMAXAREA)
[2525]	2055	meshDict['regions'] = regionlist
	2056	meshDict['region_tags'] = regiontaglist
	2057	meshDict['region_max_areas'] = regionmaxarealist
	2058	#print "(("
	2059	#print meshDict
	2060	#print meshDict['regionlist']
	2061	#print "(("
	2062	return meshDict
	2063
	2064	def IOTriangulation2Mesh(self, genDict):
	2065	"""
	2066	Set the mesh attributes given an tsh IO dictionary
	2067	"""
	2068	#Clear the current generated mesh values
[4901]	2069	self.tri_mesh = None
	2070
	2071	self.tri_mesh = Rigid_triangulation(
	2072	genDict['triangles']
	2073	,genDict['segments']
	2074	,genDict['vertices']
	2075	,genDict['triangle_tags']
	2076	,genDict['triangle_neighbors']
	2077	,genDict['segment_tags']
[4902]	2078	,genDict['vertex_attributes']
	2079	,genDict['vertex_attribute_titles']
[4901]	2080	)
[4902]	2081	self.attributeTitles = genDict['vertex_attribute_titles']
	2082	self.maxVertexIndex = len(genDict['vertices'])
	2083	#print "self.maxVertexIndex ", self.maxVertexIndex
	2084
[2525]	2085	def IOOutline2Mesh(self, genDict):
	2086	"""
	2087	Set the outline (user Mesh attributes) given a IO tsh dictionary
	2088
	2089	mesh is an instance of a mesh object
	2090	"""
	2091	#Clear the current user mesh values
	2092	self.clearUserSegments()
	2093	self.userVertices=[]
	2094	self.Holes=[]
	2095	self.Regions=[]
	2096
	2097	#print "mesh.IOOutline2Mesh@#@#@#"
	2098	#print "genDict",genDict
	2099	#print "@#@#@#"
	2100
	2101	#index = 0
	2102	for point in genDict['points']:
	2103	v=Vertex(point[0], point[1])
	2104	#v.index = index
	2105	#index +=1
	2106	self.userVertices.append(v)
	2107
	2108	#index = 0
	2109	for seg,tag in map(None,genDict['outline_segments'],
	2110	genDict['outline_segment_tags']):
[2778]	2111
	2112	segObject = Segment( self.userVertices[int(seg[0])],
	2113	self.userVertices[int(seg[1])], tag = tag )
[2525]	2114	#segObject.index = index
	2115	#index +=1
	2116	self.userSegments.append(segObject)
	2117
	2118	# Remove the loading of attribute info.
	2119	# Have attribute info added using least_squares in pyvolution
	2120	# index = 0
	2121	# for att in genDict['point_attributes']:
	2122	# if att == None:
	2123	# self.userVertices[index].setAttributes([])
	2124	# else:
	2125	# self.userVertices[index].setAttributes(att)
	2126	# index += 1
	2127
	2128	#index = 0
	2129	for point in genDict['holes']:
	2130	h=Hole(point[0], point[1])
	2131	#h.index = index
	2132	#index +=1
	2133	self.holes.append(h)
	2134
	2135	#index = 0
	2136	for reg,att,maxArea in map(None,
	2137	genDict['regions'],
	2138	genDict['region_tags'],
	2139	genDict['region_max_areas']):
	2140	if maxArea > 0: # maybe I should ref NOMAXAREA? Prob' not though
	2141	Object = Region( reg[0],
	2142	reg[1],
	2143	tag = att,
	2144	maxArea = maxArea)
	2145	else:
	2146	Object = Region( reg[0],
	2147	reg[1],
	2148	tag = att)
	2149
	2150	#Object.index = index
	2151	#index +=1
	2152	self.regions.append(Object)
	2153
	2154
	2155	def importMeshFromFile(ofile):
[4663]	2156	"""returns a mesh object, made from a points file or .tsh/.msh file
[2525]	2157	Often raises IOError,RuntimeError
	2158	"""
	2159	newmesh = None
[4663]	2160	if (ofile[-4:]== ".pts" or ofile[-4:]== ".txt" or \
[4165]	2161	ofile[-4:]== ".csv"):
	2162	geospatial = Geospatial_data(ofile)
	2163	dict = {}
	2164	dict['points'] = geospatial.get_data_points(absolute=False)
[2525]	2165	dict['outline_segments'] = []
	2166	dict['outline_segment_tags'] = []
	2167	dict['regions'] = []
	2168	dict['region_tags'] = []
	2169	dict['region_max_areas'] = []
	2170	dict['holes'] = []
[4165]	2171	newmesh= Mesh(geo_reference = geospatial.geo_reference)
[2525]	2172	newmesh.IOOutline2Mesh(dict)
	2173	counter = newmesh.removeDuplicatedUserVertices()
	2174	if (counter >0):
[7317]	2175	log.critical("%i duplicate vertices removed from dataset" % counter)
[2525]	2176	elif (ofile[-4:]== ".tsh" or ofile[-4:]== ".msh"):
[4182]	2177	dict = import_mesh_file(ofile)
[2525]	2178	#print "********"
	2179	#print "zq mesh.dict",dict
	2180	#print "********"
	2181	newmesh= Mesh()
	2182	newmesh.IOOutline2Mesh(dict)
	2183	newmesh.IOTriangulation2Mesh(dict)
	2184	else:
	2185	raise RuntimeError
	2186
	2187	if dict.has_key('geo_reference') and not dict['geo_reference'] == None:
	2188	newmesh.geo_reference = dict['geo_reference']
	2189	return newmesh
	2190
	2191	def loadPickle(currentName):
	2192	fd = open(currentName)
	2193	mesh = pickle.load(fd)
	2194	fd.close()
	2195	return mesh
	2196
	2197	def square_outline(side_length = 1,up = "top", left = "left", right = "right",
	2198	down = "bottom", regions = False):
	2199
	2200	a = Vertex (0,0)
	2201	b = Vertex (0,side_length)
	2202	c = Vertex (side_length,0)
	2203	d = Vertex (side_length,side_length)
	2204
	2205	s2 = Segment(b,d, tag = up)
	2206	s3 = Segment(b,a, tag = left)
	2207	s4 = Segment(d,c, tag = right)
	2208	s5 = Segment(a,c, tag = down)
	2209
	2210	if regions:
	2211	e = Vertex (side_length/2,side_length/2)
	2212	s6 = Segment(a,e, tag = down + left)
	2213	s7 = Segment(b,e, tag = up + left)
	2214	s8 = Segment(c,e, tag = down + right)
	2215	s9 = Segment(d,e, tag = up + right)
	2216	r1 = Region(side_length/2,3.*side_length/4, tag = up)
	2217	r2 = Region(1.*side_length/4,side_length/2, tag = left)
	2218	r3 = Region(3.*side_length/4,side_length/2, tag = right)
	2219	r4 = Region(side_length/2,1.*side_length/4, tag = down)
	2220	mesh = Mesh(userVertices=[a,b,c,d,e],
	2221	userSegments=[s2,s3,s4,s5,s6,s7,s8,s9],
	2222	regions = [r1,r2,r3,r4])
	2223	else:
	2224	mesh = Mesh(userVertices=[a,b,c,d],
	2225	userSegments=[s2,s3,s4,s5])
	2226
	2227	return mesh
	2228
	2229
	2230
	2231	def region_strings2ints(region_list):
	2232	"""Given a list of (x_int,y_int,tag_string) lists it returns a list of
	2233	(x_int,y_int,tag_int) and a list to convert the tag_int's back to
	2234	the tag_strings
	2235	"""
	2236	# Make sure "" has an index of 0
	2237	region_list.reverse()
	2238	region_list.append((1.0,2.0,""))
	2239	region_list.reverse()
	2240	convertint2string = []
	2241	for i in xrange(len(region_list)):
	2242	convertint2string.append(region_list[i][2])
	2243	if len(region_list[i]) == 4: # there's an area value
	2244	region_list[i] = (region_list[i][0],
	2245	region_list[i][1],i,region_list[i][3])
	2246	elif len(region_list[i]) == 3: # no area value
	2247	region_list[i] = (region_list[i][0],region_list[i][1],i)
	2248	else:
[7317]	2249	log.critical("The region list has a bad size")
[2525]	2250	# raise an error ..
	2251	raise Error
	2252
	2253	#remove "" from the region_list
	2254	region_list.pop(0)
	2255
	2256	return [region_list, convertint2string]
	2257
	2258	def region_ints2strings(region_list,convertint2string):
	2259	"""Reverses the transformation of region_strings2ints
	2260	"""
[4899]	2261	#print 'region_ints2strings region_list', region_list
	2262
	2263	returned_region_list = []
	2264	# may not need (not region_list[0] == [])
	2265	# or region_list[0] == [0.0]
	2266	if (not region_list[0] == []): # or region_list[0] == [0.0]:
	2267	#print "in loop"
[2525]	2268	for i in xrange(len(region_list)):
[4899]	2269	temp = region_list[i]
[4901]	2270	returned_region_list.append(convertint2string[int(temp[0])])
[4899]	2271	return returned_region_list
[2525]	2272
	2273	def segment_ints2strings(intlist, convertint2string):
	2274	"""Reverses the transformation of segment_strings2ints """
	2275	stringlist = []
	2276	for x in intlist:
[4955]	2277	stringlist.append(convertint2string[int(x)])
[2525]	2278	return stringlist
	2279
	2280	def segment_strings2ints(stringlist, preset):
	2281	"""Given a list of strings return a list of 0 to n ints which represent
	2282	the strings and a converting list of the strings, indexed by 0 to n ints.
	2283	Also, input an initial converting list of the strings
	2284	Note, the converting list starts off with
	2285	["internal boundary", "external boundary", "internal boundary"]
	2286	example input and output
	2287	input ["a","b","a","c"],["c"]
	2288	output [[2, 1, 2, 0], ['c', 'b', 'a']]
	2289
	2290	the first element in the converting list will be
	2291	overwritten with "".
	2292	?This will become the third item in the converting list?
	2293
	2294	# note, order the initial converting list is important,
	2295	since the index = the int tag
	2296	"""
	2297	nodups = unique(stringlist)
	2298	# note, order is important, the index = the int tag
	2299	#preset = ["internal boundary", "external boundary"]
	2300	#Remove the preset tags from the list with no duplicates
	2301	nodups = [x for x in nodups if x not in preset]
	2302
	2303	try:
	2304	nodups.remove("") # this has to go to zero
	2305	except ValueError:
	2306	pass
	2307
	2308	# Add the preset list at the beginning of no duplicates
	2309	preset.reverse()
	2310	nodups.extend(preset)
	2311	nodups.reverse()
	2312
	2313
	2314	convertstring2int = {}
	2315	convertint2string = []
	2316	index = 0
	2317	for x in nodups:
	2318	convertstring2int[x] = index
	2319	convertint2string.append(x)
	2320	index += 1
	2321	convertstring2int[""] = 0
	2322
	2323	intlist = []
	2324	for x in stringlist:
	2325	intlist.append(convertstring2int[x])
	2326	return [intlist, convertint2string]
	2327
	2328
	2329	def unique(s):
	2330	"""Return a list of the elements in s, but without duplicates.
	2331
	2332	For example, unique([1,2,3,1,2,3]) is some permutation of [1,2,3],
	2333	unique("abcabc") some permutation of ["a", "b", "c"], and
	2334	unique(([1, 2], [2, 3], [1, 2])) some permutation of
	2335	[[2, 3], [1, 2]].
	2336
	2337	For best speed, all sequence elements should be hashable. Then
	2338	unique() will usually work in linear time.
	2339
	2340	If not possible, the sequence elements should enjoy a total
	2341	ordering, and if list(s).sort() doesn't raise TypeError it's
	2342	assumed that they do enjoy a total ordering. Then unique() will
	2343	usually work in O(N*log2(N)) time.
	2344
	2345	If that's not possible either, the sequence elements must support
	2346	equality-testing. Then unique() will usually work in quadratic
	2347	time.
	2348	"""
	2349
	2350	n = len(s)
	2351	if n == 0:
	2352	return []
	2353
	2354	# Try using a dict first, as that's the fastest and will usually
	2355	# work. If it doesn't work, it will usually fail quickly, so it
	2356	# usually doesn't cost much to try it. It requires that all the
	2357	# sequence elements be hashable, and support equality comparison.
	2358	u = {}
	2359	try:
	2360	for x in s:
	2361	u[x] = 1
	2362	except TypeError:
	2363	del u # move on to the next method
	2364	else:
	2365	return u.keys()
	2366
	2367	# We can't hash all the elements. Second fastest is to sort,
	2368	# which brings the equal elements together; then duplicates are
	2369	# easy to weed out in a single pass.
	2370	# NOTE: Python's list.sort() was designed to be efficient in the
	2371	# presence of many duplicate elements. This isn't true of all
	2372	# sort functions in all languages or libraries, so this approach
	2373	# is more effective in Python than it may be elsewhere.
	2374	try:
	2375	t = list(s)
	2376	t.sort()
	2377	except TypeError:
	2378	del t # move on to the next method
	2379	else:
	2380	assert n > 0
	2381	last = t[0]
	2382	lasti = i = 1
	2383	while i < n:
	2384	if t[i] != last:
	2385	t[lasti] = last = t[i]
	2386	lasti += 1
	2387	i += 1
	2388	return t[:lasti]
	2389
	2390	# Brute force is all that's left.
	2391	u = []
	2392	for x in s:
	2393	if x not in u:
	2394	u.append(x)
	2395	return u
	2396
	2397
	2398	# FIXME (DSG-DSG)
	2399	# To do-
	2400	# Create a clear interface. eg
	2401	# have the interface methods more at the top of this file and add comments
	2402	# for the interface functions/methods, use function_name (not functionName),
	2403
	2404	#Currently
	2405	#function_name methods assume absolute values. Geo-refs can be passed in.
	2406	#
	2407
	2408	# instead of functionName
	2409	if __name__ == "__main__":
[4901]	2410	pass

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

Download in other formats: