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

Last change on this file since 4955 was 4955, checked in by steve, 16 years ago

Fixed problems with entries of integer arrays not being integers, but
arrays. So needed to coerce to int in a couple of places.

Also mkstemp was being imported from the wrong module in test_system_tools

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