source: branches/source_numpy_conversion/anuga/pmesh/mesh.py @ 6982

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

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