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

Last change on this file since 5421 was 5209, checked in by duncan, 16 years ago

comments

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