#!/usr/bin/env python # """General 2D triangular classes for triangular mesh generation. Note: A .index attribute is added to objects such as vertices and segments, often when reading and writing to files. This information should not be used as persistant information. It is not the 'index' of an element in a list. Copyright 2003/2004 Ole Nielsen, Stephen Roberts, Duncan Gray, Christopher Zoppou Geoscience Australia """ import sys import math import re import os import pickle import types import exceptions import numpy as num from anuga.coordinate_transforms.geo_reference import Geo_reference, \ DEFAULT_ZONE from anuga.load_mesh.loadASCII import NOMAXAREA, export_mesh_file, \ import_mesh_file import anuga.alpha_shape.alpha_shape from anuga.geospatial_data.geospatial_data import Geospatial_data, \ ensure_geospatial, ensure_absolute, ensure_numeric from anuga.mesh_engine.mesh_engine import generate_mesh import anuga.utilities.log as log try: import kinds except ImportError: # Hand-built mockup of the things we need from the kinds package, since it # was recently removed from the standard numeric distro. Some users may # not have it by default. class _bunch: pass class _kinds(_bunch): default_float_kind = _bunch() default_float_kind.MIN = 2.2250738585072014e-308 #smallest +ve number default_float_kind.MAX = 1.7976931348623157e+308 kinds = _kinds() # 1st and third values must be the same # FIXME: maybe make this a switch that the user can change? - DSG initialconversions = ['', 'exterior', ''] SEG_COLOUR = 'blue' class MeshObject: """ An abstract superclass for the basic mesh objects, eg vertex, segment, triangle. """ def __init__(self): pass class Point(MeshObject): """ Define a point in a 2D space. """ def __init__(self,X,Y): __slots__ = ['x','y'] self.x=X self.y=Y def DistanceToPoint(self, OtherPoint): """ Returns the distance from this point to another """ SumOfSquares = ((self.x - OtherPoint.x)**2) + \ ((self.y - OtherPoint.y)**2) return math.sqrt(SumOfSquares) def IsInsideCircle(self, Center, Radius): """ Return 1 if this point is inside the circle, 0 otherwise """ if (self.DistanceToPoint(Center) point.x: return 1 else: if self.y < point.y: return -1 elif self.y > point.y: return 1 else: return 0 def same_x_y(self, point): if self.x == point.x and self.y == point.y: return True else: return False class Vertex(Point): """A point on the mesh. Object attributes based on the Triangle program """ VERTEXSQUARESIDELENGTH = 6 def __init__(self, X, Y, attributes=None): __slots__ = ['x','y','attributes'] # we don't care what we get, as long as we can get float *value* self.x = float(X) self.y = float(Y) self.attributes = [] if not attributes is None: self.attributes = attributes def setAttributes(self,attributes): """attributes is a list. """ self.attributes = attributes def draw(self, canvas, tags, colour = 'black',scale = 1, xoffset = 0, yoffset =0, ): x = scale*(self.x + xoffset) y = -1*scale*(self.y + yoffset) # - since for a canvas - is up cornerOffset= self.VERTEXSQUARESIDELENGTH/2 # A hack to see the vert tags # note: there will be many tags, since tags will not be removed #when zooming #canvas.create_text(x+ 2*cornerOffset, # y+ 2*cornerOffset, # text=tags) # This gives points info. It is a mess though, since numbers are # not deleted when zaooming in. #canvas.create_text(x+ 2*cornerOffset, # y+ 2*cornerOffset, # text=str(x)+','+str(y)) return canvas.create_rectangle(x-cornerOffset, y-cornerOffset, x+cornerOffset, y+cornerOffset, tags = tags, outline=colour, fill = 'white') #return tags def __repr__(self): return "[(%f,%f),%r]" % (self.x,self.y,self.attributes) class Hole(Point): """A region of the mesh were no triangles are generated. Defined by a point in the hole enclosed by segments. """ HOLECORNERLENGTH = 6 def draw(self, canvas, tags, colour = 'purple',scale = 1, xoffset = 0, yoffset =0, ): x = scale*(self.x + xoffset) y = -1*scale*(self.y + yoffset) # - since for a canvas - is up cornerOffset= self.HOLECORNERLENGTH/2 return canvas.create_oval(x-cornerOffset, y-cornerOffset, x+cornerOffset, y+cornerOffset, tags = tags, outline=colour, fill = 'white') class Region(Point): """A region of the mesh. Defined by a point in the region enclosed by segments. Used to tag areas. """ CROSSLENGTH = 6 TAG = 0 MAXAREA = 1 def __init__(self, X, Y, tag=None, maxArea=None): """Precondition: tag is a string and maxArea is a real""" self.x=X self.y=Y self.attributes = [] # index 0 is the tag string # optional index 1 is the max triangle area # NOTE the size of this attribute is assumed # to be 1 or 2 in regionstrings2int if tag is None: self.attributes.append("") else: self.attributes.append(tag) #this is a string if maxArea is not None: self.setMaxArea(maxArea) # maxArea is a number def getTag(self,): return self.attributes[self.TAG] def setTag(self,tag): self.attributes[self.TAG] = tag def getMaxArea(self): """ Returns the Max Area of a Triangle or None, if the Max Area has not been set. """ if self.isMaxArea(): return self.attributes[1] else: return None def setMaxArea(self,MaxArea): if MaxArea is not None: if self.isMaxArea(): self.attributes[self.MAXAREA] = float(MaxArea) else: self.attributes.append( float(MaxArea) ) def deleteMaxArea(self): if self.isMaxArea(): self.attributes.pop(self.MAXAREA) def isMaxArea(self): return len(self.attributes)> 1 def draw(self, canvas, tags, scale=1, xoffset = 0, yoffset =0, colour = "black"): """ Draw a black cross, returning the objectID """ x = scale*(self.x + xoffset) y = -1*scale*(self.y + yoffset) cornerOffset= self.CROSSLENGTH/2 return canvas.create_polygon(x, y-cornerOffset, x, y, x+cornerOffset, y, x, y, x, y+cornerOffset, x, y, x-cornerOffset, y, x, y, tags = tags, outline = colour,fill = '') def __repr__(self): if self.isMaxArea(): area = self.getMaxArea() return "(%f,%f,%s,%f)" % (self.x,self.y, self.getTag(), self.getMaxArea()) else: return "(%f,%f,%s)" % (self.x,self.y, self.getTag()) class Segment(MeshObject): """ Segments are edges whose presence in the triangulation is enforced. """ def __init__(self, vertex1, vertex2, tag = None ): """ Each segment is specified by listing the vertices of its endpoints The vertices are Vertex objects """ assert(vertex1 != vertex2) self.vertices = [vertex1,vertex2 ] if tag is None: self.tag = self.__class__.default else: self.tag = tag #this is a string def __repr__(self): return "[%s,%s]" % (self.vertices,self.tag) def draw(self, canvas, tags,scale=1, xoffset=0, yoffset=0, colour=SEG_COLOUR): x=[] y=[] for end in self.vertices: #end.draw(canvas,scale, xoffset, yoffset ) # draw the vertices x.append(scale*(end.x + xoffset)) y.append(-1*scale*(end.y + yoffset)) # - since for a canvas - is up return canvas.create_line(x[0], y[0], x[1], y[1], tags = tags,fill=colour) def set_tag(self,tag): self.tag = tag # Class methods def set_default_tag(cls, default): cls.default = default def get_default_tag(cls): return cls.default set_default_tag = classmethod(set_default_tag) get_default_tag = classmethod(get_default_tag) Segment.set_default_tag("") class Rigid_triangulation: """ This is a triangulation that can't have triangles added or taken away. It just represents the triangulation, not the mesh outline needed to build the triangulation. This is the guts of the data structure; generated vertex list: [(x1,y1),(x2,y2),...] (Tuples of doubles) generated segment list: [(point1,point2),(p3,p4),...] (Tuples of integers) generated segment tag list: [tag,tag,...] list of strings generated triangle list: [(p1,p2,p3), (p4,p5,p6),....] tuple of points generated triangle attribute list: [s1,s2,...] list of strings generated triangle neighbor list: [(t1,t2,t3), (t4,t5,t6),....] tuple of triangles Should I do the basic structure like triangle, general mesh or something else? How about like triangle, since that's where the info is from, and it'll fit easier into this file.. Removing the loners is difficult, since all the vert's after it must be removed. This happens in set_triangulation. """ def __init__(self, triangles, segments, vertices, triangle_tags, triangle_neighbors, segment_tags, vertex_attributes, vertex_attribute_titles=None ): self.triangles = ensure_numeric(triangles) self.triangle_neighbors = ensure_numeric(triangle_neighbors) self.triangle_tags = triangle_tags # list of strings self.segments = ensure_numeric(segments) self.segment_tags = segment_tags # list of strings self.vertices = ensure_numeric(vertices) # This is needed for __cmp__ if vertex_attributes is None: self.vertex_attributes = [] else: self.vertex_attributes = ensure_numeric(vertex_attributes) if vertex_attribute_titles is None: self.vertex_attribute_titles = [] else: self.vertex_attribute_titles = vertex_attribute_titles def draw_triangulation(self, canvas, scale=1, xoffset=0, yoffset=0, colour="green"): """ Draw a triangle, returning the objectID """ # FIXME(DSG-DSG) This could be a data structure that is # remembered and doesn't have any duplicates. # but I wouldn't be able to use create_polygon. # regard it as drawing a heap of segments for tri in self.triangles: vertices = [] for v_index in range(3): vertices.append(self.vertices[tri[v_index]]) objectID = canvas.create_polygon( scale*(vertices[1][0] + xoffset), scale*-1*(vertices[1][1] + yoffset), scale*(vertices[0][0] + xoffset), scale*-1*(vertices[0][1] + yoffset), scale*(vertices[2][0] + xoffset), scale*-1*(vertices[2][1] + yoffset), outline = colour,fill = '' ) def calc_mesh_area(self): area = 0 for tri in self.triangles: vertices = [] for v_index in range(3): vertices.append(self.vertices[tri[v_index]]) ax = vertices[0][0] ay = vertices[0][1] bx = vertices[1][0] by = vertices[1][1] cx = vertices[2][0] cy = vertices[2][1] area += abs((bx*ay-ax*by)+(cx*by-bx*cy)+(ax*cy-cx*ay))/2 return area class Mesh: """ Representation of a 2D triangular mesh. User attributes describe the mesh region/segments/vertices/attributes mesh attributes describe the mesh that is produced eg triangles and vertices. All point information is relative to the geo_reference passed in """ def __repr__(self): return """ mesh Triangles: %s mesh Attribute Titles: %s mesh Segments: %s mesh Vertices: %s user Segments: %s user Vertices: %s holes: %s regions: %s""" % (self.meshTriangles, self.attributeTitles, self.meshSegments, self.meshVertices, self.getUserSegments(), self.userVertices, self.holes, self.regions) def __init__(self, userSegments=None, userVertices=None, holes=None, regions=None, geo_reference=None): self.meshTriangles=[] self.attributeTitles=[] self.meshSegments=[] self.meshVertices=[] # Rigid self.tri_mesh=None self.visualise_graph = True if userSegments is None: self.userSegments=[] else: self.userSegments=userSegments self.alphaUserSegments=[] if userVertices is None: self.userVertices=[] else: self.userVertices=userVertices if holes is None: self.holes=[] else: self.holes=holes if regions is None: self.regions=[] else: self.regions=regions if geo_reference is None: self.geo_reference = Geo_reference(DEFAULT_ZONE,0,0) else: self.geo_reference = geo_reference self.shape = None def __cmp__(self,other): # A dic for the initial m dic = self.Mesh2triangList() dic_mesh = self.Mesh2MeshList() for element in dic_mesh.keys(): dic[element] = dic_mesh[element] for element in dic.keys(): dic[element].sort() # A dic for the exported/imported m dic_other = other.Mesh2triangList() dic_mesh = other.Mesh2MeshList() for element in dic_mesh.keys(): dic_other[element] = dic_mesh[element] for element in dic.keys(): dic_other[element].sort() #print "dsg************************8" #print "dic ",dic #print "*******8" #print "mesh",dic_other #print "dic.__cmp__(dic_o)",dic.__cmp__(dic_other) #print "dsg************************8" return (dic.__cmp__(dic_other)) def addUserPoint(self, pointType, x,y): if pointType == Vertex: point = self.addUserVertex(x,y) if pointType == Hole: point = self._addHole(x,y) if pointType == Region: point = self._addRegion(x,y) return point def addUserVertex(self, x,y): v=Vertex(x, y) self.userVertices.append(v) return v def _addHole(self, x,y): h=Hole(x, y) self.holes.append(h) return h def add_hole(self, x,y, geo_reference=None): """ adds a point, which represents a hole. The point data can have it's own geo_refernece. If geo_reference is None the data is asumed to be absolute """ [[x,y]] = self.geo_reference.change_points_geo_ref([x,y], points_geo_ref=geo_reference) return self._addHole(x, y) def _addRegion(self, x,y): h=Region(x, y) self.regions.append(h) return h def add_region(self, x,y, geo_reference=None, tag=None): """ adds a point, which represents a region. The point data can have it's own geo_refernece. If geo_reference is None the data is asumed to be absolute """ #FIXME: have the user set the tag and resolution here, # but still return the instance, just in case. [[x,y]] = self.geo_reference.change_points_geo_ref([x,y], points_geo_ref=geo_reference) region = self._addRegion(x, y) if tag is not None: region.setTag(tag) return region def build_grid(self, vert_rows, vert_columns): """ Build a grid with vert_rows number of vertex rows and vert_columns number if vertex columns Grid spacing of 1, the origin is the lower left hand corner. FIXME(DSG-DSG) no test. """ for i in range(vert_rows): for j in range(vert_columns): self.addUserVertex(j,i) self.auto_segment() self.generateMesh(mode = "Q", minAngle=20.0) def add_vertices(self, point_data): """ Add user vertices. The point_data can be a list of (x,y) values, a numeric array or a geospatial_data instance. """ point_data = ensure_geospatial(point_data) #print "point_data",point_data # get points relative to the mesh geo_ref points = point_data.get_data_points(geo_reference=self.geo_reference) for point in points: v=Vertex(point[0], point[1]) self.userVertices.append(v) def add_hole_from_polygon(self, polygon, segment_tags=None, geo_reference=None): """ Add a polygon with tags to the current mesh, as a region. The maxArea of the region can be specified. If a geo_reference of the polygon points is given, this is used. If not; The x,y info is assumed to be Easting and Northing, absolute, for the meshes zone. polygon a list of points, in meters that describe the polygon (e.g. [[x1,y1],[x2,y2],...] tags (e.g.{'wall':[0,1,3],'ocean':[2]}) This returns the region instance, so if the user whats to modify it they can. """ return self._add_area_from_polygon(polygon, segment_tags=segment_tags, hole=True, geo_reference=geo_reference ) def add_region_from_polygon(self, polygon, segment_tags=None, max_triangle_area=None, geo_reference=None, region_tag=None): """ Add a polygon with tags to the current mesh, as a region. The maxArea of the region can be specified. If a geo_reference of the polygon points is given, this is used. If not; The x,y info is assumed to be Easting and Northing, absolute, for the meshes zone. polygon a list of points, in meters that describe the polygon (e.g. [[x1,y1],[x2,y2],...] segment_tags (e.g.{'wall':[0,1,3],'ocean':[2]}) add tags to the segements of the region region_tags - add a tag to all of the triangles in the region. This returns the region instance (if a max_triangle_area is given), so if the user whats to modify it they can. """ if max_triangle_area is None and region_tag is None: create_region = False else: create_region = True region = self._add_area_from_polygon(polygon, segment_tags=segment_tags, geo_reference=geo_reference, region=create_region) if max_triangle_area is not None: region.setMaxArea(max_triangle_area) if region_tag is not None: region.setTag(region_tag) return region def _add_area_from_polygon(self, polygon, segment_tags=None, geo_reference=None, hole=False, region=False): """ Add a polygon with tags to the current mesh, as a region. The maxArea of the region can be specified. If a geo_reference of the polygon points is given, this is used. If not; The x,y info is assumed to be Easting and Northing, absolute, for the meshes zone. polygon a list of points, in meters that describe the polygon (e.g. [[x1,y1],[x2,y2],...] segment_tags (e.g.{'wall':[0,1,3],'ocean':[2]})add tags to the segements of the region. This returns the region instance, so if the user whats to modify it they can. """ # Only import this if necessary. # Trying to get pmesh working in an uncompiled environment from anuga.geometry.polygon import point_in_polygon #get absolute values if geo_reference is not None: polygon = geo_reference.get_absolute(polygon) # polygon is now absolute #print "polygon should be absolute",polygon #create points, segs and tags region_dict = {} region_dict['points'] = polygon #Create segments #E.g. [[0,1], [1,2], [2,3], [3,0]] #from polygon #[0,1,2,3] segments = [] N = len(polygon) for i in range(N): lo = i hi = (lo + 1) % N segments.append( [lo, hi] ) region_dict['segments'] = segments region_dict['segment_tags'] = self._tag_dict2list(segment_tags, N) self.addVertsSegs(region_dict) #this is passing absolute values if region is True: #get inner point - absolute values inner_point = point_in_polygon(polygon) inner = self.add_region(inner_point[0], inner_point[1], geo_reference=None) elif hole is True: #get inner point - absolute values inner_point = point_in_polygon(polygon) inner = self.add_hole(inner_point[0], inner_point[1], geo_reference=None) else: inner = None return inner def _tag_dict2list(self, tags, number_of_segs): """ Convert a tag dictionary from this sort of format; #{'wall':[0,3],'ocean':[2]} To a list format # ['wall', '', 'ocean', 'wall'] Note: '' is a default value of nothing """ # FIXME (DSG-DSG): Using '' as a default isn't good. #Try None. # Due to this default this method is too connected to # _add_area_from_polygon segment_tags = ['']*number_of_segs if tags is not None: for key in tags: indices = tags[key] for i in indices: segment_tags[i] = key return segment_tags def add_circle(self, center, radius, segment_count=100, center_geo_reference=None, tag = None, region=False, hole=False): """ center is a point, in absulute or relative to center_geo_ref radius is the radius of the circle segment_count is the number of segments used to represent the circle tag is a string name, given to the segments. If region is True a region object is returned. If hole is True a hole object is returned. (Don't have them both True.) """ # convert center and radius to a polygon cuts = [] factor = 2* math.pi/segment_count for cut in range(segment_count): cuts.append(cut*factor) polygon = [] for cut in cuts: x = center[0] + radius * math.cos(cut) y = center[1] + radius * math.sin(cut) polygon.append([x,y]) # build the tags tags = {tag:range(segment_count)} return self._add_area_from_polygon(polygon, segment_tags=tags, region=region, hole=hole, geo_reference=center_geo_reference) def auto_set_geo_reference(self): """ Automatically set the georeference, based in the minimum x and y user vertex values. Not to be used with the graphical interface Not implemented. Don't implement now. using the new georeferenced points class will change this? """ #to do # find the lower left hand corner [xmin, ymin, xmax, ymax] = self.boxsize() # set all points to that lower left hand corner. # use change_geo_reference new_geo = Geo_reference(self.geo_reference.get_zone(), xmin, ymin) self.change_geo_reference(new_geo) def change_geo_reference(self, new_geo_reference): """ Change from one geo_reference to another. Not to be used with the graphical interface """ # FIXME # change georeference of; #self.userVertices = self.geo_reference.change_points_geo_ref( \ #self.userVertices) #self.holes = self.geo_reference.change_points_geo_ref(self.holes) #self.regions = self.geo_reference.change_points_geo_ref(self.regions) # The above will not work. # since userVertices (etc) is a list of point objects, #not a list of lists. # add a method to the points class to fix this up. def add_segment(self, v1, v2, tag): """ Don't do this function. what will the v's be objects? How is the user suppost to get them? Indexes? If add vertstosegs or add_region_from_polygon is called more than once then the actual index is not obvious. Making this function confusing. """ pass def add_points_and_segments(self, points, segments=None, segment_tags = None): """ Add an outline of the mesh. Points is a list of points a standard representation of points. Segments is a list of tuples of integers. Each tuple defines the start and end of the segment by it's point index. segment_tags is an optional dictionary which is used to add tags to the segments. The key is the tag name, value is the list of segment indexes the tag will apply to. eg. {'wall':[0,3],'ocean':[2]} """ # make sure the points are absolute # Since addVertsSegs will deal with georeferencing. points = ensure_absolute(points) if segments is None: segments = [] for i in range(len(points)-1): segments.append([i, i+1]) #create points, segs and tags region_dict = {} region_dict['points'] = points region_dict['segments'] = segments region_dict['segment_tags'] = self._tag_dict2list(segment_tags, len(segments)) self.addVertsSegs(region_dict) def addVertsSegs(self, outlineDict): """ Add out-line (user Mesh) attributes given a dictionary of the lists points: [(x1,y1),(x2,y2),...] (Tuples of doubles) segments: [(point1,point2),(p3,p4),...] (Tuples of integers) segment_tags: [S1Tag, S2Tag, ...] (list of strings) Assume the values are in Eastings and Northings, with no reference point. eg absolute """ if not outlineDict.has_key('segment_tags'): outlineDict['segment_tags'] = [] for i in range(len(outlineDict['segments'])): outlineDict['segment_tags'].append('') #print "outlineDict['segment_tags']",outlineDict['segment_tags'] #print "outlineDict['points']",outlineDict['points'] #print "outlineDict['segments']",outlineDict['segments'] i_offset = len(self.userVertices) #print "self.userVertices",self.userVertices #print "index_offset",index_offset for point in outlineDict['points']: v=Vertex(point[0]-self.geo_reference.xllcorner, point[1]-self.geo_reference.yllcorner) self.userVertices.append(v) for seg,seg_tag in map(None,outlineDict['segments'], outlineDict['segment_tags']): segObject = Segment(self.userVertices[int(seg[0])+i_offset], self.userVertices[int(seg[1])+i_offset] ) if not seg_tag == '': segObject.set_tag(seg_tag) self.userSegments.append(segObject) def get_triangle_count(self): return len(self.getTriangulation()) def getUserVertices(self): """ Note: The x,y values will be relative to the mesh geo_reference This returns a list of vertex objects """ return self.userVertices def get_user_vertices(self, absolute=True): """ This returns a list of (x,y) values (maybe it should return a geospatical object? It makes it a bit confusing though.) """ pointlist=[] for vertex in self.userVertices: pointlist.append([vertex.x,vertex.y]) spat = Geospatial_data(pointlist, geo_reference=self.geo_reference) return spat.get_data_points(absolute=absolute) def getUserSegments(self): allSegments = self.userSegments + self.alphaUserSegments #print "self.userSegments",self.userSegments #print "self.alphaUserSegments",self.alphaUserSegments #print "allSegments",allSegments return allSegments def deleteUserSegments(self,seg): if self.userSegments.count(seg) == 0: self.alphaUserSegments.remove(seg) pass else: self.userSegments.remove(seg) def clearUserSegments(self): self.userSegments = [] self.alphaUserSegments = [] #FIXME see where this is used. return an array instead def getTriangulation(self): #return self.meshTriangles return self.tri_mesh.triangles.tolist() def getMeshVertices(self): #return self.meshVertices return self.tri_mesh.vertices def getMeshVerticeAttributes(self): #return self.meshVertices return self.tri_mesh.vertex_attributes def getMeshSegments(self): #return self.meshSegments return self.tri_mesh.segments def getMeshSegmentTags(self): #return self.meshSegments return self.tri_mesh.segment_tags def getHoles(self): return self.holes def getRegions(self): return self.regions def isTriangulation(self): if self.meshVertices == []: return False else: return True def addUserSegment(self, v1,v2): """ PRECON: A segment between the two vertices is not already present. Check by calling isUserSegmentNew before calling this function. """ s=Segment( v1,v2) self.userSegments.append(s) return s def generate_mesh(self, maximum_triangle_area="", minimum_triangle_angle=28.0, verbose=True): if verbose is True: silent = '' else: silent = 'Q' self.generateMesh(mode = silent +"pzq"+str(minimum_triangle_angle) +"a"+str(maximum_triangle_area) +"a") #The last a is so areas for regions will be used def generateMesh(self, mode = None, maxArea = None, minAngle= None, isRegionalMaxAreas = True): """ Based on the current user vaules, holes and regions generate a new mesh mode is a string that sets conditions on the mesh generations see triangle_instructions.txt for a definition of the commands PreCondition: maxArea is a double between 1e-20 and 1e30 or is a string. """ #print "mode ",mode if mode == None: self.mode = "" else: self.mode = mode if self.mode.find('p') < 0: self.mode += 'p' #p - read a planar straight line graph. #there must be segments to use this switch # TODO throw an aception if there aren't seg's # it's more comlex than this. eg holes if self.mode.find('z') < 0: self.mode += 'z' # z - Number all items starting from zero # (rather than one) if self.mode.find('n'): self.mode += 'n' # n - output a list of neighboring triangles if self.mode.find('A') < 0: self.mode += 'A' # A - output region attribute list for triangles if not self.mode.find('V') and not self.mode.find('Q'): self.mode += 'V' # V - output info about what Triangle is doing if self.mode.find('q') < 0 and minAngle is not None: # print "**********8minAngle******** ",minAngle min_angle = 'q' + str(minAngle) self.mode += min_angle # z - Number all items starting from zero # (rather than one) if maxArea != None: self.mode += 'a' + str(maxArea) try: self.mode += 'a' + '%20.20f' %maxArea except TypeError: self.mode += 'a' + str(maxArea) #print "self.mode", self.mode #FIXME (DSG-DSG) This isn't explained. if isRegionalMaxAreas: self.mode += 'a' #print "mesh#generateMesh# self.mode",self.mode meshDict = self.Mesh2triangList() #FIXME (DSG-DSG) move below section into generate_mesh.py # & 4 functions eg segment_strings2ints # Actually, because of region_list.append((1.0,2.0,"")) # don't move it, without careful thought #print "*************************!@!@ This is going to triangle !@!@" #print meshDict #print "************************!@!@ This is going to triangle !@!@" #print "meshDict['segmenttaglist']", meshDict['segmenttaglist'] [meshDict['segmenttaglist'], segconverter] = segment_strings2ints(meshDict['segmenttaglist'], initialconversions) #print "regionlist",meshDict['regionlist'] [meshDict['regionlist'], regionconverter] = region_strings2ints(meshDict['regionlist']) #print "%%%%%%%%%%%%%%%%%%%%%%%%%%%regionlist",meshDict['regionlist'] #print "meshDict['segmenttaglist']", meshDict['segmenttaglist' #print "self.mode", self.mode generatedMesh = generate_mesh( meshDict['pointlist'], meshDict['segmentlist'], meshDict['holelist'], meshDict['regionlist'], meshDict['pointattributelist'], meshDict['segmenttaglist'], self.mode, meshDict['pointlist']) #print "%%%%%%%%%%%%%%%%%%%%%%%%%%%generated",generatedMesh generatedMesh['qaa'] = 1 if generatedMesh['generatedsegmentmarkerlist'] is not None: generatedMesh['generatedsegmentmarkerlist'] = \ segment_ints2strings(generatedMesh['generatedsegmentmarkerlist'], segconverter) #print "processed gen",generatedMesh['generatedsegmentmarkerlist'] #print "pmesh mesh generatedMesh['generatedtriangleattributelist']", generatedMesh['generatedtriangleattributelist'] if generatedMesh['generatedtriangleattributelist'] is not None: generatedMesh['generatedtriangleattributelist'] = \ region_ints2strings(generatedMesh['generatedtriangleattributelist'], regionconverter) #print "pmesh mesh generatedMesh['generatedtriangleattributelist']", generatedMesh['generatedtriangleattributelist'] #FIXME (DSG-DSG) move above section into generate_mesh.py if generatedMesh['generatedpointattributelist'] is None or \ generatedMesh['generatedpointattributelist'].shape[1] ==0: self.attributeTitles = [] generatedMesh['generatedpointattributetitlelist']= \ self.attributeTitles #print "################ FROM TRIANGLE" #print "generatedMesh",generatedMesh #print "##################" self.setTriangulation(generatedMesh) def clearTriangulation(self): #Clear the current generated mesh values self.meshTriangles=[] self.meshSegments=[] self.meshVertices=[] def removeDuplicatedUserVertices(self): """Pre-condition: There are no user segments This function will keep the first duplicate """ assert self.getUserSegments() == [] self.userVertices, counter = self.removeDuplicatedVertices( self.userVertices) return counter def removeDuplicatedVertices(self, Vertices): """ This function will keep the first duplicate, remove all others Precondition: Each vertex has a dupindex, which is the list index. Note: this removes vertices that have the same x,y values, not duplicate instances in the Vertices list. """ remove = [] index = 0 for v in Vertices: v.dupindex = index index += 1 t = list(Vertices) t.sort(Point.cmp_xy) length = len(t) behind = 0 ahead = 1 counter = 0 while ahead < length: b = t[behind] ah = t[ahead] if (b.y == ah.y and b.x == ah.x): remove.append(ah.dupindex) behind += 1 ahead += 1 # remove the duplicate vertices remove.sort() remove.reverse() for i in remove: Vertices.pop(i) pass #Remove the attribute that this function added for v in Vertices: del v.dupindex return Vertices,counter # FIXME (DSG-DSG) Move this to geospatial def thinoutVertices(self, delta): """Pre-condition: There are no user segments This function will keep the first duplicate """ assert self.getUserSegments() == [] #t = self.userVertices #self.userVertices =[] boxedVertices = {} thinnedUserVertices =[] delta = round(delta,1) for v in self.userVertices : # tag is the center of the boxes tag = (round(v.x/delta,0)*delta,round(v.y/delta,0)*delta) #this creates a dict of lists of faces, indexed by tag boxedVertices.setdefault(tag,[]).append(v) for [tag,verts] in boxedVertices.items(): min = delta bestVert = None tagVert = Vertex(tag[0],tag[1]) for v in verts: dist = v.DistanceToPoint(tagVert) if (dist yrange: min,max = xmin, xmax else: min,max = ymin, ymax for obj in self.getUserVertices(): obj.x = (obj.x - xmin)/(max- min)*scale + offset obj.y = (obj.y - ymin)/(max- min)*scale + offset if len(obj.attributes) > 0 and attmin0 != attmax0: obj.attributes[0] = (obj.attributes[0]-attmin0)/ \ (attmax0-attmin0)*height_scale if len(obj.attributes) > 1 and attmin1 != attmax1: obj.attributes[1] = (obj.attributes[1]-attmin1)/ \ (attmax1-attmin1)*height_scale for obj in self.getMeshVertices(): obj.x = (obj.x - xmin)/(max- min)*scale + offset obj.y = (obj.y - ymin)/(max- min)*scale + offset if len(obj.attributes) > 0 and attmin0 != attmax0: obj.attributes[0] = (obj.attributes[0]-attmin0)/ \ (attmax0-attmin0)*height_scale if len(obj.attributes) > 1 and attmin1 != attmax1: obj.attributes[1] = (obj.attributes[1]-attmin1)/ \ (attmax1-attmin1)*height_scale for obj in self.getHoles(): obj.x = (obj.x - xmin)/(max- min)*scale + offset obj.y = (obj.y - ymin)/(max- min)*scale + offset for obj in self.getRegions(): obj.x = (obj.x - xmin)/(max- min)*scale + offset obj.y = (obj.y - ymin)/(max- min)*scale + offset [xmin, ymin, xmax, ymax] = self.boxsize() #print [xmin, ymin, xmax, ymax] def boxsizeVerts(self): """ Returns a list of verts denoting a box or triangle that contains verts on the xmin, ymin, xmax and ymax axis. Structure: list of verts """ large = kinds.default_float_kind.MAX xmin= large xmax=-large ymin= large ymax=-large for vertex in self.userVertices: if vertex.x < xmin: xmin = vertex.x xminVert = vertex if vertex.x > xmax: xmax = vertex.x xmaxVert = vertex if vertex.y < ymin: ymin = vertex.y yminVert = vertex if vertex.y > ymax: ymax = vertex.y ymaxVert = vertex verts, count = self.removeDuplicatedVertices([xminVert, xmaxVert, yminVert, ymaxVert]) return verts def boxsize(self): """ Returns a list denoting a box that contains the entire structure of vertices Structure: [xmin, ymin, xmax, ymax] """ large = kinds.default_float_kind.MAX xmin= large xmax=-large ymin= large ymax=-large for vertex in self.userVertices: if vertex.x < xmin: xmin = vertex.x if vertex.x > xmax: xmax = vertex.x if vertex.y < ymin: ymin = vertex.y if vertex.y > ymax: ymax = vertex.y return [xmin, ymin, xmax, ymax] def maxMinVertAtt(self, iatt): """ Returns a list denoting a box that contains the entire structure of vertices Structure: [xmin, ymin, xmax, ymax] """ large = kinds.default_float_kind.MAX min= large max=-large for vertex in self.userVertices: if len(vertex.attributes) > iatt: if vertex.attributes[iatt] < min: min = vertex.attributes[iatt] if vertex.attributes[iatt] > max: max = vertex.attributes[iatt] for vertex in self.meshVertices: if len(vertex.attributes) > iatt: if vertex.attributes[iatt] < min: min = vertex.attributes[iatt] if vertex.attributes[iatt] > max: max = vertex.attributes[iatt] return [min, max] def scaleoffset(self, WIDTH, HEIGHT): """ Returns a list denoting the scale and offset terms that need to be applied when converting mesh co-ordinates onto grid co-ordinates Structure: [scale, xoffset, yoffset] """ OFFSET = 0.05*min([WIDTH, HEIGHT]) [xmin, ymin, xmax, ymax] = self.boxsize() SCALE = min([0.9*WIDTH, 0.9*HEIGHT])/max([xmax-xmin, ymax-ymin]) if SCALE*xmin < OFFSET: xoffset = abs(SCALE*xmin) + OFFSET if SCALE*xmax > WIDTH - OFFSET: xoffset= -(SCALE*xmax - WIDTH + OFFSET) if SCALE*ymin < OFFSET: b = abs(SCALE*ymin)+OFFSET if SCALE*ymax > HEIGHT-OFFSET: b = -(SCALE*ymax - HEIGHT + OFFSET) yoffset = HEIGHT - b return [SCALE, xoffset, yoffset] def exportASCIIobj(self,ofile): """ export a file, ofile, with the format lines: v f (of the triangles) """ fd = open(ofile,'w') self.writeASCIIobj(fd) fd.close() def writeASCIIobj(self,fd): fd.write(" # Triangulation as an obj file\n") numVert = str(len(self.meshVertices)) index1 = 1 for vert in self.meshVertices: vert.index1 = index1 index1 += 1 fd.write("v " + str(vert.x) + " " + str(vert.y) + " " + str(vert.attributes[0]) + "\n") for tri in self.meshTriangles: fd.write("f " + str(tri.vertices[0].index1) + " " + str(tri.vertices[1].index1) + " " + str(tri.vertices[2].index1) + "\n") def exportASCIIsegmentoutlinefile(self,ofile): """Write the boundary user mesh info, eg vertices that are connected to segments, segments """ verts = {} for seg in self.getUserSegments(): verts[seg.vertices[0]] = seg.vertices[0] verts[seg.vertices[1]] = seg.vertices[1] meshDict = self.Mesh2IOOutlineDict(userVertices=verts.values()) export_mesh_file(ofile,meshDict) # exportASCIImeshfile - this function is used def export_mesh_file(self,ofile): """ export a file, ofile, with the format """ dict = self.Mesh2IODict() export_mesh_file(ofile,dict) # FIXME(DSG-DSG):Break this into two functions. #One for the outline points. #One for the mesh points. # Note: this function is not in the gui def exportPointsFile(self,ofile): """ export a points file, ofile. """ mesh_dict = self.Mesh2IODict() #point_dict = {} #point_dict['attributelist'] = {} #this will need to be expanded.. # if attributes are brought back in. #point_dict['geo_reference'] = self.geo_reference if mesh_dict['vertices'] == []: #point_dict['pointlist'] = mesh_dict['points'] geo = Geospatial_data(mesh_dict['points'], geo_reference=self.geo_reference) else: #point_dict['pointlist'] = mesh_dict['vertices'] geo = Geospatial_data(mesh_dict['vertices'], geo_reference=self.geo_reference) geo.export_points_file(ofile, absolute=True) def import_ungenerate_file(self,ofile, tag=None, region_tag=None): """ Imports an ungenerate file, from arcGIS into mesh. This file describes many polygons. ofile is the name of the ungenerated file. Tag is a string name to be taggged on each segment. region_tag is the tag applied to the polygon regions. if it is a string the one value will be assigned to all regions if it is a list the first value in the list will be applied to the first polygon etc. WARNING: size of list and number of polygons isn't checked WARNING values are assumed to be absolute. geo-refs are not taken into account.. """ dict = importUngenerateFile(ofile) default_tag = Segment.get_default_tag() if tag is not None: Segment.set_default_tag(str(tag)) if region_tag is None: self.addVertsSegs(dict) else: if not isinstance(region_tag, list): region_tag = [region_tag]*len(dict['polygons']) for a_tag,polygon in map(None, region_tag, dict['polygons']): segment_tags = {tag:range(len(polygon))} self.add_region_from_polygon(polygon,segment_tags=segment_tags, region_tag=a_tag) Segment.set_default_tag(default_tag) # change the tag back to it's default ########### IO CONVERTERS ################## """ The dict fromat for IO with .tsh files is; (the triangulation) vertices: [[x1,y1],[x2,y2],...] (lists of doubles) vertex_attributes: [[a11,a12,...],[a21,a22],...] (lists of doubles) vertex_attribute_titles:[A1Title, A2Title ...] (A list of strings) segments: [[v1,v2],[v3,v4],...] (lists of integers) segment_tags : [tag,tag,...] list of strings triangles : [(v1,v2,v3), (v4,v5,v6),....] lists of points triangle tags: [s1,s2,...] A list of strings triangle neighbors: [[t1,t2,t3], [t4,t5,t6],..] lists of triangles (the outline) points: [[x1,y1],[x2,y2],...] (lists of doubles) point_attributes: [[a11,a12,...],[a21,a22],...] (lists of doubles) outline_segments: [[point1,point2],[p3,p4],...] (lists of integers) outline_segment_tags : [tag1,tag2,...] list of strings holes : [[x1,y1],...](List of doubles, one inside each hole region) regions : [ [x1,y1],...] (List of 4 doubles) region_tags : [tag1,tag2,...] (list of strings) region_max_areas: [ma1,ma2,...] (A list of doubles) {Convension: A -ve max area means no max area} """ def Mesh2IODict(self): """ Convert the triangulation and outline info of a mesh to a dictionary structure """ dict = self.Mesh2IOTriangulationDict() dict_mesh = self.Mesh2IOOutlineDict() for element in dict_mesh.keys(): dict[element] = dict_mesh[element] # add the geo reference dict['geo_reference'] = self.geo_reference return dict def Mesh2IOTriangulationDict(self): """ Convert the Mesh to a dictionary of lists describing the triangulation variables; Used to produce .tsh file """ meshDict = {} if self.tri_mesh is not None: meshDict['triangles'] = self.tri_mesh.triangles meshDict['triangle_tags'] = self.tri_mesh.triangle_tags #print "mesh meshDict['triangle_tags']", meshDict['triangle_tags'] meshDict['triangle_neighbors'] = self.tri_mesh.triangle_neighbors meshDict['vertices'] = self.tri_mesh.vertices meshDict['vertex_attributes'] = self.tri_mesh.vertex_attributes meshDict['vertex_attribute_titles'] = \ self.tri_mesh.vertex_attribute_titles meshDict['segments'] = self.tri_mesh.segments meshDict['segment_tags'] = self.tri_mesh.segment_tags else: meshDict['triangles'] = [] meshDict['triangle_tags'] = [] meshDict['triangle_neighbors'] = [] meshDict['vertices'] = [] meshDict['vertex_attributes'] = [] meshDict['vertex_attribute_titles'] = [] meshDict['segments'] = [] meshDict['segment_tags'] = [] #print "mesh.Mesh2IOTriangulationDict*)*)" #print meshDict #print "mesh.Mesh2IOTriangulationDict*)*)" return meshDict def Mesh2IOOutlineDict(self, userVertices=None, userSegments=None, holes=None, regions=None): """ Convert the mesh outline to a dictionary of the lists needed for the triang module; Note, this adds an index attribute to the user Vertex objects. Used to produce .tsh file and output to triangle """ if userVertices is None: userVertices = self.getUserVertices() if userSegments is None: userSegments = self.getUserSegments() if holes is None: holes = self.getHoles() if regions is None: regions = self.getRegions() meshDict = {} #print "userVertices",userVertices #print "userSegments",userSegments pointlist=[] pointattributelist=[] index = 0 for vertex in userVertices: vertex.index = index pointlist.append([vertex.x,vertex.y]) pointattributelist.append(vertex.attributes) index += 1 meshDict['points'] = pointlist meshDict['point_attributes'] = pointattributelist segmentlist=[] segmenttaglist=[] for seg in userSegments: segmentlist.append([seg.vertices[0].index,seg.vertices[1].index]) segmenttaglist.append(seg.tag) meshDict['outline_segments'] =segmentlist meshDict['outline_segment_tags'] =segmenttaglist holelist=[] for hole in holes: holelist.append([hole.x,hole.y]) meshDict['holes'] = holelist regionlist=[] regiontaglist = [] regionmaxarealist = [] for region in regions: regionlist.append([region.x,region.y]) regiontaglist.append(region.getTag()) if (region.getMaxArea() != None): regionmaxarealist.append(region.getMaxArea()) else: regionmaxarealist.append(NOMAXAREA) meshDict['regions'] = regionlist meshDict['region_tags'] = regiontaglist meshDict['region_max_areas'] = regionmaxarealist #print "*(*(" #print meshDict #print meshDict['regionlist'] #print "*(*(" return meshDict def IOTriangulation2Mesh(self, genDict): """ Set the mesh attributes given an tsh IO dictionary """ #Clear the current generated mesh values self.tri_mesh = None self.tri_mesh = Rigid_triangulation( genDict['triangles'] ,genDict['segments'] ,genDict['vertices'] ,genDict['triangle_tags'] ,genDict['triangle_neighbors'] ,genDict['segment_tags'] ,genDict['vertex_attributes'] ,genDict['vertex_attribute_titles'] ) self.attributeTitles = genDict['vertex_attribute_titles'] self.maxVertexIndex = len(genDict['vertices']) #print "self.maxVertexIndex ", self.maxVertexIndex def IOOutline2Mesh(self, genDict): """ Set the outline (user Mesh attributes) given a IO tsh dictionary mesh is an instance of a mesh object """ #Clear the current user mesh values self.clearUserSegments() self.userVertices=[] self.Holes=[] self.Regions=[] #print "mesh.IOOutline2Mesh@#@#@#" #print "genDict",genDict #print "@#@#@#" #index = 0 for point in genDict['points']: v=Vertex(point[0], point[1]) #v.index = index #index +=1 self.userVertices.append(v) #index = 0 for seg,tag in map(None,genDict['outline_segments'], genDict['outline_segment_tags']): segObject = Segment( self.userVertices[int(seg[0])], self.userVertices[int(seg[1])], tag = tag ) #segObject.index = index #index +=1 self.userSegments.append(segObject) # Remove the loading of attribute info. # Have attribute info added using least_squares in pyvolution # index = 0 # for att in genDict['point_attributes']: # if att == None: # self.userVertices[index].setAttributes([]) # else: # self.userVertices[index].setAttributes(att) # index += 1 #index = 0 for point in genDict['holes']: h=Hole(point[0], point[1]) #h.index = index #index +=1 self.holes.append(h) #index = 0 for reg,att,maxArea in map(None, genDict['regions'], genDict['region_tags'], genDict['region_max_areas']): if maxArea > 0: # maybe I should ref NOMAXAREA? Prob' not though Object = Region( reg[0], reg[1], tag = att, maxArea = maxArea) else: Object = Region( reg[0], reg[1], tag = att) #Object.index = index #index +=1 self.regions.append(Object) def importUngenerateFile(ofile): """ import a file, ofile, with the format [poly] poly format: First line: <# of vertices> Following lines: last line: "END" Note: These are clockwise. """ fd = open(ofile,'r') Dict = readUngenerateFile(fd) fd.close() return Dict def readUngenerateFile(fd): """ import a file, ofile, with the format [poly] poly format: First line: <# of polynomial> Following lines: last line: "END" """ END_DELIMITER = 'END' points = [] segments = [] polygons = [] isEnd = False line = fd.readline() #not used <# of polynomial> while not isEnd: poly = [] line = fd.readline() fragments = line.split() x = float(fragments.pop(0)) y = float(fragments.pop(0)) points.append([x,y]) poly.append([x,y]) PreviousVertIndex = len(points)-1 firstVertIndex = PreviousVertIndex line = fd.readline() #Read the next line while not line.startswith(END_DELIMITER): #print "line >" + line + "<" fragments = line.split() x = float(fragments.pop(0)) y = float(fragments.pop(0)) points.append([x,y]) poly.append([x,y]) thisVertIndex = len(points)-1 segment = [PreviousVertIndex,thisVertIndex] segments.append(segment) PreviousVertIndex = thisVertIndex line = fd.readline() #Read the next line # If the last and first segments are the same, # Remove the last segment and the last vertex # then add a segment from the second last vert to the 1st vert thisVertIndex = len(points)-1 firstVert = points[firstVertIndex] thisVert = points[thisVertIndex] #print "firstVert",firstVert #print "thisVert",thisVert if (firstVert[0] == thisVert[0] and firstVert[1] == thisVert[1]): points.pop() segments.pop() poly.pop() thisVertIndex = len(points)-1 segments.append([thisVertIndex, firstVertIndex]) line = fd.readline() # read <# of polynomial> OR END #print "line >>" + line + "<<" # do poly stuff here polygons.append(poly) if line.startswith(END_DELIMITER): isEnd = True #print "points", points #print "segments", segments ungenerated_dict = {} ungenerated_dict['points'] = points ungenerated_dict['segments'] = segments ungenerated_dict['polygons'] = polygons return ungenerated_dict def importMeshFromFile(ofile): """returns a mesh object, made from a points file or .tsh/.msh file Often raises IOError,RuntimeError """ newmesh = None if (ofile[-4:]== ".pts" or ofile[-4:]== ".txt" or \ ofile[-4:]== ".csv"): geospatial = Geospatial_data(ofile) dict = {} dict['points'] = geospatial.get_data_points(absolute=False) dict['outline_segments'] = [] dict['outline_segment_tags'] = [] dict['regions'] = [] dict['region_tags'] = [] dict['region_max_areas'] = [] dict['holes'] = [] newmesh= Mesh(geo_reference = geospatial.geo_reference) newmesh.IOOutline2Mesh(dict) counter = newmesh.removeDuplicatedUserVertices() if (counter >0): log.critical("%i duplicate vertices removed from dataset" % counter) elif (ofile[-4:]== ".tsh" or ofile[-4:]== ".msh"): dict = import_mesh_file(ofile) #print "********" #print "zq mesh.dict",dict #print "********" newmesh= Mesh() newmesh.IOOutline2Mesh(dict) newmesh.IOTriangulation2Mesh(dict) else: raise RuntimeError if dict.has_key('geo_reference') and not dict['geo_reference'] == None: newmesh.geo_reference = dict['geo_reference'] return newmesh def loadPickle(currentName): fd = open(currentName) mesh = pickle.load(fd) fd.close() return mesh def square_outline(side_length = 1,up = "top", left = "left", right = "right", down = "bottom", regions = False): a = Vertex (0,0) b = Vertex (0,side_length) c = Vertex (side_length,0) d = Vertex (side_length,side_length) s2 = Segment(b,d, tag = up) s3 = Segment(b,a, tag = left) s4 = Segment(d,c, tag = right) s5 = Segment(a,c, tag = down) if regions: e = Vertex (side_length/2,side_length/2) s6 = Segment(a,e, tag = down + left) s7 = Segment(b,e, tag = up + left) s8 = Segment(c,e, tag = down + right) s9 = Segment(d,e, tag = up + right) r1 = Region(side_length/2,3.*side_length/4, tag = up) r2 = Region(1.*side_length/4,side_length/2, tag = left) r3 = Region(3.*side_length/4,side_length/2, tag = right) r4 = Region(side_length/2,1.*side_length/4, tag = down) mesh = Mesh(userVertices=[a,b,c,d,e], userSegments=[s2,s3,s4,s5,s6,s7,s8,s9], regions = [r1,r2,r3,r4]) else: mesh = Mesh(userVertices=[a,b,c,d], userSegments=[s2,s3,s4,s5]) return mesh def region_strings2ints(region_list): """Given a list of (x_int,y_int,tag_string) lists it returns a list of (x_int,y_int,tag_int) and a list to convert the tag_int's back to the tag_strings """ # Make sure "" has an index of 0 region_list.reverse() region_list.append((1.0,2.0,"")) region_list.reverse() convertint2string = [] for i in xrange(len(region_list)): convertint2string.append(region_list[i][2]) if len(region_list[i]) == 4: # there's an area value region_list[i] = (region_list[i][0], region_list[i][1],i,region_list[i][3]) elif len(region_list[i]) == 3: # no area value region_list[i] = (region_list[i][0],region_list[i][1],i) else: log.critical("The region list has a bad size") # raise an error .. raise Error #remove "" from the region_list region_list.pop(0) return [region_list, convertint2string] def region_ints2strings(region_list,convertint2string): """Reverses the transformation of region_strings2ints """ #print 'region_ints2strings region_list', region_list returned_region_list = [] # may not need (not region_list[0] == []) # or region_list[0] == [0.0] if (not region_list[0] == []): # or region_list[0] == [0.0]: #print "in loop" for i in xrange(len(region_list)): temp = region_list[i] returned_region_list.append(convertint2string[int(temp[0])]) return returned_region_list def segment_ints2strings(intlist, convertint2string): """Reverses the transformation of segment_strings2ints """ stringlist = [] for x in intlist: stringlist.append(convertint2string[int(x)]) return stringlist def segment_strings2ints(stringlist, preset): """Given a list of strings return a list of 0 to n ints which represent the strings and a converting list of the strings, indexed by 0 to n ints. Also, input an initial converting list of the strings Note, the converting list starts off with ["internal boundary", "external boundary", "internal boundary"] example input and output input ["a","b","a","c"],["c"] output [[2, 1, 2, 0], ['c', 'b', 'a']] the first element in the converting list will be overwritten with "". ?This will become the third item in the converting list? # note, order the initial converting list is important, since the index = the int tag """ nodups = unique(stringlist) # note, order is important, the index = the int tag #preset = ["internal boundary", "external boundary"] #Remove the preset tags from the list with no duplicates nodups = [x for x in nodups if x not in preset] try: nodups.remove("") # this has to go to zero except ValueError: pass # Add the preset list at the beginning of no duplicates preset.reverse() nodups.extend(preset) nodups.reverse() convertstring2int = {} convertint2string = [] index = 0 for x in nodups: convertstring2int[x] = index convertint2string.append(x) index += 1 convertstring2int[""] = 0 intlist = [] for x in stringlist: intlist.append(convertstring2int[x]) return [intlist, convertint2string] def unique(s): """Return a list of the elements in s, but without duplicates. For example, unique([1,2,3,1,2,3]) is some permutation of [1,2,3], unique("abcabc") some permutation of ["a", "b", "c"], and unique(([1, 2], [2, 3], [1, 2])) some permutation of [[2, 3], [1, 2]]. For best speed, all sequence elements should be hashable. Then unique() will usually work in linear time. If not possible, the sequence elements should enjoy a total ordering, and if list(s).sort() doesn't raise TypeError it's assumed that they do enjoy a total ordering. Then unique() will usually work in O(N*log2(N)) time. If that's not possible either, the sequence elements must support equality-testing. Then unique() will usually work in quadratic time. """ n = len(s) if n == 0: return [] # Try using a dict first, as that's the fastest and will usually # work. If it doesn't work, it will usually fail quickly, so it # usually doesn't cost much to *try* it. It requires that all the # sequence elements be hashable, and support equality comparison. u = {} try: for x in s: u[x] = 1 except TypeError: del u # move on to the next method else: return u.keys() # We can't hash all the elements. Second fastest is to sort, # which brings the equal elements together; then duplicates are # easy to weed out in a single pass. # NOTE: Python's list.sort() was designed to be efficient in the # presence of many duplicate elements. This isn't true of all # sort functions in all languages or libraries, so this approach # is more effective in Python than it may be elsewhere. try: t = list(s) t.sort() except TypeError: del t # move on to the next method else: assert n > 0 last = t[0] lasti = i = 1 while i < n: if t[i] != last: t[lasti] = last = t[i] lasti += 1 i += 1 return t[:lasti] # Brute force is all that's left. u = [] for x in s: if x not in u: u.append(x) return u # FIXME (DSG-DSG) # To do- # Create a clear interface. eg # have the interface methods more at the top of this file and add comments # for the interface functions/methods, use function_name (not functionName), #Currently #function_name methods assume absolute values. Geo-refs can be passed in. # # instead of functionName if __name__ == "__main__": pass