source: anuga_core/source/anuga/pmesh/mesh.py @ 6982

Last change on this file since 6982 was 6156, checked in by rwilson, 16 years ago

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