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

Last change on this file since 7753 was 7711, checked in by hudson, 14 years ago

Refactored geometry classes to live in their own folder.

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