source: inundation/pmesh/mesh.py @ 3158

Last change on this file since 3158 was 3100, checked in by duncan, 19 years ago

added import_ungenerate_file method
File size: 125.1 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 percistant 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"""
15import load_mesh.loadASCII
16import alpha_shape.alpha_shape
17
18import sys
19import math
20import mesh_engine.triang as triang
21import re
22import os
23import pickle
24
25import types
26import exceptions
27
28import load_mesh
29from coordinate_transforms.geo_reference import Geo_reference,DEFAULT_ZONE
30from utilities.polygon import point_in_polygon
31
32SET_COLOUR='red'
33
34#FIXME: this is not tested.
35from utilities.numerical_tools import gradient
36
37
38
39# 1st and third values must be the same
40# FIXME: maybe make this a switch that the user can change? - DSG
41initialconversions = ['', 'exterior', '']
42
43#from os import sep
44#sys.path.append('..'+sep+'pmesh')
45#print "sys.path",sys.path
46
47class MeshObject:
48    """
49    An abstract superclass for the basic mesh objects, eg vertex, segment,
50    triangle.
51    """
52    def __init__(self):
53        pass
54   
55class Point(MeshObject): 
56    """
57    Define a point in a 2D space.
58    """
59    def __init__(self,X,Y):
60        __slots__ = ['x','y']
61        self.x=X
62        self.y=Y
63       
64    def DistanceToPoint(self, OtherPoint):
65        """
66        Returns the distance from this point to another
67        """
68        SumOfSquares = ((self.x - OtherPoint.x)**2) + ((self.y - OtherPoint.y)**2)
69        return math.sqrt(SumOfSquares)
70       
71    def IsInsideCircle(self, Center, Radius):
72        """
73        Return 1 if this point is inside the circle,
74        0 otherwise
75        """
76       
77        if (self.DistanceToPoint(Center)<Radius):
78            return 1
79        else:
80            return 0
81       
82    def __repr__(self):
83        return "(%f,%f)" % (self.x,self.y) 
84
85    def cmp_xy(self, point):
86        if self.x < point.x:
87            return -1
88        elif self.x > point.x:
89            return 1
90        else:           
91            if self.y < point.y:
92                return -1
93            elif self.y > point.y:
94                return 1
95            else:
96                return 0
97       
98    def same_x_y(self, point):
99        if self.x == point.x and self.y == point.y:
100            return True
101        else:
102            return False
103       
104           
105
106class Vertex(Point):
107    """
108    A point on the mesh.
109    Object attributes based on the Triangle program
110    """
111    def __init__(self,X,Y, attributes = None):
112        __slots__ = ['x','y','attributes']
113       
114        assert (type(X) == types.FloatType or type(X) == types.IntType)
115        assert (type(Y) == types.FloatType or type(Y) == types.IntType)
116        self.x=X
117        self.y=Y       
118        self.attributes=[] 
119       
120        if attributes is None:
121            self.attributes=[]
122        else:
123            self.attributes=attributes
124   
125
126    def setAttributes(self,attributes):
127        """
128        attributes is a list.
129        """
130        self.attributes = attributes
131       
132    VERTEXSQUARESIDELENGTH = 6
133    def draw(self, canvas, tags, colour = 'black',scale = 1, xoffset = 0, yoffset =0, ):
134        x =  scale*(self.x + xoffset)
135        y = -1*scale*(self.y + yoffset)  # - since for a canvas - is up
136        #print "draw x:", x
137        #print "draw y:", y
138        cornerOffset= self.VERTEXSQUARESIDELENGTH/2
139
140        # A hack to see the vert tags
141        # note: there will be many tags, since tags will not be removed
142        #when zooming
143        #canvas.create_text(x+ 2*cornerOffset,
144        #                   y+ 2*cornerOffset,
145        #                        text=tags)
146       
147        return canvas.create_rectangle(x-cornerOffset,
148                                       y-cornerOffset,
149                                       x+cornerOffset,
150                                       y+cornerOffset,
151                                       tags = tags,
152                                       outline=colour,
153                                       fill = 'white')
154   
155        #return tags
156     
157    def __repr__(self):
158        return "[(%f,%f),%r]" % (self.x,self.y,self.attributes)
159   
160class Hole(Point):
161    """
162    A region of the mesh were no triangles are generated.
163    Defined by a point in the hole enclosed by segments.
164    """
165
166    HOLECORNERLENGTH = 6
167   
168    def draw(self, canvas, tags, colour = 'purple',scale = 1, xoffset = 0, yoffset =0, ):
169        x =  scale*(self.x + xoffset)
170        y = -1*scale*(self.y + yoffset)  # - since for a canvas - is up
171        #print "draw x:", x
172        #print "draw y:", y
173        cornerOffset= self.HOLECORNERLENGTH/2
174        return canvas.create_oval(x-cornerOffset,
175                                       y-cornerOffset,
176                                       x+cornerOffset,
177                                       y+cornerOffset,
178                                       tags = tags,
179                                       outline=colour,
180                                       fill = 'white')
181   
182class Region(Point):
183    """
184    A region of the mesh, defined by a point in the region
185    enclosed by segments. Used to tag areas.
186    """
187    CROSSLENGTH = 6
188    TAG = 0
189    MAXAREA = 1
190   
191    def __init__(self,X,Y, tag = None, maxArea = None):
192        """Precondition: tag is a string and maxArea is a real
193        """
194        # This didn't work. 
195        #super(Region,self)._init_(self,X,Y)
196        self.x=X
197        self.y=Y   
198        self.attributes =[] # index 0 is the tag string
199                            #optoinal index 1 is the max triangle area
200                            #NOTE the size of this attribute is assumed
201                            # to be 1 or 2 in regionstrings2int
202        if tag is None:
203            self.attributes.append("")
204        else:
205            self.attributes.append(tag) #this is a string
206           
207        if maxArea is not None:
208            self.setMaxArea(maxArea) # maxArea is a number
209           
210    def getTag(self,):
211        return self.attributes[self.TAG]
212   
213    def setTag(self,tag):
214        self.attributes[self.TAG] = tag
215       
216    def getMaxArea(self):
217        """ Returns the Max Area of a Triangle or
218        None, if the Max Area has not been set.
219        """
220        if self.isMaxArea():
221            return self.attributes[1]
222        else:
223            return None
224   
225    def setMaxArea(self,MaxArea):
226        if self.isMaxArea(): 
227            self.attributes[self.MAXAREA] = float(MaxArea)
228        else:
229            self.attributes.append( float(MaxArea) )
230   
231    def deleteMaxArea(self):
232        if self.isMaxArea():
233            self.attributes.pop(self.MAXAREA)
234           
235    def isMaxArea(self):
236        return len(self.attributes)> 1
237   
238    def draw(self, canvas, tags, scale=1, xoffset = 0, yoffset =0, colour = "black"):
239        """
240        Draw a black cross, returning the objectID
241        """
242        x =  scale*(self.x + xoffset)
243        y = -1*scale*(self.y + yoffset) 
244        cornerOffset= self.CROSSLENGTH/2
245        return canvas.create_polygon(x,
246                                     y-cornerOffset,
247                                     x,
248                                     y,
249                                     x+cornerOffset,
250                                     y,
251                                     x,
252                                     y,
253                                     x,
254                                     y+cornerOffset,
255                                     x,
256                                     y,
257                                     x-cornerOffset,
258                                     y,
259                                     x,
260                                     y,
261                                     tags = tags,
262                                     outline = colour,fill = '')
263   
264    def __repr__(self):
265        if self.isMaxArea():
266            area = self.getMaxArea() 
267            return "(%f,%f,%s,%f)" % (self.x,self.y,
268                                      self.getTag(), self.getMaxArea())
269        else:
270            return "(%f,%f,%s)" % (self.x,self.y,
271                                   self.getTag())
272       
273class Triangle(MeshObject):
274    """
275    A triangle element, defined by 3 vertices.
276    Attributes based on the Triangle program.
277    """
278
279    def __init__(self, vertex1, vertex2, vertex3, attribute = None, neighbors = None ):
280        """
281        Vertices, the initial arguments, are listed in counterclockwise order.
282        """
283        self.vertices= [vertex1,vertex2, vertex3 ]
284       
285        if attribute is None:
286            self.attribute =""
287        else:
288            self.attribute = attribute #this is a string
289           
290        if neighbors is None:
291            self.neighbors=[]
292        else:
293            self.neighbors=neighbors
294
295    def replace(self,new_triangle):
296        self = new_triangle
297
298    def longestSideID(self):
299        ax = self.vertices[0].x
300        ay = self.vertices[0].y
301       
302        bx = self.vertices[1].x
303        by = self.vertices[1].y
304       
305        cx = self.vertices[2].x
306        cy = self.vertices[2].y
307
308        lenA = ((cx-bx)**2+(cy-by)**2)**0.5
309        lenB = ((ax-cx)**2+(ay-cy)**2)**0.5
310        lenC = ((bx-ax)**2+(by-ay)**2)**0.5
311 
312        len = [lenA,lenB,lenC]
313        return len.index(max(len))
314
315    def rotate(self,offset):
316        """
317        permute the order of the sides of the triangle
318        offset must be 0,1 or 2
319        """
320
321        if offset == 0:
322            pass
323        else:
324            if offset == 1:
325                self.vertices = [self.vertices[1],self.vertices[2],self.vertices[0]]
326                self.neighbors = [self.neighbors[1],self.neighbors[2],self.neighbors[0]]
327            if offset == 2:
328                self.vertices = [self.vertices[2],self.vertices[0],self.vertices[1]]
329                self.neighbors = [self.neighbors[2],self.neighbors[0],self.neighbors[1]]
330
331    def rotate_longest_side(self):
332        self.rotate(self.longestSideID())
333
334    def getVertices(self):
335        return self.vertices
336   
337    def get_vertices(self):
338        """
339        Return a list of the vertices.  The x and y values will be relative
340        Easting and Northings for the zone of the current geo_ref.
341        """
342        return self.vertices
343   
344    def calcArea(self):
345        ax = self.vertices[0].x
346        ay = self.vertices[0].y
347       
348        bx = self.vertices[1].x
349        by = self.vertices[1].y
350       
351        cx = self.vertices[2].x
352        cy = self.vertices[2].y
353       
354        return abs((bx*ay-ax*by)+(cx*by-bx*cy)+(ax*cy-cx*ay))/2
355   
356    def calcP(self):
357        #calculate the perimeter
358        ax = self.vertices[0].x
359        ay = self.vertices[0].y
360       
361        bx = self.vertices[1].x
362        by = self.vertices[1].y
363       
364        cx = self.vertices[2].x
365        cy = self.vertices[2].y
366
367        a =  ((cx-bx)**2+(cy-by)**2)**0.5 
368        b =  ((ax-cx)**2+(ay-cy)**2)**0.5 
369        c =  ((bx-ax)**2+(by-ay)**2)**0.5 
370
371        return a+b+c
372           
373    def setNeighbors(self,neighbor1 = None, neighbor2 = None, neighbor3 = None):
374        """
375        neighbor1 is the triangle opposite vertex1 and so on.
376        Null represents no neighbor
377        """
378        self.neighbors = [neighbor1, neighbor2, neighbor3]
379
380    def setAttribute(self,attribute):
381        """
382        neighbor1 is the triangle opposite vertex1 and so on.
383        Null represents no neighbor
384        """
385        self.attribute = attribute #this is a string
386       
387    def __repr__(self):
388        return "[%s,%s]" % (self.vertices,self.attribute)
389       
390
391    def draw(self, canvas, tags, scale=1, xoffset = 0, yoffset =0, colour = "green"):
392        """
393        Draw a triangle, returning the objectID
394        """
395        return canvas.create_polygon(scale*(self.vertices[1].x + xoffset),
396                                     scale*-1*(self.vertices[1].y + yoffset),
397                                     scale*(self.vertices[0].x + xoffset),
398                                     scale*-1*(self.vertices[0].y + yoffset),
399                                     scale*(self.vertices[2].x + xoffset),
400                                     scale*-1*(self.vertices[2].y + yoffset),
401                                     tags = tags,
402                                     outline = colour,fill = '')
403
404class Segment(MeshObject):
405    """
406    Segments are edges whose presence in the triangulation is enforced.
407   
408    """
409    def __init__(self, vertex1, vertex2, tag = None ):
410        """
411        Each segment is specified by listing the vertices of its endpoints
412        The vertices are Vertex objects
413        """
414        assert(vertex1 != vertex2)
415        self.vertices = [vertex1,vertex2 ]
416       
417        if tag is None:
418            self.tag = self.__class__.default
419        else:
420            self.tag = tag #this is a string
421       
422    def __repr__(self):
423        return "[%s,%s]" % (self.vertices,self.tag)
424           
425       
426    def draw(self, canvas, tags,scale=1 , xoffset=0 , yoffset=0,colour='blue' ):
427        x=[]
428        y=[]
429        for end in self.vertices:
430            #end.draw(canvas,scale, xoffset, yoffset ) # draw the vertices
431            x.append(scale*(end.x + xoffset))
432            y.append(-1*scale*(end.y + yoffset)) # - since for a canvas - is up
433       
434        return canvas.create_line(x[0], y[0], x[1], y[1],
435                                  tags = tags,fill=colour)
436    def set_tag(self,tag):
437        self.tag = tag
438       
439    # Class methods
440    def set_default_tag(cls, default):
441        cls.default = default
442   
443    def get_default_tag(cls):
444        return cls.default
445   
446    set_default_tag = classmethod(set_default_tag) 
447    get_default_tag = classmethod(get_default_tag)
448
449Segment.set_default_tag("")       
450
451class Mesh:
452    """
453    Representation of a 2D triangular mesh.
454    User attributes describe the mesh region/segments/vertices/attributes
455
456    mesh attributes describe the mesh that is produced eg triangles and vertices.
457    All point information is relative to the geo_reference passed in
458   
459   
460    """
461
462    def __repr__(self):
463        return """
464        mesh Triangles: %s 
465        mesh Attribute Titles: %s 
466        mesh Segments: %s 
467        mesh Vertices: %s 
468        user Segments: %s 
469        user Vertices: %s 
470        holes: %s 
471        regions: %s""" % (self.meshTriangles,
472                                self.attributeTitles,
473                                self.meshSegments,
474                                self.meshVertices,
475                                self.getUserSegments(),
476                                self.userVertices,
477                                self.holes,
478                                self.regions) 
479   
480    def __init__(self,
481                 userSegments=None,
482                 userVertices=None,
483                 holes=None,
484                 regions=None,
485                 geo_reference=None):
486        self.meshTriangles=[] 
487        self.attributeTitles=[] 
488        self.meshSegments=[]
489        self.meshVertices=[]
490
491        self.setID={}
492        #a dictionary of names.
493        #multiple sets are allowed, but the gui does not yet
494        #support this
495       
496        self.setID['None']=0
497        #contains the names of the sets pointing to the indexes
498        #in the list.
499       
500        self.sets=[[]]
501        #Contains the lists of triangles (triangle sets)
502
503       
504        self.visualise_graph = True
505
506        if userSegments is None:
507            self.userSegments=[]
508        else:
509            self.userSegments=userSegments
510        self.alphaUserSegments=[]
511           
512        if userVertices is None:
513            self.userVertices=[]
514        else:
515            self.userVertices=userVertices
516           
517        if holes is None:
518            self.holes=[]
519        else:
520            self.holes=holes
521           
522        if regions is None:
523            self.regions=[]
524        else:
525            self.regions=regions
526
527        if geo_reference is None:
528            self.geo_reference = Geo_reference(DEFAULT_ZONE,0,0) 
529        else:
530            self.geo_reference = geo_reference
531           
532    def __cmp__(self,other):
533       
534        # A dic for the initial m
535        dic = self.Mesh2triangList()
536        dic_mesh = self.Mesh2MeshList()
537        for element in dic_mesh.keys():
538            dic[element] = dic_mesh[element]
539        for element in dic.keys():
540            dic[element].sort()
541           
542        # A dic for the exported/imported m
543        dic_other = other.Mesh2triangList()
544        dic_mesh = other.Mesh2MeshList()
545        for element in dic_mesh.keys():
546            dic_other[element] = dic_mesh[element]
547        for element in dic.keys():
548            dic_other[element].sort()
549
550        #print "dsg************************8"
551        #print "dic ",dic
552        #print "*******8"
553        #print "mesh",dic_other
554        #print "dic.__cmp__(dic_o)",dic.__cmp__(dic_other)
555        #print "dsg************************8"
556       
557        return (dic.__cmp__(dic_other))
558   
559    def addUserPoint(self, pointType, x,y):
560        if pointType == Vertex:
561            point = self.addUserVertex(x,y)
562        if pointType == Hole:
563            point = self._addHole(x,y)
564        if pointType == Region:
565            point = self._addRegion(x,y)
566        return point
567   
568    def addUserVertex(self, x,y):
569        v=Vertex(x, y)
570        self.userVertices.append(v)
571        return v
572
573    def _addHole(self, x,y):
574        h=Hole(x, y)
575        self.holes.append(h)
576        return h
577   
578    def add_hole(self, x,y, geo_reference=None):
579        """
580        adds a point, which represents a hole.
581
582        The point data can have it's own geo_refernece.
583        If geo_reference is None the data is asumed to be absolute
584        """
585        [[x,y]] = self.geo_reference.change_points_geo_ref([x,y],
586                                                 points_geo_ref=geo_reference)
587        return self._addHole(x, y)
588
589    def _addRegion(self, x,y):
590        h=Region(x, y)
591        self.regions.append(h)
592        return h
593   
594    def add_region(self, x,y, geo_reference=None):
595        """
596        adds a point, which represents a region.
597
598        The point data can have it's own geo_refernece.
599        If geo_reference is None the data is asumed to be absolute
600        """
601        [[x,y]] = self.geo_reference.change_points_geo_ref([x,y],
602                                                 points_geo_ref=geo_reference)
603        return self._addRegion(x, y)
604
605    # Depreciated
606    def addRegionEN(self, x,y):
607        print "depreciated, use add_region"
608        return self.add_region(x,y)
609
610   
611   
612    def add_hole_from_polygon(self, polygon, tags=None, geo_reference=None):
613        """
614        Add a polygon with tags to the current mesh, as a region.
615        The maxArea of the region can be specified.
616
617        If a geo_reference of the polygon points is given, this is used.
618        If not;
619        The x,y info is assumed to be Easting and Northing, absolute,
620        for the meshes zone.
621
622        polygon a list of points, in meters that describe the polygon
623             (e.g. [[x1,y1],[x2,y2],...]
624        tags (e.g.{'wall':[0,1,3],'ocean':[2]})
625
626        This returns the region instance, so if the user whats to modify
627        it they can.       
628        """
629        return self._add_area_from_polygon(polygon, tags=tags,
630                                           hole=True,
631                                           geo_reference=geo_reference
632                                           )
633
634       
635    def add_region_from_polygon(self, polygon, tags=None,
636                                max_triangle_area=None, geo_reference=None):
637        """
638        Add a polygon with tags to the current mesh, as a region.
639        The maxArea of the region can be specified.
640
641        If a geo_reference of the polygon points is given, this is used.
642        If not;
643        The x,y info is assumed to be Easting and Northing, absolute,
644        for the meshes zone.
645
646        polygon a list of points, in meters that describe the polygon
647             (e.g. [[x1,y1],[x2,y2],...]
648        tags (e.g.{'wall':[0,1,3],'ocean':[2]})
649
650        This returns the region instance (if a max_triangle_area is given),
651        so if the user whats to modify it they can.     
652        """
653        if max_triangle_area is None:
654            create_region = False
655        else:
656            create_region = True
657           
658        region = self._add_area_from_polygon(polygon, tags=tags,
659                                             geo_reference=geo_reference,
660                                             region=create_region)
661        if max_triangle_area is not None:
662            region.setMaxArea(max_triangle_area)
663        return region
664   
665       
666    def _add_area_from_polygon(self, polygon, tags=None,
667                               geo_reference=None,
668                               hole=False,
669                               region=False):
670        """
671        Add a polygon with tags to the current mesh, as a region.
672        The maxArea of the region can be specified.
673
674        If a geo_reference of the polygon points is given, this is used.
675        If not;
676        The x,y info is assumed to be Easting and Northing, absolute,
677        for the meshes zone.
678
679        polygon a list of points, in meters that describe the polygon
680             (e.g. [[x1,y1],[x2,y2],...]
681        tags (e.g.{'wall':[0,1,3],'ocean':[2]})
682
683        This returns the region instance, so if the user whats to modify
684        it they can.
685       
686        """
687        #get absolute values
688        if geo_reference is not None:
689            polygon = geo_reference.get_absolute(polygon)
690        # polygon is now absolute
691        #print "polygon  should be absolute",polygon
692       
693        #create points, segs and tags
694        region_dict = {}
695        region_dict['points'] = polygon
696       
697        #Create segments
698        #E.g. [[0,1], [1,2], [2,3], [3,0]]
699        #from polygon
700        #[0,1,2,3]
701        segments = []
702        N = len(polygon)
703        for i in range(N):
704            lo = i
705            hi = (lo + 1) % N
706            segments.append( [lo, hi] ) 
707        region_dict['segments'] = segments
708
709
710        #Create tags
711        #E.g. ['wall', 'wall', 'ocean', 'wall']
712        # from a dic
713        #{'wall':[0,1,3],'ocean':[2]}
714        segment_tags = ['']*N
715        if tags is not None:
716            for key in tags:
717                indices = tags[key]
718                for i in indices:
719                    segment_tags[i] = key
720        region_dict['segment_tags'] = segment_tags
721   
722        self.addVertsSegs(region_dict) #this is passing absolute values
723
724        if region is True:
725            #get inner point - absolute values
726            inner_point = point_in_polygon(polygon)
727            inner = self.add_region(inner_point[0], inner_point[1],
728                                    geo_reference=None) 
729        elif hole is True:
730            #get inner point - absolute values
731            inner_point = point_in_polygon(polygon)
732            inner = self.add_hole(inner_point[0], inner_point[1],
733                                    geo_reference=None)
734        else:
735            inner = None
736           
737        return inner
738
739    def add_circle(self, center, radius, segment_count=100,
740                   center_geo_reference=None, tag = None,
741                   region=False, hole=False):
742        """
743        center is a point, in absulute or relative to center_geo_ref
744        radius is the radius of the circle
745        segment_count is the number of segments used to represent the circle
746        tag is a string name, given to the segments.
747        If region is True a region object is returned.
748        If hole is True a hole object is returned.
749           (Don't have them both True.)
750
751           
752        """
753        pass
754        # convert center and radius to a polygon
755        ### pi = math.pi
756        ### num_of_cuts = 100
757        cuts = []
758        factor = 2* math.pi/segment_count
759        for cut in range(segment_count):
760             cuts.append(cut*factor)
761
762        polygon = []
763        for cut in cuts:
764           
765            x = center[0] + radius * math.cos(cut)
766            y = center[1] + radius * math.sin(cut)
767            polygon.append([x,y])
768        # build the tags
769        tags = {tag:range(segment_count)}
770       
771        return self._add_area_from_polygon(polygon, tags=tags,
772                                           region=region, hole=hole,
773                                           geo_reference=center_geo_reference)
774
775    def auto_set_geo_reference(self):
776        """
777        Automatically set the georeference, based in the minimum x and y
778        user vertex values.
779
780        Not to be used with the graphical interface
781       
782        Not implemented.
783        Don't implement now.  using the new georeferenced points class
784        will change this?
785        """
786        #to do
787        # find the lower left hand corner
788        [xmin, ymin, xmax, ymax] = self.boxsize()
789
790        # set all points to that lower left hand corner.
791        # use change_geo_reference
792        new_geo = Geo_reference(self.geo_reference.get_zone(), xmin, ymin)
793        self.change_geo_reference(new_geo)
794       
795    def change_geo_reference(self, new_geo_reference):
796        """
797        Change from one geo_reference to another.
798        Not to be used with the graphical interface
799        """
800        # FIXME
801        # change georeference of;
802        #self.userVertices = self.geo_reference.change_points_geo_ref( \
803        #self.userVertices)
804        #self.holes = self.geo_reference.change_points_geo_ref(self.holes)
805        #self.regions = self.geo_reference.change_points_geo_ref(self.regions)
806        # The above will not work.
807        # since userVertices (etc) is a list of point objects,
808        #not a list of lists.
809        # add a method to the points class to fix this up.
810       
811    def addVertsSegs(self, outlineDict):
812        """
813        Add out-line (user Mesh) attributes given a dictionary of the lists
814        points: [(x1,y1),(x2,y2),...] (Tuples of doubles) 
815        segments: [(point1,point2),(p3,p4),...] (Tuples of integers)
816        segment_tags: [S1Tag, S2Tag, ...] (list of strings)
817
818        Assume the values are in Eastings and Northings, with no reference
819        point. eg absolute
820        """
821        if not outlineDict.has_key('segment_tags'):
822            outlineDict['segment_tags'] = []
823            for i in range(len(outlineDict['segments'])):
824                outlineDict['segment_tags'].append('')
825        #print "outlineDict['segment_tags']",outlineDict['segment_tags']
826        #print "outlineDict['points']",outlineDict['points']
827        #print "outlineDict['segments']",outlineDict['segments']
828       
829        localUserVertices = []
830        #index = 0
831        for point in outlineDict['points']:
832            v=Vertex(point[0]-self.geo_reference.xllcorner,
833                     point[1]-self.geo_reference.yllcorner)
834            #v.index = index
835            #index +=1
836            self.userVertices.append(v)
837            localUserVertices.append(v)
838           
839        #index = 0
840        for seg,seg_tag in map(None,outlineDict['segments'],
841                       outlineDict['segment_tags']):
842            segObject = Segment( localUserVertices[int(seg[0])],
843                           localUserVertices[int(seg[1])] )
844            if not seg_tag == '':
845                segObject.set_tag(seg_tag)
846            #segObject.index = index
847            #index +=1
848            self.userSegments.append(segObject)
849            #DSG!!!
850       
851    def get_triangle_count(self):
852        return len(self.getTriangulation())
853       
854    def getUserVertices(self):
855        return self.userVertices
856   
857    def getUserSegments(self):
858        allSegments = self.userSegments + self.alphaUserSegments
859        #print "self.userSegments",self.userSegments
860        #print "self.alphaUserSegments",self.alphaUserSegments
861        #print "allSegments",allSegments
862        return allSegments
863   
864    def deleteUserSegments(self,seg):
865        if self.userSegments.count(seg) == 0:
866            self.alphaUserSegments.remove(seg)
867            pass
868        else:
869            self.userSegments.remove(seg)
870           
871    def clearUserSegments(self):
872        self.userSegments = []
873        self.alphaUserSegments = []
874               
875    def getTriangulation(self):
876        return self.meshTriangles
877   
878    def getMeshVertices(self):
879        return self.meshVertices
880 
881    def getMeshSegments(self):
882        return self.meshSegments
883   
884    def getHoles(self):
885        return self.holes
886   
887    def getRegions(self):
888        return self.regions
889   
890    def isTriangulation(self):
891        if self.meshVertices == []:
892            return False 
893        else:
894            return True
895   
896    def addUserSegment(self, v1,v2):
897        """
898        PRECON: A segment between the two vertices is not already present.
899        Check by calling isUserSegmentNew before calling this function.
900       
901        """
902        s=Segment( v1,v2)
903        self.userSegments.append(s)
904        return s
905       
906    def generate_mesh(self,
907                      maximum_triangle_area=None,
908                      minimum_triangle_angle=28.0,
909                      verbose=False):
910        if verbose is True:
911            silent = ''
912        else:
913            silent = 'Q'
914        self.generateMesh(mode = silent +"pzq"+str(minimum_triangle_angle)
915                                  +"a"+str(maximum_triangle_area)
916                                  +"a")
917        #The last a is so areas for regions will be used
918       
919    def generateMesh(self, mode = None, maxArea = None, minAngle= None,
920                     isRegionalMaxAreas = True):
921        """
922        Based on the current user vaules, holes and regions
923        generate a new mesh
924        mode is a string that sets conditions on the mesh generations
925        see triangle_instructions.txt for a definition of the commands
926        PreCondition: maxArea is a double
927        """
928        #print "mode ",mode
929        if mode == None:
930            self.mode = ""
931        else:
932            self.mode = mode
933       
934        if self.mode.find('p') < 0:
935            self.mode += 'p' #p - read a planar straight line graph.
936            #there must be segments to use this switch
937            # TODO throw an aception if there aren't seg's
938            # it's more comlex than this.  eg holes
939        if self.mode.find('z') < 0:
940            self.mode += 'z' # z - Number all items starting from zero
941                             # (rather than one)
942        if self.mode.find('n'):
943            self.mode += 'n' # n - output a list of neighboring triangles
944        if self.mode.find('A') < 0:
945            self.mode += 'A' # A - output region attribute list for triangles
946
947        if not self.mode.find('V') and not self.mode.find('Q'):
948            self.mode += 'V' # V - output info about what Triangle is doing
949       
950        if self.mode.find('q') < 0 and minAngle is not None:
951            #   print "**********8minAngle******** ",minAngle
952            min_angle = 'q' + str(minAngle)
953            self.mode += min_angle # z - Number all items starting from zero
954                             # (rather than one)
955        if maxArea != None:
956            self.mode += 'a' + str(maxArea)
957        #FIXME (DSG-DSG) This isn't explained.  I don't think it works
958        # well with  maxArea = None as well.
959        if isRegionalMaxAreas:
960            self.mode += 'a'
961         
962        meshDict = self.Mesh2triangList()
963        #print "*************************!@!@ This is going to triangle   !@!@"
964        #print meshDict
965        #print "************************!@!@ This is going to triangle   !@!@"
966
967        #print "meshDict['segmenttaglist']", meshDict['segmenttaglist']
968        [meshDict['segmenttaglist'],
969         segconverter] =  segment_strings2ints(meshDict['segmenttaglist'],
970                                             initialconversions)
971        #print "regionlist",meshDict['regionlist']
972        [meshDict['regionlist'],
973         regionconverter] =  region_strings2ints(meshDict['regionlist'])
974        #print "%%%%%%%%%%%%%%%%%%%%%%%%%%%regionlist",meshDict['regionlist']
975        #print "meshDict['segmenttaglist']", meshDict['segmenttaglist'
976        generatedMesh = triang.genMesh(
977                              meshDict['pointlist'],
978                              meshDict['segmentlist'],
979                              meshDict['holelist'],
980                              meshDict['regionlist'],
981                              meshDict['pointattributelist'],
982                              meshDict['segmenttaglist'],
983                              [],  # since the trianglelist isn't used
984                              self.mode)
985        #print "%%%%%%%%%%%%%%%%%%%%%%%%%%%generated",generatedMesh
986        generatedMesh['generatedsegmentmarkerlist'] = \
987             segment_ints2strings(generatedMesh['generatedsegmentmarkerlist'],
988                                  segconverter)
989        #print "processed gen",generatedMesh['generatedsegmentmarkerlist']
990        generatedMesh['generatedtriangleattributelist'] = \
991         region_ints2strings(generatedMesh['generatedtriangleattributelist'],
992                                  regionconverter)
993
994
995        if len(generatedMesh['generatedpointattributelist'][0])==0:
996            self.attributeTitles = []
997        generatedMesh['generatedpointattributetitlelist']= \
998                                            self.attributeTitles
999        #print "################  FROM TRIANGLE"
1000        #print "generatedMesh",generatedMesh
1001        #print "##################"
1002       
1003        self.setTriangulation(generatedMesh)
1004   
1005    def clearTriangulation(self):
1006
1007        #Clear the current generated mesh values
1008        self.meshTriangles=[] 
1009        self.meshSegments=[]
1010        self.meshVertices=[]
1011
1012    def removeDuplicatedUserVertices(self):
1013        """Pre-condition: There are no user segments
1014        This function will keep the first duplicate
1015        """
1016        assert self.getUserSegments() == []
1017        self.userVertices, counter =  self.removeDuplicatedVertices(
1018            self.userVertices)
1019        return counter
1020   
1021    def removeDuplicatedVertices(self, Vertices):
1022        """
1023        This function will keep the first duplicate, remove all others
1024        Precondition: Each vertex has a dupindex, which is the list
1025        index.
1026
1027        Note: this removes vertices that have the same x,y values,
1028        not duplicate instances in the Vertices list.
1029        """
1030        remove = []
1031        index = 0
1032        for v in Vertices:
1033            v.dupindex = index
1034            index += 1
1035        t = list(Vertices)
1036        t.sort(Point.cmp_xy)
1037   
1038        length = len(t)
1039        behind = 0
1040        ahead  = 1
1041        counter = 0
1042        while ahead < length:
1043            b = t[behind]
1044            ah = t[ahead]
1045            if (b.y == ah.y and b.x == ah.x):
1046                remove.append(ah.dupindex) 
1047            behind += 1
1048            ahead += 1
1049
1050        # remove the duplicate vertices
1051        remove.sort()
1052        remove.reverse()
1053        for i in remove:
1054            Vertices.pop(i)
1055            pass
1056
1057        #Remove the attribute that this function added
1058        for v in Vertices:
1059            del v.dupindex
1060        return Vertices,counter
1061   
1062    def thinoutVertices(self, delta):
1063        """Pre-condition: There are no user segments
1064        This function will keep the first duplicate
1065        """
1066        assert self.getUserSegments() == []
1067        #t = self.userVertices
1068        #self.userVertices =[]
1069        boxedVertices = {}
1070        thinnedUserVertices =[]
1071        delta = round(delta,1)
1072       
1073        for v in self.userVertices :
1074            # tag is the center of the boxes
1075            tag = (round(v.x/delta,0)*delta,round(v.y/delta,0)*delta)
1076            #this creates a dict of lists of faces, indexed by tag
1077            boxedVertices.setdefault(tag,[]).append(v)
1078
1079        for [tag,verts] in boxedVertices.items():
1080            min = delta
1081            bestVert = None
1082            tagVert = Vertex(tag[0],tag[1])
1083            for v in verts:
1084                dist = v.DistanceToPoint(tagVert)
1085                if (dist<min):
1086                    min = dist
1087                    bestVert = v
1088            thinnedUserVertices.append(bestVert)
1089        self.userVertices =thinnedUserVertices
1090       
1091           
1092    def isUserSegmentNew(self, v1,v2):
1093        identicalSegs= [x for x in self.getUserSegments() \
1094                        if (x.vertices[0] == v1 and x.vertices[1] == v2)
1095        or (x.vertices[0] == v2 and x.vertices[1] == v1) ]
1096       
1097        return len(identicalSegs) == 0
1098
1099       
1100    def deleteSegsOfVertex(self, delVertex):
1101        """
1102        Delete this vertex and any segments that connect to it.
1103        """
1104        #Find segments that connect to delVertex
1105        deletedSegments = []
1106        for seg in self.getUserSegments():
1107            if (delVertex in seg.vertices):
1108                deletedSegments.append(seg)
1109        # Delete segments that connect to delVertex
1110        for seg in deletedSegments:
1111            self.deleteUserSegments(seg)
1112        return deletedSegments
1113
1114   
1115    def deleteMeshObject(self, MeshObject):
1116        """
1117        Returns a list of all objects that were removed
1118        """
1119        deletedObs = []
1120        if isinstance(MeshObject, Vertex ):
1121            deletedObs = self.deleteSegsOfVertex(MeshObject)
1122            deletedObs.append(MeshObject)
1123            self.userVertices.remove(MeshObject)
1124        elif isinstance(MeshObject, Segment):
1125            deletedObs.append(MeshObject)
1126            self.deleteUserSegments(MeshObject)
1127        elif isinstance(MeshObject, Hole):
1128            deletedObs.append(MeshObject)
1129            self.holes.remove(MeshObject)
1130        elif isinstance(MeshObject, Region):
1131            deletedObs.append(MeshObject)
1132            self.regions.remove(MeshObject)         
1133        return deletedObs
1134                                                 
1135    def Mesh2triangList(self, userVertices=None,
1136                        userSegments=None,
1137                        holes=None,
1138                        regions=None):
1139        """
1140        Convert the Mesh to a dictionary of the lists needed for the
1141        triang module
1142        points list: [(x1,y1),(x2,y2),...] (Tuples of doubles)
1143        pointattributelist: [(a11,a12,...),(a21,a22),...] (Tuples of doubles)
1144        segment list: [(point1,point2),(p3,p4),...] (Tuples of integers)
1145        hole list: [(x1,y1),...](Tuples of doubles, one inside each hole)
1146        regionlist: [ (x1,y1,tag, max area),...] (Tuple of 3-4 doubles)
1147       
1148        Note, this adds an index attribute to the user Vertex objects.
1149
1150        Used to produce output to triangle
1151        """
1152        if userVertices is None:
1153            userVertices = self.getUserVertices()
1154        if userSegments is None:
1155            userSegments = self.getUserSegments()
1156        if holes is None:
1157            holes = self.getHoles()
1158        if regions is None:
1159            regions = self.getRegions()
1160           
1161        meshDict = {}
1162       
1163        pointlist=[]
1164        pointattributelist=[]
1165        index = 0
1166        for vertex in userVertices:
1167            vertex.index = index
1168            pointlist.append((vertex.x,vertex.y))
1169            pointattributelist.append((vertex.attributes))
1170           
1171            index += 1
1172        meshDict['pointlist'] = pointlist
1173        meshDict['pointattributelist'] = pointattributelist
1174
1175        segmentlist=[]
1176        segmenttaglist=[]
1177        for seg in userSegments:
1178            segmentlist.append((seg.vertices[0].index,seg.vertices[1].index))
1179            segmenttaglist.append(seg.tag)
1180        meshDict['segmentlist'] =segmentlist
1181        meshDict['segmenttaglist'] =segmenttaglist
1182       
1183        holelist=[]
1184        for hole in holes:
1185            holelist.append((hole.x,hole.y)) 
1186        meshDict['holelist'] = holelist
1187       
1188        regionlist=[]
1189        for region in regions:
1190            if (region.getMaxArea() != None): 
1191                regionlist.append((region.x,region.y,region.getTag(),
1192                               region.getMaxArea()))
1193            else:
1194                regionlist.append((region.x,region.y,region.getTag()))
1195        meshDict['regionlist'] = regionlist
1196        #print "*(*("
1197        #print meshDict
1198        #print meshDict['regionlist']
1199        #print "*(*("
1200        return meshDict
1201                                               
1202    def Mesh2MeshList(self):
1203        """
1204        Convert the Mesh to a dictionary of lists describing the
1205        triangulation variables;
1206        generated point list: [(x1,y1),(x2,y2),...] (Tuples of doubles)
1207        generated point attribute list: [(a11,a12,...),(a21,a22),...]
1208            (Tuples of doubles)
1209        generated point attribute title list:[A1Title, A2Title ...]
1210            (list of strings)
1211        generated segment list: [(point1,point2),(p3,p4),...]
1212            (Tuples of integers)
1213        generated segment tag list: [tag,tag,...] list of strings
1214
1215        generated triangle list: [(p1,p2,p3), (p4,p5,p6),....] tuple of points
1216
1217        generated triangle attribute list: [s1,s2,...] list of strings
1218
1219        generated triangle neighbor list: [(t1,t2,t3), (t4,t5,t6),....]
1220            tuple of triangles
1221       
1222        Used to produce .tsh file
1223        """
1224       
1225        meshDict = {}       
1226        pointlist=[]
1227        pointattributelist=[]
1228
1229
1230        self.maxVertexIndex=0
1231        for vertex in self.meshVertices:
1232            vertex.index = self.maxVertexIndex
1233            pointlist.append((vertex.x,vertex.y))
1234            pointattributelist.append((vertex.attributes))           
1235            self.maxVertexIndex += 1
1236
1237        meshDict['generatedpointlist'] = pointlist
1238        meshDict['generatedpointattributelist'] = pointattributelist
1239        meshDict['generatedpointattributetitlelist'] = self.attributeTitles
1240        #segments
1241        segmentlist=[]
1242        segmenttaglist=[]
1243        for seg in self.meshSegments:
1244            segmentlist.append((seg.vertices[0].index,seg.vertices[1].index))
1245            segmenttaglist.append(seg.tag)
1246        meshDict['generatedsegmentlist'] =segmentlist
1247        meshDict['generatedsegmenttaglist'] =segmenttaglist
1248
1249        # Make sure that the indexation is correct
1250        index = 0
1251        for tri in self.meshTriangles:
1252            tri.index = index           
1253            index += 1
1254
1255        trianglelist = []
1256        triangleattributelist = []
1257        triangleneighborlist = []
1258        for tri in self.meshTriangles: 
1259            trianglelist.append((tri.vertices[0].index,tri.vertices[1].index,
1260                                 tri.vertices[2].index)) 
1261            triangleattributelist.append([tri.attribute])
1262            neighborlist = [-1,-1,-1]
1263            for neighbor,index in map(None,tri.neighbors,
1264                                      range(len(tri.neighbors))):
1265                if neighbor:   
1266                    neighborlist[index] = neighbor.index
1267            triangleneighborlist.append(neighborlist)
1268       
1269        meshDict['generatedtrianglelist'] = trianglelist
1270        meshDict['generatedtriangleattributelist'] = triangleattributelist
1271        meshDict['generatedtriangleneighborlist'] = triangleneighborlist
1272       
1273        #print "mesh.Mesh2MeshList*)*)"
1274        #print meshDict
1275        #print "mesh.Mesh2MeshList*)*)"
1276
1277        return meshDict
1278
1279                               
1280    def Mesh2MeshDic(self):
1281        """
1282        Convert the user and generated info of a mesh to a dictionary
1283        structure
1284        """
1285        dic = self.Mesh2triangList()
1286        dic_mesh = self.Mesh2MeshList()
1287        for element in dic_mesh.keys():
1288            dic[element] = dic_mesh[element]
1289        return dic
1290   
1291    def setTriangulation(self, genDict):
1292        """
1293        Set the mesh attributes given a dictionary of the lists
1294        returned from the triang module       
1295        generated point list: [(x1,y1),(x2,y2),...] (Tuples of doubles) 
1296        generated point attribute list:[(P1att1,P1attt2, ...),
1297            (P2att1,P2attt2,...),...]
1298        generated point attribute title list:[A1Title, A2Title ...]
1299            (list of strings)
1300        generated segment list: [(point1,point2),(p3,p4),...]
1301            (Tuples of integers)
1302        generated segment marker list: [S1Tag, S2Tag, ...] (list of ints)
1303        triangle list:  [(point1,point2, point3),(p5,p4, p1),...]
1304            (Tuples of integers)
1305        triangle neighbor list: [(triangle1,triangle2, triangle3),
1306            (t5,t4, t1),...] (Tuples of integers)
1307            -1 means there's no triangle neighbor
1308        triangle attribute list: [(T1att), (T2att), ...]
1309            (list of a list of strings)
1310        lone point list:[point1, ...] (list of integers)
1311        """
1312        #Clear the current generated mesh values
1313        self.meshTriangles=[]
1314        self.attributeTitles=[]
1315        self.meshSegments=[]
1316        self.meshVertices=[]
1317
1318        #print "mesh.setTriangulation@#@#@#"
1319        #print genDict
1320        #print "@#@#@#"
1321
1322        self.maxVertexIndex = 0
1323        for point in genDict['generatedpointlist']:
1324            v=Vertex(point[0], point[1])
1325            v.index =  self.maxVertexIndex
1326            self.maxVertexIndex +=1
1327            self.meshVertices.append(v)
1328
1329        self.attributeTitles = genDict['generatedpointattributetitlelist']
1330
1331        index = 0
1332        for seg,marker in map(None,genDict['generatedsegmentlist'],
1333                              genDict['generatedsegmentmarkerlist']):
1334            segObject = Segment( self.meshVertices[int(seg[0])],
1335                           self.meshVertices[int(seg[1])], tag = marker )
1336            segObject.index = index
1337            index +=1
1338            self.meshSegments.append(segObject)
1339
1340        index = 0
1341        for triangle in genDict['generatedtrianglelist']:
1342            tObject =Triangle( self.meshVertices[int(triangle[0])],
1343                        self.meshVertices[int(triangle[1])],
1344                        self.meshVertices[int(triangle[2])] )
1345            tObject.index = index
1346            index +=1
1347            self.meshTriangles.append(tObject)
1348
1349        index = 0
1350        for att in genDict['generatedtriangleattributelist']:
1351            if att == []:
1352                self.meshTriangles[index].setAttribute("")
1353            else:
1354                self.meshTriangles[index].setAttribute(att[0])
1355            index += 1
1356           
1357        index = 0
1358        for att in genDict['generatedpointattributelist']:
1359            if att == None:
1360                self.meshVertices[index].setAttributes([])
1361            else:
1362                self.meshVertices[index].setAttributes(att)
1363            index += 1
1364   
1365        index = 0
1366        for triangle in genDict['generatedtriangleneighborlist']:
1367            # Build a list of triangle object neighbors
1368            ObjectNeighbor = []
1369            for neighbor in triangle:
1370                if ( neighbor != -1):
1371                    ObjectNeighbor.append(self.meshTriangles[int(neighbor)])
1372                else:
1373                    ObjectNeighbor.append(None)
1374            self.meshTriangles[index].setNeighbors(ObjectNeighbor[0],
1375                                                   ObjectNeighbor[1],
1376                                                   ObjectNeighbor[2])
1377            index += 1
1378        genDict['lonepointlist'].sort() 
1379        genDict['lonepointlist'].reverse()       
1380        for loner in genDict['lonepointlist']:
1381            # Remove the loner vertex
1382            #print "Removing the loner", loner
1383            self.meshVertices.pop(loner)
1384
1385    def setMesh(self, genDict):
1386        """
1387        Set the user Mesh attributes given a dictionary of the lists
1388        point list: [(x1,y1),(x2,y2),...] (Tuples of doubles) 
1389        point attribute list:[(P1att1,P1attt2, ...),(P2att1,P2attt2,...),...]
1390        segment list: [(point1,point2),(p3,p4),...] (Tuples of integers)
1391        segment tag list: [S1Tag, S2Tag, ...] (list of ints)
1392        region list: [(x1,y1),(x2,y2),...] (Tuples of doubles)
1393        region attribute list: ["","reservoir",""] list of strings
1394        region max area list:[real, None, Real,...] list of None and reals
1395       
1396        mesh is an instance of a mesh object
1397        """
1398        #Clear the current user mesh values
1399        self.clearUserSegments()
1400        self.userVertices=[]
1401        self.Holes=[]
1402        self.Regions=[]
1403
1404        #print "mesh.setMesh@#@#@#"
1405        #print genDict
1406        #print "@#@#@#"
1407       
1408        #index = 0
1409        for point in genDict['pointlist']:
1410            v=Vertex(point[0], point[1])
1411            #v.index = index
1412            #index +=1
1413            self.userVertices.append(v)
1414
1415        #index = 0
1416        for seg,tag in map(None,genDict['segmentlist'],
1417                           genDict['segmenttaglist']):
1418            segObject = Segment( self.userVertices[int(seg[0])],
1419                           self.userVertices[int(seg[1])], tag = tag )
1420            #segObject.index = index
1421            #index +=1
1422            self.userSegments.append(segObject)
1423
1424# Remove the loading of attribute info.
1425# Have attribute info added using least_squares in pyvolution
1426#         index = 0
1427#         for att in genDict['pointattributelist']:
1428#             if att == None:
1429#                 self.userVertices[index].setAttributes([])
1430#             else:
1431#                 self.userVertices[index].setAttributes(att)
1432#            index += 1
1433       
1434        #index = 0
1435        for point in genDict['holelist']:
1436            h=Hole(point[0], point[1])
1437            #h.index = index
1438            #index +=1
1439            self.holes.append(h)
1440
1441        #index = 0
1442        for reg,att,maxArea in map(None,
1443                                   genDict['regionlist'],
1444                                   genDict['regionattributelist'],
1445                                   genDict['regionmaxarealist']):
1446            Object = Region( reg[0],
1447                             reg[1],
1448                             tag = att,
1449                             maxArea = maxArea)
1450            #Object.index = index
1451            #index +=1
1452            self.regions.append(Object)
1453           
1454    def TestautoSegment(self):
1455        newsegs = []
1456        s1 = Segment(self.userVertices[0],
1457                               self.userVertices[1])
1458        s2 = Segment(self.userVertices[0],
1459                               self.userVertices[2])
1460        s3 = Segment(self.userVertices[2],
1461                               self.userVertices[1])
1462        if self.isUserSegmentNew(s1.vertices[0],s1.vertices[1]):
1463            newsegs.append(s1)
1464        if self.isUserSegmentNew(s2.vertices[0],s2.vertices[1]):
1465            newsegs.append(s2)
1466        if self.isUserSegmentNew(s3.vertices[0],s3.vertices[1]):
1467            newsegs.append(s3)
1468        #DSG!!!
1469        self.userSegments.extend(newsegs)
1470        return newsegs
1471
1472   
1473    def savePickle(self, currentName):
1474        fd = open(currentName, 'w')
1475        pickle.dump(self,fd)
1476        fd.close()
1477
1478    def autoSegmentHull(self):
1479        """
1480        initially work by running an executable
1481        Later compile the c code with a python wrapper.
1482
1483        Precon: There must be 3 or more vertices in the userVertices structure
1484        """
1485        newsegs = []
1486        inputfile = 'hull_in.txt'
1487        outputfile = inputfile + '-alf'
1488        #write vertices to file
1489        fd = open(inputfile,'w')
1490        for v in self.userVertices:
1491            fd.write(str(v.x))
1492            fd.write(' ')
1493            fd.write(str(v.y))
1494            fd.write('\n')
1495        fd.close()
1496       
1497        #run hull executable
1498        #warning need to compile hull for the current operating system
1499        command = 'hull.exe -A -i ' + inputfile
1500        os.system(command)
1501       
1502        #read results into this object
1503        fd = open(outputfile)
1504        lines = fd.readlines()
1505        fd.close()
1506        #print "(*(*(*("
1507        #print lines
1508        #print "(*(*(*("
1509        lines.pop(0) #remove the first (title) line
1510        for line in lines:
1511            vertindexs = line.split()
1512            #print 'int(vertindexs[0])', int(vertindexs[0])
1513            #print 'int(vertindexs[1])', int(vertindexs[1])
1514            #print 'self.userVertices[int(vertindexs[0])]' ,self.userVertices[int(vertindexs[0])]
1515            #print 'self.userVertices[int(vertindexs[1])]' ,self.userVertices[int(vertindexs[1])]
1516            v1 = self.userVertices[int(vertindexs[0])]
1517            v2 = self.userVertices[int(vertindexs[1])]
1518           
1519            if self.isUserSegmentNew(v1,v2):
1520                newseg = Segment(v1, v2)
1521                newsegs.append(newseg)
1522            #DSG!!!
1523        self.userSegments.extend(newsegs)
1524        return newsegs
1525    def autoSegmentFilter(self,raw_boundary=True,
1526                    remove_holes=False,
1527                    smooth_indents=False,
1528                    expand_pinch=False):
1529        """
1530        Precon: There is a self.shape
1531        """
1532        #FIXME remove the precon.  Internally check
1533        return self._boundary2mesh(raw_boundary=raw_boundary,
1534                    remove_holes=remove_holes,
1535                    smooth_indents=smooth_indents,
1536                    expand_pinch=expand_pinch)
1537   
1538       
1539   
1540    def autoSegment(self, alpha = None,
1541                    raw_boundary=True,
1542                    remove_holes=False,
1543                    smooth_indents=False,
1544                    expand_pinch=False): 
1545        """
1546        Precon: There must be 3 or more vertices in the userVertices structure
1547        """
1548        self._createBoundary(alpha=alpha)
1549        return self._boundary2mesh(raw_boundary=raw_boundary,
1550                    remove_holes=remove_holes,
1551                    smooth_indents=smooth_indents,
1552                    expand_pinch=expand_pinch)
1553
1554    def _createBoundary(self,alpha=None):
1555        """
1556        """
1557        points=[]
1558        for vertex in self.getUserVertices():
1559            points.append((vertex.x,vertex.y))
1560        self.shape = alpha_shape.alpha_shape.Alpha_Shape(points, alpha = alpha)
1561
1562
1563    def _boundary2mesh(self, raw_boundary=True,
1564                    remove_holes=False,
1565                    smooth_indents=False,
1566                    expand_pinch=False):
1567        """
1568        Precon there must be a shape object.
1569        """
1570        self.shape.set_boundary_type(raw_boundary=raw_boundary,
1571                                 remove_holes=remove_holes,
1572                                 smooth_indents=smooth_indents,
1573                                 expand_pinch=expand_pinch)
1574        boundary_segs = self.shape.get_boundary()
1575        #print "boundary_segs",boundary_segs
1576        segs2delete = self.alphaUserSegments
1577        #FIXME(DSG-DSG) this algorithm needs comments
1578        new_segs = {}
1579        #alpha_segs = []
1580        #user_segs = []
1581        for seg in boundary_segs:
1582            v1 = self.userVertices[int(seg[0])]
1583            v2 = self.userVertices[int(seg[1])]
1584            boundary_seg = Segment(v1, v2)
1585            new_segs[(v1,v2)] = boundary_seg
1586
1587        for user_seg in self.userSegments:
1588            if new_segs.has_key((user_seg.vertices[0],
1589                                user_seg.vertices[1])):
1590                del new_segs[user_seg.vertices[0],
1591                                user_seg.vertices[1]]
1592            elif new_segs.has_key((user_seg.vertices[1],
1593                                user_seg.vertices[0])):
1594                del new_segs[user_seg.vertices[1],
1595                                user_seg.vertices[0]]
1596               
1597        optimum_alpha = self.shape.get_alpha()
1598        alpha_segs_no_user_segs  = new_segs.values()
1599        self.alphaUserSegments = alpha_segs_no_user_segs
1600        return alpha_segs_no_user_segs, segs2delete, optimum_alpha
1601   
1602    def _boundary2mesh_old(self, raw_boundary=True,
1603                    remove_holes=False,
1604                    smooth_indents=False,
1605                    expand_pinch=False):
1606        """
1607        Precon there must be a shape object.
1608        """
1609        self.shape.set_boundary_type(raw_boundary=raw_boundary,
1610                                 remove_holes=remove_holes,
1611                                 smooth_indents=smooth_indents,
1612                                 expand_pinch=expand_pinch)
1613        boundary_segs = self.shape.get_boundary()
1614        #print "boundary_segs",boundary_segs
1615        segs2delete = self.alphaUserSegments
1616
1617        #FIXME(DSG-DSG) this algorithm needs comments
1618        #FIXME(DSG-DSG) can it be sped up?  It's slow
1619        new_segs = []
1620        alpha_segs = []
1621        user_segs = []
1622        for seg in boundary_segs:
1623            v1 = self.userVertices[int(seg[0])]
1624            v2 = self.userVertices[int(seg[1])]
1625            alpha_seg = self.representedAlphaUserSegment(v1, v2)
1626            user_seg = self.representedUserSegment(v1, v2)
1627            #DSG!!!
1628            assert not(not (alpha_seg == None) and not (user_seg == None))
1629            if not alpha_seg == None:
1630                alpha_segs.append(alpha_seg)
1631            elif not user_seg  == None:
1632                user_segs.append(user_seg)
1633            else:
1634                unique_seg = Segment(v1, v2)
1635                new_segs.append(unique_seg)
1636               
1637            for seg in alpha_segs:
1638                try:
1639                    segs2delete.remove(seg)
1640                except:
1641                    pass
1642       
1643        self.alphaUserSegments = []
1644        self.alphaUserSegments.extend(new_segs)
1645        self.alphaUserSegments.extend(alpha_segs)
1646
1647        optimum_alpha = self.shape.get_alpha()
1648        # need to draw newsegs
1649        return new_segs, segs2delete, optimum_alpha
1650   
1651    def representedAlphaUserSegment(self, v1,v2):
1652        identicalSegs= [x for x in self.alphaUserSegments \
1653                        if (x.vertices[0] == v1 and x.vertices[1] == v2)
1654        or (x.vertices[0] == v2 and x.vertices[1] == v1) ]
1655
1656        if identicalSegs == []:
1657            return None
1658        else:
1659            # Only return the first one.
1660            return identicalSegs[0]
1661   
1662    def representedUserSegment(self, v1,v2):
1663        identicalSegs= [x for x in self.userSegments \
1664                        if (x.vertices[0] == v1 and x.vertices[1] == v2)
1665        or (x.vertices[0] == v2 and x.vertices[1] == v1) ]
1666
1667        if identicalSegs == []:
1668            return None
1669        else:
1670            # Only return the first one.
1671            return identicalSegs[0]
1672       
1673    def joinVertices(self):
1674        """
1675        Return list of segments connecting all userVertices
1676        in the order they were given
1677       
1678        Precon: There must be 3 or more vertices in the userVertices structure
1679        """
1680
1681        newsegs = []
1682       
1683        v1 = self.userVertices[0]
1684        for v2 in self.userVertices[1:]:
1685            if self.isUserSegmentNew(v1,v2):           
1686                newseg = Segment(v1, v2)
1687                newsegs.append(newseg)
1688            v1 = v2
1689
1690        #Connect last point to the first
1691        v2 = self.userVertices[0]       
1692        if self.isUserSegmentNew(v1,v2):           
1693                newseg = Segment(v1, v2)
1694                newsegs.append(newseg)
1695           
1696
1697        #Update list of user segments       
1698        #DSG!!!
1699        self.userSegments.extend(newsegs)               
1700        return newsegs
1701   
1702    def normaliseMesh(self,scale, offset, height_scale):
1703        [xmin, ymin, xmax, ymax] = self.boxsize()
1704        [attmin0, attmax0] = self.maxMinVertAtt(0)
1705        #print "[attmin0, attmax0]" ,[attmin0, attmax0]
1706        [attmin1, attmax1] = self.maxMinVertAtt(1)
1707        #print [xmin, ymin, xmax, ymax]
1708        xrange = xmax - xmin
1709        yrange = ymax - ymin
1710        if xrange > yrange:
1711            min,max = xmin, xmax
1712        else:
1713            min,max = ymin, ymax
1714           
1715        for obj in self.getUserVertices():
1716            obj.x = (obj.x - xmin)/(max- min)*scale + offset
1717            obj.y = (obj.y - ymin)/(max- min)*scale + offset
1718            if len(obj.attributes)  > 0 and attmin0 != attmax0:
1719                obj.attributes[0] = (obj.attributes[0]-attmin0)/ \
1720                                    (attmax0-attmin0)*height_scale
1721            if len(obj.attributes) > 1 and attmin1 != attmax1:
1722                obj.attributes[1] = (obj.attributes[1]-attmin1)/ \
1723                                    (attmax1-attmin1)*height_scale
1724           
1725        for obj in self.getMeshVertices():
1726            obj.x = (obj.x - xmin)/(max- min)*scale + offset
1727            obj.y = (obj.y - ymin)/(max- min)*scale + offset
1728            if len(obj.attributes)  > 0 and attmin0 != attmax0:
1729                obj.attributes[0] = (obj.attributes[0]-attmin0)/ \
1730                                    (attmax0-attmin0)*height_scale
1731            if len(obj.attributes) > 1 and attmin1 != attmax1:
1732                obj.attributes[1] = (obj.attributes[1]-attmin1)/ \
1733                                    (attmax1-attmin1)*height_scale
1734               
1735        for obj in self.getHoles():
1736            obj.x = (obj.x - xmin)/(max- min)*scale + offset
1737            obj.y = (obj.y - ymin)/(max- min)*scale + offset
1738        for obj in self.getRegions():
1739            obj.x = (obj.x - xmin)/(max- min)*scale + offset
1740            obj.y = (obj.y - ymin)/(max- min)*scale + offset
1741        [xmin, ymin, xmax, ymax] = self.boxsize()
1742        #print [xmin, ymin, xmax, ymax]
1743     
1744    def boxsizeVerts(self):
1745        """
1746        Returns a list of verts denoting a box or triangle that contains
1747        verts on the xmin, ymin, xmax and ymax axis.
1748        Structure: list of verts
1749        """
1750        # FIXME dsg!!! large is a hack
1751        #You want the kinds package, part of Numeric:
1752        #In [2]: import kinds
1753       
1754        #In [3]: kinds.default_float_kind.M
1755        #kinds.default_float_kind.MAX         kinds.default_float_kind.MIN
1756    #kinds.default_float_kind.MAX_10_EXP  kinds.default_float_kind.MIN_10_EXP
1757        #kinds.default_float_kind.MAX_EXP     kinds.default_float_kind.MIN_EXP
1758
1759        #In [3]: kinds.default_float_kind.MIN
1760        #Out[3]: 2.2250738585072014e-308
1761
1762        large = 1e100
1763        xmin= large
1764        xmax=-large
1765        ymin= large
1766        ymax=-large
1767        for vertex in self.userVertices:
1768            if vertex.x < xmin:
1769                xmin = vertex.x
1770                xminVert = vertex
1771            if vertex.x > xmax:
1772                xmax = vertex.x
1773                xmaxVert = vertex
1774               
1775            if vertex.y < ymin:
1776                ymin = vertex.y
1777                yminVert = vertex
1778            if vertex.y > ymax:
1779                ymax = vertex.y
1780                ymaxVert = vertex
1781        verts, count = self.removeDuplicatedVertices([xminVert,
1782                                                      xmaxVert,
1783                                                      yminVert,
1784                                                      ymaxVert])
1785         
1786        return verts
1787   
1788    def boxsize(self):
1789        """
1790        Returns a list denoting a box that contains the entire
1791        structure of vertices
1792        Structure: [xmin, ymin, xmax, ymax]
1793        """
1794        # FIXME dsg!!! large is a hack
1795        #You want the kinds package, part of Numeric:
1796        #In [2]: import kinds
1797       
1798        #In [3]: kinds.default_float_kind.M
1799        #kinds.default_float_kind.MAX         kinds.default_float_kind.MIN
1800    #kinds.default_float_kind.MAX_10_EXP  kinds.default_fltesting oat_kind.MIN_10_EXP
1801        #kinds.default_float_kind.MAX_EXP     kinds.default_float_kind.MIN_EXP
1802
1803        #In [3]: kinds.default_float_kind.MIN
1804        #Out[3]: 2.2250738585072014e-308
1805
1806        large = 1e100
1807        xmin= large
1808        xmax=-large
1809        ymin= large
1810        ymax=-large
1811        for vertex in self.userVertices:
1812            if vertex.x < xmin:
1813                xmin = vertex.x
1814            if vertex.x > xmax:
1815                xmax = vertex.x
1816               
1817            if vertex.y < ymin:
1818                ymin = vertex.y
1819            if vertex.y > ymax:
1820                ymax = vertex.y
1821        return [xmin, ymin, xmax, ymax]
1822 
1823    def maxMinVertAtt(self, iatt):
1824        """
1825        Returns a list denoting a box that contains the entire structure
1826        of vertices
1827        Structure: [xmin, ymin, xmax, ymax]
1828        """
1829        # FIXME dsg!!! large is a hacktesting
1830        #You want the kinds package, part of Numeric:
1831        #In [2]: import kinds
1832       
1833        #In [3]: kinds.default_float_kind.M
1834        #kinds.default_float_kind.MAX         kinds.default_float_kind.MIN
1835    #kinds.default_float_kind.MAX_10_EXP  kinds.default_float_kind.MIN_10_EXP
1836        #kinds.default_float_kind.MAX_EXP     kinds.default_float_kind.MIN_EXP
1837
1838        #In [3]: kinds.default_float_kind.MIN
1839        #Out[3]: 2.2250738585072014e-308
1840
1841        large = 1e100
1842        min= large
1843        max=-large
1844        for vertex in self.userVertices:
1845            if len(vertex.attributes) > iatt:
1846                if vertex.attributes[iatt] < min:
1847                    min = vertex.attributes[iatt]
1848                if vertex.attributes[iatt] > max:
1849                    max = vertex.attributes[iatt] 
1850        for vertex in self.meshVertices:
1851            if len(vertex.attributes) > iatt:
1852                if vertex.attributes[iatt] < min:
1853                    min = vertex.attributes[iatt]
1854                if vertex.attributes[iatt] > max:
1855                    max = vertex.attributes[iatt] 
1856        return [min, max]
1857   
1858    def scaleoffset(self, WIDTH, HEIGHT):
1859        """
1860        Returns a list denoting the scale and offset terms that need to be
1861        applied when converting  mesh co-ordinates onto grid co-ordinates
1862        Structure: [scale, xoffset, yoffset]
1863        """   
1864        OFFSET = 0.05*min([WIDTH, HEIGHT])
1865        [xmin, ymin, xmax, ymax] = self.boxsize()
1866        SCALE = min([0.9*WIDTH, 0.9*HEIGHT])/max([xmax-xmin, ymax-ymin])
1867       
1868        if SCALE*xmin < OFFSET:
1869            xoffset = abs(SCALE*xmin) + OFFSET
1870        if SCALE*xmax > WIDTH - OFFSET:
1871            xoffset= -(SCALE*xmax - WIDTH + OFFSET)
1872        if SCALE*ymin < OFFSET:
1873            b = abs(SCALE*ymin)+OFFSET
1874        if SCALE*ymax > HEIGHT-OFFSET:
1875            b = -(SCALE*ymax - HEIGHT + OFFSET)
1876        yoffset = HEIGHT - b
1877        return [SCALE, xoffset, yoffset]
1878
1879         
1880    def exportASCIIobj(self,ofile):
1881        """
1882        export a file, ofile, with the format
1883         lines:  v <x> <y> <first attribute>
1884        f <vertex #>  <vertex #> <vertex #> (of the triangles)
1885        """
1886        fd = open(ofile,'w')
1887        self.writeASCIIobj(fd)   
1888        fd.close()
1889
1890
1891    def writeASCIIobj(self,fd):
1892        fd.write(" # Triangulation as an obj file\n")
1893        numVert = str(len(self.meshVertices))
1894       
1895        index1 = 1 
1896        for vert in self.meshVertices:
1897            vert.index1 = index1
1898            index1 += 1
1899           
1900            fd.write("v "
1901                     + str(vert.x) + " "
1902                     + str(vert.y) + " "
1903                     + str(vert.attributes[0]) + "\n")
1904           
1905        for tri in self.meshTriangles:
1906            fd.write("f "
1907                     + str(tri.vertices[0].index1) + " " 
1908                     + str(tri.vertices[1].index1) + " " 
1909                     + str(tri.vertices[2].index1) + "\n")
1910           
1911    def exportASCIIsegmentoutlinefile(self,ofile):
1912        """Write the boundary user mesh info, eg
1913         vertices that are connected to segments,
1914         segments
1915        """
1916       
1917        verts = {}
1918        for seg in self.getUserSegments():
1919            verts[seg.vertices[0]] = seg.vertices[0]
1920            verts[seg.vertices[1]] = seg.vertices[1]
1921        meshDict = self.Mesh2IOOutlineDict(userVertices=verts.values())
1922        load_mesh.loadASCII.export_mesh_file(ofile,meshDict)
1923       
1924        # exportASCIImeshfile   - this function is used
1925    def export_mesh_file(self,ofile):
1926        """
1927        export a file, ofile, with the format
1928        """
1929       
1930        dict = self.Mesh2IODict()
1931        load_mesh.loadASCII.export_mesh_file(ofile,dict)
1932
1933    # FIXME(DSG-DSG):Break this into two functions.
1934    #One for the outline points.
1935    #One for the mesh points.
1936    # Note: this function is not in the gui
1937    def exportPointsFile(self,ofile):
1938        """
1939        export a points (.xya or .pts)  file, ofile.
1940       
1941        """
1942       
1943        mesh_dict = self.Mesh2IODict()
1944        point_dict = {}
1945        point_dict['attributelist'] = {} #this will need to be expanded..
1946                                         # if attributes are brought back in.
1947        point_dict['geo_reference'] = self.geo_reference
1948        if mesh_dict['vertices'] == []:
1949            point_dict['pointlist'] = mesh_dict['points']
1950        else:
1951            point_dict['pointlist'] = mesh_dict['vertices']
1952
1953        load_mesh.loadASCII.export_points_file(ofile,point_dict)
1954
1955
1956    def import_ungenerate_file(self,ofile, tag=None):
1957        """
1958        Imports an ungenerate file, from arcGIS into mesh.
1959
1960        ofile is the name of the ungenerated file.
1961        Tag is a string name to be taggged on each segment.
1962
1963        WARNING values are assumed to be absolute.
1964        geo-refs are not taken into account..
1965        """
1966   
1967        dict = importUngenerateFile(ofile)
1968        default_tag = Segment.get_default_tag()
1969        if tag is not None:
1970            Segment.set_default_tag(str(tag))
1971        self.addVertsSegs(dict)
1972        Segment.set_default_tag(default_tag)
1973
1974        # change the tag back to  it's default
1975   
1976       
1977########### IO CONVERTERS ##################
1978        """
1979        The dict fromat for IO with .tsh files is;
1980        (the triangulation)
1981        vertices: [[x1,y1],[x2,y2],...] (lists of doubles)
1982        vertex_attributes: [[a11,a12,...],[a21,a22],...] (lists of doubles)
1983        vertex_attribute_titles:[A1Title, A2Title ...] (A list of strings)
1984        segments: [[v1,v2],[v3,v4],...] (lists of integers)
1985        segment_tags : [tag,tag,...] list of strings
1986        triangles : [(v1,v2,v3), (v4,v5,v6),....] lists of points
1987        triangle tags: [s1,s2,...] A list of strings
1988        triangle neighbors: [[t1,t2,t3], [t4,t5,t6],..] lists of triangles
1989       
1990        (the outline)   
1991        points: [[x1,y1],[x2,y2],...] (lists of doubles)
1992        point_attributes: [[a11,a12,...],[a21,a22],...] (lists of doubles)
1993        outline_segments: [[point1,point2],[p3,p4],...] (lists of integers)
1994        outline_segment_tags : [tag1,tag2,...] list of strings
1995        holes : [[x1,y1],...](List of doubles, one inside each hole region)
1996        regions : [ [x1,y1],...] (List of 4 doubles)
1997        region_tags : [tag1,tag2,...] (list of strings)
1998        region_max_areas: [ma1,ma2,...] (A list of doubles)
1999        {Convension: A -ve max area means no max area}
2000       
2001        """
2002     
2003
2004                               
2005    def Mesh2IODict(self):
2006        """
2007        Convert the triangulation and outline info of a mesh to a dictionary
2008        structure
2009        """
2010        dict = self.Mesh2IOTriangulationDict()
2011        dict_mesh = self.Mesh2IOOutlineDict()
2012        for element in dict_mesh.keys():
2013            dict[element] = dict_mesh[element]
2014
2015        # add the geo reference
2016        dict['geo_reference'] = self.geo_reference
2017        return dict
2018   
2019    def Mesh2IOTriangulationDict(self):
2020        """
2021        Convert the Mesh to a dictionary of lists describing the
2022        triangulation variables;
2023       
2024        Used to produce .tsh file
2025        """
2026       
2027        meshDict = {}       
2028        vertices=[]
2029        vertex_attributes=[]
2030
2031        self.maxVertexIndex=0
2032        for vertex in self.meshVertices:
2033            vertex.index = self.maxVertexIndex
2034            vertices.append([vertex.x,vertex.y])
2035            vertex_attributes.append(vertex.attributes)           
2036            self.maxVertexIndex += 1
2037
2038        meshDict['vertices'] = vertices
2039        meshDict['vertex_attributes'] = vertex_attributes
2040        meshDict['vertex_attribute_titles'] = self.attributeTitles
2041        #segments
2042        segments=[]
2043        segment_tags=[]
2044        for seg in self.meshSegments:
2045            segments.append([seg.vertices[0].index,seg.vertices[1].index])
2046            segment_tags.append(seg.tag)
2047        meshDict['segments'] =segments
2048        meshDict['segment_tags'] =segment_tags
2049
2050        # Make sure that the indexation is correct
2051        index = 0
2052        for tri in self.meshTriangles:
2053            tri.index = index           
2054            index += 1
2055
2056        triangles = []
2057        triangle_tags = []
2058        triangle_neighbors = []
2059        for tri in self.meshTriangles: 
2060            triangles.append([tri.vertices[0].index,
2061                              tri.vertices[1].index,
2062                              tri.vertices[2].index]) 
2063            triangle_tags.append(tri.attribute)
2064            neighborlist = [-1,-1,-1]
2065            for neighbor,index in map(None,tri.neighbors,
2066                                      range(len(tri.neighbors))):
2067                if neighbor:   
2068                    neighborlist[index] = neighbor.index
2069            triangle_neighbors.append(neighborlist)
2070       
2071        meshDict['triangles'] = triangles
2072        meshDict['triangle_tags'] = triangle_tags
2073        meshDict['triangle_neighbors'] = triangle_neighbors
2074       
2075        #print "mesh.Mesh2IOTriangulationDict*)*)"
2076        #print meshDict
2077        #print "mesh.Mesh2IOTriangulationDict*)*)"
2078
2079        return meshDict
2080
2081                                                     
2082    def Mesh2IOOutlineDict(self, userVertices=None,
2083                        userSegments=None,
2084                        holes=None,
2085                        regions=None):
2086        """
2087        Convert the mesh outline to a dictionary of the lists needed for the
2088        triang module;
2089       
2090        Note, this adds an index attribute to the user Vertex objects.
2091
2092        Used to produce .tsh file and output to triangle
2093        """
2094        if userVertices is None:
2095            userVertices = self.getUserVertices()
2096        if userSegments is None:
2097            userSegments = self.getUserSegments()
2098        if holes is None:
2099            holes = self.getHoles()
2100        if regions is None:
2101            regions = self.getRegions()
2102           
2103        meshDict = {}
2104        #print "userVertices",userVertices
2105        #print "userSegments",userSegments
2106        pointlist=[]
2107        pointattributelist=[]
2108        index = 0
2109        for vertex in userVertices:
2110            vertex.index = index
2111            pointlist.append([vertex.x,vertex.y])
2112            pointattributelist.append(vertex.attributes)
2113           
2114            index += 1
2115        meshDict['points'] = pointlist
2116        meshDict['point_attributes'] = pointattributelist
2117
2118        segmentlist=[]
2119        segmenttaglist=[]
2120        for seg in userSegments:
2121            segmentlist.append([seg.vertices[0].index,seg.vertices[1].index])
2122            segmenttaglist.append(seg.tag)
2123        meshDict['outline_segments'] =segmentlist
2124        meshDict['outline_segment_tags'] =segmenttaglist
2125       
2126        holelist=[]
2127        for hole in holes:
2128            holelist.append([hole.x,hole.y]) 
2129        meshDict['holes'] = holelist
2130       
2131        regionlist=[]
2132        regiontaglist = []
2133        regionmaxarealist = []
2134        for region in regions:
2135            regionlist.append([region.x,region.y])
2136            regiontaglist.append(region.getTag())
2137           
2138            if (region.getMaxArea() != None): 
2139                regionmaxarealist.append(region.getMaxArea())
2140            else:
2141                regionmaxarealist.append( load_mesh.loadASCII.NOMAXAREA)
2142        meshDict['regions'] = regionlist
2143        meshDict['region_tags'] = regiontaglist
2144        meshDict['region_max_areas'] = regionmaxarealist
2145        #print "*(*("
2146        #print meshDict
2147        #print meshDict['regionlist']
2148        #print "*(*("
2149        return meshDict
2150
2151    def IOTriangulation2Mesh(self, genDict):
2152        """
2153        Set the mesh attributes given an tsh IO dictionary
2154        """
2155        #Clear the current generated mesh values
2156        self.meshTriangles=[]
2157        self.attributeTitles=[]
2158        self.meshSegments=[]
2159        self.meshVertices=[]
2160
2161        #print "mesh.setTriangulation@#@#@#"
2162        #print genDict
2163        #print "@#@#@#"
2164
2165        self.maxVertexIndex = 0
2166        for point in genDict['vertices']:
2167            v=Vertex(point[0], point[1])
2168            v.index =  self.maxVertexIndex
2169            self.maxVertexIndex +=1
2170            self.meshVertices.append(v)
2171
2172        self.attributeTitles = genDict['vertex_attribute_titles']
2173
2174        index = 0
2175        for seg,tag in map(None,genDict['segments'],genDict['segment_tags']):
2176            segObject = Segment( self.meshVertices[int(seg[0])],
2177                           self.meshVertices[int(seg[1])], tag = tag )
2178            segObject.index = index
2179            index +=1
2180            self.meshSegments.append(segObject)
2181
2182        index = 0
2183        for triangle in genDict['triangles']:
2184            tObject =Triangle( self.meshVertices[int(triangle[0])],
2185                        self.meshVertices[int(triangle[1])],
2186                        self.meshVertices[int(triangle[2])] )
2187            tObject.index = index
2188            index +=1
2189            self.meshTriangles.append(tObject)
2190
2191        index = 0
2192        for att in genDict['triangle_tags']:
2193            if att == []:
2194                self.meshTriangles[index].setAttribute("")
2195            else:
2196                self.meshTriangles[index].setAttribute(att)
2197            index += 1
2198           
2199        index = 0
2200        for att in genDict['vertex_attributes']:
2201            if att == None:
2202                self.meshVertices[index].setAttributes([])
2203            else:
2204                self.meshVertices[index].setAttributes(att)
2205            index += 1
2206   
2207        index = 0
2208        for triangle in genDict['triangle_neighbors']:
2209            # Build a list of triangle object neighbors
2210            ObjectNeighbor = []
2211            for neighbor in triangle:
2212                if ( neighbor != -1):
2213                    ObjectNeighbor.append(self.meshTriangles[int(neighbor)])
2214                else:
2215                    ObjectNeighbor.append(None)
2216            self.meshTriangles[index].setNeighbors(ObjectNeighbor[0],
2217                                                   ObjectNeighbor[1],
2218                                                   ObjectNeighbor[2])
2219            index += 1
2220
2221
2222    def IOOutline2Mesh(self, genDict):
2223        """
2224        Set the outline (user Mesh attributes) given a IO tsh dictionary
2225       
2226        mesh is an instance of a mesh object
2227        """
2228        #Clear the current user mesh values
2229        self.clearUserSegments()
2230        self.userVertices=[]
2231        self.Holes=[]
2232        self.Regions=[]
2233
2234        #print "mesh.IOOutline2Mesh@#@#@#"
2235        #print "genDict",genDict
2236        #print "@#@#@#"
2237       
2238        #index = 0
2239        for point in genDict['points']:
2240            v=Vertex(point[0], point[1])
2241            #v.index = index
2242            #index +=1
2243            self.userVertices.append(v)
2244
2245        #index = 0
2246        for seg,tag in map(None,genDict['outline_segments'],
2247                           genDict['outline_segment_tags']):
2248
2249            segObject = Segment( self.userVertices[int(seg[0])],
2250                           self.userVertices[int(seg[1])], tag = tag )
2251            #segObject.index = index
2252            #index +=1
2253            self.userSegments.append(segObject)
2254
2255# Remove the loading of attribute info.
2256# Have attribute info added using least_squares in pyvolution
2257#         index = 0
2258#         for att in genDict['point_attributes']:
2259#             if att == None:
2260#                 self.userVertices[index].setAttributes([])
2261#             else:
2262#                 self.userVertices[index].setAttributes(att)
2263#            index += 1
2264       
2265        #index = 0
2266        for point in genDict['holes']:
2267            h=Hole(point[0], point[1])
2268            #h.index = index
2269            #index +=1
2270            self.holes.append(h)
2271
2272        #index = 0
2273        for reg,att,maxArea in map(None,
2274                                   genDict['regions'],
2275                                   genDict['region_tags'],
2276                                   genDict['region_max_areas']):
2277            if maxArea > 0:  # maybe I should ref NOMAXAREA? Prob' not though
2278                Object = Region( reg[0],
2279                                 reg[1],
2280                                 tag = att,
2281                                 maxArea = maxArea)
2282            else:
2283                Object = Region( reg[0],
2284                                 reg[1],
2285                                 tag = att)
2286               
2287            #Object.index = index
2288            #index +=1
2289            self.regions.append(Object)
2290 
2291############################################
2292
2293       
2294    def refineSet(self,setName):
2295        Triangles = self.sets[self.setID[setName]]
2296        Refine(self,Triangles)
2297
2298    def selectAllTriangles(self):
2299        A=[]
2300        A.extend(self.meshTriangles)
2301        if not('All' in self.setID.keys()):
2302            self.setID['All']=len(self.sets)
2303            self.sets.append(A)
2304        else:
2305            self.sets[self.setID['All']]=A
2306        return 'All'
2307        # and objectIDs
2308
2309
2310    def clearSelection(self):
2311        A = []
2312        if not('None' in self.setID.keys()):
2313            self.setID['None']=len(self.sets)
2314            self.sets.append(A)
2315        return 'None'
2316
2317    def drawSet(self,canvas,setName,SCALE,colour=SET_COLOUR):
2318    #FIXME Draws over previous triangles - may bloat canvas
2319        Triangles = self.sets[self.setID[setName]]
2320        for triangle in Triangles:
2321            triangle.draw(canvas,1,
2322                          scale = SCALE,
2323                          colour = colour)
2324           
2325    def undrawSet(self,canvas,setName,SCALE,colour='green'):
2326    #FIXME Draws over previous lines - may bloat canvas
2327        Triangles = self.sets[self.setID[setName]]
2328        for triangle in Triangles:
2329            triangle.draw(canvas,1,
2330                          scale = SCALE,
2331                          colour = colour)
2332
2333    def weed(self,Vertices,Segments):
2334        #Depreciated
2335        #weed out existing duplicates
2336        print 'len(self.getUserSegments())'
2337        print len(self.getUserSegments())
2338        print 'len(self.getUserVertices())'
2339        print len(self.getUserVertices())
2340
2341        point_keys = {}
2342        for vertex in Vertices:
2343            point = (vertex.x,vertex.y)
2344            point_keys[point]=vertex
2345        #inlined would looks very ugly
2346
2347        line_keys = {}
2348        for segment in Segments:
2349            vertex1 = segment.vertices[0]
2350            vertex2 = segment.vertices[1]
2351            point1 = (vertex1.x,vertex1.y)
2352            point2 = (vertex2.x,vertex2.y)
2353            segment.vertices[0]=point_keys[point1]
2354            segment.vertices[1]=point_keys[point2]
2355            vertex1 = segment.vertices[0]
2356            vertex2 = segment.vertices[1]
2357            point1 = (vertex1.x,vertex1.y)
2358            point2 = (vertex2.x,vertex2.y)
2359            line1 = (point1,point2)
2360            line2 = (point2,point1)
2361            if not (line_keys.has_key(line1) \
2362                 or line_keys.has_key(line2)):
2363                 line_keys[line1]=segment
2364        Vertices=point_keys.values()
2365        Segments=line_keys.values()
2366        return Vertices,Segments
2367
2368    def segs_to_dict(self,segments):
2369        dict={}
2370        for segment in segments:
2371            vertex1 = segment.vertices[0]
2372            vertex2 = segment.vertices[1]
2373            point1 = (vertex1.x,vertex1.y)
2374            point2 = (vertex2.x,vertex2.y)
2375            line = (point1,point2)
2376            dict[line]=segment
2377        return dict
2378
2379    def seg2line(self,s):
2380        return ((s.vertices[0].x,s.vertices[0].y,)\
2381                (s.vertices[1].x,s.vertices[1].y))
2382
2383    def line2seg(self,line,tag=None):
2384        point0 = self.point2ver(line[0])
2385        point1 = self.point2ver(line[1])
2386        return Segment(point0,point1,tag=tag)
2387
2388    def ver2point(self,vertex):
2389        return (vertex.x,vertex.y)
2390
2391    def point2ver(self,point):
2392        return Vertex(point[0],point[1])
2393
2394    def smooth_polySet(self,min_radius=0.05):
2395        #for all pairs of connecting segments:
2396        #    propose a new segment that replaces the 2
2397
2398        #    If the difference between the new segment
2399        #    and the old lines is small: replace the
2400        #    old lines.
2401
2402        seg2line = self.seg2line
2403        ver2point= self.ver2point
2404        line2seg = self.line2seg
2405        point2ver= self.point2ver
2406
2407        #create dictionaries of lines -> segments
2408        userSegments = self.segs_to_dict(self.userSegments)
2409        alphaSegments = self.segs_to_dict(self.alphaUserSegments)
2410
2411        #lump user and alpha segments
2412        for key in alphaSegments.keys():
2413            userSegments[key]=alphaSegments[key]
2414
2415        #point_keys = tuple -> vertex
2416        #userVertices = vertex -> [line,line] - lines from that node
2417        point_keys = {}
2418        userVertices={}
2419        for vertex in self.getUserVertices():
2420            point = ver2point(vertex)
2421            if not point_keys.has_key(point):
2422                point_keys[point]=vertex
2423                userVertices[vertex]=[]
2424        for key in userSegments.keys():
2425            line = key
2426            point_0 = key[0]
2427            point_1 = key[1]
2428            userVertices[point_keys[point_0]].append(line)
2429            userVertices[point_keys[point_1]].append(line)
2430
2431        for point in point_keys.keys():
2432            try:
2433            #removed keys can cause keyerrors
2434                vertex = point_keys[point]
2435                lines = userVertices[vertex]
2436   
2437                #if there are 2 lines on the node
2438                if len(lines)==2:
2439                    line_0 = lines[0]
2440                    line_1 = lines[1]
2441   
2442                    #if the tags are the the same on the 2 lines
2443                    if userSegments[line_0].tag == userSegments[line_1].tag:
2444                        tag = userSegments[line_0].tag
2445   
2446                        #point_a is one of the next nodes, point_b is the other
2447                        if point==line_0[0]:
2448                            point_a = line_0[1]
2449                        if point==line_0[1]:
2450                            point_a = line_0[0]
2451                        if point==line_1[0]:
2452                            point_b = line_1[1]
2453                        if point==line_1[1]:
2454                            point_b = line_1[0]
2455   
2456   
2457                        #line_2 is proposed
2458                        line_2 = (point_a,point_b)
2459
2460                        #calculate the area of the triangle between
2461                        #the two existing segments and the proposed
2462                        #new segment
2463                        ax = point_a[0]
2464                        ay = point_a[1]
2465                        bx = point_b[0]
2466                        by = point_b[1]
2467                        cx = point[0]
2468                        cy = point[1]
2469                        area=abs((bx*ay-ax*by)+(cx*by-bx*cy)+(ax*cy-cx*ay))/2
2470
2471                        #calculate the perimeter
2472                        len_a =  ((cx-bx)**2+(cy-by)**2)**0.5 
2473                        len_b =  ((ax-cx)**2+(ay-cy)**2)**0.5 
2474                        len_c =  ((bx-ax)**2+(by-ay)**2)**0.5 
2475                        perimeter = len_a+len_b+len_c
2476
2477                        #calculate the radius
2478                        r = area/(2*perimeter)
2479
2480                        #if the radius is small: then replace the existing
2481                        #segments with the new one
2482                        if r < min_radius:
2483                            if len_c < min_radius: append = False
2484                            else: append = True
2485                            #if the new seg is also time, don't add it
2486                            if append:
2487                                segment = self.line2seg(line_2,tag=tag)
2488
2489                            list_a=userVertices[point_keys[point_a]]
2490                            list_b=userVertices[point_keys[point_b]]
2491
2492                            if line_0 in list_a:
2493                                list_a.remove(line_0)
2494                            else:
2495                                list_a.remove(line_1)
2496
2497                            if line_0 in list_b:
2498                                list_b.remove(line_0)
2499                            else:
2500                                list_b.remove(line_1)
2501
2502                            if append:
2503                                list_a.append(line_2)
2504                                list_b.append(line_2)
2505                            else:
2506                                if len(list_a)==0:
2507                                    userVertices.pop(point_keys[point_a])
2508                                    point_keys.pop(point_a)
2509                                if len(list_b)==0:
2510                                    userVertices.pop(point_keys[point_b])
2511                                    point_keys.pop(point_b)
2512
2513                            userVertices.pop(point_keys[point])
2514                            point_keys.pop(point)
2515                            userSegments.pop(line_0)
2516                            userSegments.pop(line_1)
2517
2518                            if append:
2519                                userSegments[line_2]=segment
2520            except:
2521                pass
2522
2523        #self.userVerticies = userVertices.keys()
2524        #self.userSegments = []
2525        #for key in userSegments.keys():
2526        #    self.userSegments.append(userSegments[key])
2527        #self.alphaUserSegments = []
2528
2529        self.userVerticies = []
2530        self.userSegments = []
2531        self.alphaUserSegments = []
2532
2533        return userVertices,userSegments,alphaSegments
2534
2535    def triangles_to_polySet(self,setName):
2536        #self.smooth_polySet()
2537
2538        seg2line = self.seg2line
2539        ver2point= self.ver2point
2540        line2seg = self.line2seg
2541        point2ver= self.point2ver
2542
2543        from Numeric import array,allclose
2544        #turn the triangles into a set
2545        Triangles = self.sets[self.setID[setName]]
2546        Triangles_dict = {}
2547        for triangle in Triangles:
2548            Triangles_dict[triangle]=None
2549 
2550
2551        #create a dict of points to vertexes (tuple -> object)
2552        #also create a set of vertexes (object -> True)
2553        point_keys = {}
2554        userVertices={}
2555        for vertex in self.getUserVertices():
2556            point = ver2point(vertex)
2557            if not point_keys.has_key(point):
2558                point_keys[point]=vertex
2559                userVertices[vertex]=True
2560
2561        #create a dict of lines to segments (tuple -> object)
2562        userSegments = self.segs_to_dict(self.userSegments)
2563        #append the userlines in an affine linespace
2564        affine_lines = Affine_Linespace()
2565        for line in userSegments.keys():
2566            affine_lines.append(line)
2567        alphaSegments = self.segs_to_dict(self.alphaUserSegments)
2568        for line in alphaSegments.keys():
2569            affine_lines.append(line)
2570       
2571        for triangle in Triangles:
2572            for i in (0,1,2):
2573                #for every triangles neighbour:
2574                if not Triangles_dict.has_key(triangle.neighbors[i]):
2575                #if the neighbour is not in the set:
2576                    a = triangle.vertices[i-1]
2577                    b = triangle.vertices[i-2]
2578                    #Get possible matches:
2579                    point_a = ver2point(a)
2580                    point_b = ver2point(b)
2581                    midpoint = ((a.x+b.x)/2,(a.y+b.y)/2)
2582                    line = (point_a,point_b)
2583                    tag = None
2584
2585
2586                    #this bit checks for matching lines
2587                    possible_lines = affine_lines[line] 
2588                    possible_lines = unique(possible_lines)
2589                    found = 0                           
2590                    for user_line in possible_lines:
2591                        if self.point_on_line(midpoint,user_line):
2592                            found+=1
2593                            assert found<2
2594                            if userSegments.has_key(user_line):
2595                                parent_segment = userSegments.pop(user_line)
2596                            if alphaSegments.has_key(user_line):
2597                                parent_segment = alphaSegments.pop(user_line)
2598                            tag = parent_segment.tag
2599                            offspring = [line]
2600                            offspring.extend(self.subtract_line(user_line,line))
2601                            affine_lines.remove(user_line)
2602                            for newline in offspring:
2603                                line_vertices = []
2604                                for point in newline:
2605                                    if point_keys.has_key(point):
2606                                        vert = point_keys[point]
2607                                    else:
2608                                        vert = Vertex(point[0],point[1])
2609                                        userVertices[vert]=True
2610                                        point_keys[point]=vert
2611                                    line_vertices.append(vert)
2612                                segment = Segment(line_vertices[0],line_vertices[1],tag)
2613                                userSegments[newline]=segment
2614                                affine_lines.append(newline)
2615                            #break
2616                    assert found<2
2617
2618
2619
2620                    #if no matching lines
2621                    if not found:
2622                        line_vertices = []
2623                        for point in line:
2624                            if point_keys.has_key(point):
2625                                vert = point_keys[point]
2626                            else:
2627                                vert = Vertex(point[0],point[1])
2628                                userVertices[vert]=True
2629                                point_keys[point]=vert
2630                            line_vertices.append(vert)
2631                        segment = Segment(line_vertices[0],line_vertices[1],tag)
2632                        userSegments[line]=segment
2633                        affine_lines.append(line)
2634       
2635        self.userVerticies = []
2636        self.userSegments = []
2637        self.alphaUserSegments = []
2638
2639        return userVertices,userSegments,alphaSegments
2640
2641    def subtract_line(self,parent,child):
2642        #Subtracts child from parent
2643        #Requires that the child is a
2644        #subline of parent to work.
2645
2646        from Numeric import allclose,dot,array
2647        A= parent[0]
2648        B= parent[1]
2649        a = child[0]
2650        b = child[1]
2651
2652        A_array = array(parent[0])
2653        B_array = array(parent[1])
2654        a_array   = array(child[0])
2655        b_array   = array(child[1])
2656
2657        assert not A == B
2658        assert not a == b
2659
2660        answer = []
2661
2662        #if the new line does not share a
2663        #vertex with the old one
2664        if not (allclose(A_array,a_array)\
2665             or allclose(B_array,b_array)\
2666             or allclose(A_array,b_array)\
2667             or allclose(a_array,B_array)):
2668            if dot(A_array-a_array,A_array-a_array) \
2669            < dot(A_array-b_array,A_array-b_array):
2670                sibling1 = (A,a)
2671                sibling2 = (B,b)
2672                return [sibling1,sibling2]
2673            else:
2674                sibling1 = (A,b)
2675                sibling2 = (B,a)
2676                return [sibling1,sibling2]
2677
2678        elif allclose(A_array,a_array):
2679            if allclose(B_array,b_array):
2680                return []
2681            else:
2682                sibling = (b,B)
2683                return [sibling]
2684        elif allclose(B_array,b_array):
2685            sibling = (a,A)
2686            return [sibling]
2687
2688        elif allclose(A_array,b_array):
2689            if allclose(B,a):
2690                return []
2691            else:
2692                sibling = (a,B)
2693                return [sibling]
2694        elif allclose(a_array,B_array):
2695            sibling = (b,A)
2696            return [sibling]
2697
2698    def point_on_line(self,point,line):
2699        #returns true within a tolerance of 3 degrees
2700        x=point[0]
2701        y=point[1]
2702        x0=line[0][0]
2703        x1=line[1][0]
2704        y0=line[0][1]
2705        y1=line[1][1]
2706        from Numeric import array, dot, allclose
2707        from math import sqrt
2708        tol = 3. #DEGREES
2709        tol = tol*3.1415/180
2710
2711        a = array([x - x0, y - y0]) 
2712        a_normal = array([a[1], -a[0]])     
2713        len_a_normal = sqrt(sum(a_normal**2)) 
2714
2715        b = array([x1 - x0, y1 - y0])         
2716        len_b = sqrt(sum(b**2)) 
2717   
2718        if abs(dot(a_normal, b)/(len_b*len_a_normal))< tol:
2719            #Point is somewhere on the infinite extension of the line
2720
2721            len_a = sqrt(sum(a**2))                                     
2722            if dot(a, b) >= 0 and len_a <= len_b:
2723               return True
2724            else:   
2725               return False
2726        else:
2727          return False
2728
2729    def line_length(self,line):
2730        x0=line[0][0]
2731        x1=line[1][0]
2732        y0=line[0][1]
2733        y1=line[1][1]
2734        return ((x1-x0)**2-(y1-y0)**2)**0.5     
2735
2736    def threshold(self,setName,min=None,max=None,attribute_name = 'elevation'):
2737        """
2738        threshold using  d
2739        """
2740        triangles = self.sets[self.setID[setName]]
2741        A = []
2742
2743        if attribute_name in self.attributeTitles:
2744            i = self.attributeTitles.index(attribute_name)
2745        else: i = -1#no attribute
2746        if not max == None:
2747            for t in triangles:
2748                if (min<self.av_att(t,i)<max):
2749                    A.append(t)
2750        else:
2751            for t in triangles:
2752                if (min<self.av_att(t,i)):
2753                    A.append(t)
2754        self.sets[self.setID[setName]] = A
2755
2756    def general_threshold(self,setName,min=None,max=None\
2757              ,attribute_name = 'elevation',function=None):
2758        """
2759        Thresholds the triangles
2760        """
2761        from visual.graph import arange,ghistogram,color as colour
2762        triangles = self.sets[self.setID[setName]]
2763        A = []
2764        data=[]
2765        #data is for the graph
2766
2767        if attribute_name in self.attributeTitles:
2768            i = self.attributeTitles.index(attribute_name)
2769        else: i = -1
2770        if not max == None:
2771            for t in triangles:
2772                value=function(t,i)
2773                if (min<value<max):
2774                    A.append(t)
2775                data.append(value)
2776        else:
2777            for t in triangles:
2778                value=function(t,i)
2779                if (min<value):
2780                    A.append(t)
2781                data.append(value)
2782        self.sets[self.setID[setName]] = A
2783
2784        if self.visualise_graph:
2785            if len(data)>0:
2786                max=data[0]
2787                min=data[0]
2788                for value in data:
2789                    if value > max:
2790                        max = value
2791                    if value < min:
2792                        min = value
2793
2794                inc = (max-min)/100
2795
2796                histogram = ghistogram(bins=arange(min,max,inc),\
2797                             color = colour.red)
2798                histogram.plot(data=data)
2799       
2800    def av_att(self,triangle,i):
2801        if i==-1: return 1
2802        else:
2803            #evaluates the average attribute of the vertices of a triangle.
2804            V = triangle.getVertices()
2805            a0 = (V[0].attributes[i])
2806            a1 = (V[1].attributes[i])
2807            a2 = (V[2].attributes[i])
2808            return (a0+a1+a2)/3
2809
2810    def Courant_ratio(self,triangle,index):
2811        """
2812        Uses the courant threshold
2813        """
2814        e = self.av_att(triangle,index)
2815        A = triangle.calcArea()
2816        P = triangle.calcP()
2817        r = A/(2*P)
2818        e = max(0.1,abs(e))
2819        return r/e**0.5
2820
2821    def Gradient(self,triangle,index):
2822        V = triangle.vertices
2823        x0, y0, x1, y1, x2, y2, q0, q1, q2 = V[0].x,V[0].y,V[1].x,V[1].y,V[2].x,V[2].y,V[0].attributes[index],V[1].attributes[index],V[2].attributes[index]
2824        grad_x,grad_y = gradient(x0, y0, x1, y1, x2, y2, q0, q1, q2)
2825        if ((grad_x**2)+(grad_y**2))**(0.5)<0:
2826            print ((grad_x**2)+(grad_y**2))**(0.5)
2827        return ((grad_x**2)+(grad_y**2))**(0.5)
2828   
2829
2830    def append_triangle(self,triangle):
2831        self.meshTriangles.append(triangle)
2832
2833    def replace_triangle(self,triangle,replacement):
2834        i = self.meshTriangles.index(triangle)
2835        self.meshTriangles[i]=replacement
2836        assert replacement in self.meshTriangles
2837
2838def importUngenerateFile(ofile):
2839    """
2840    import a file, ofile, with the format
2841    [poly]
2842    poly format:
2843    First line:  <# of vertices> <x centroid> <y centroid>
2844    Following lines: <x> <y>
2845    last line:  "END"
2846
2847    Note: These are clockwise.
2848    """
2849    fd = open(ofile,'r')
2850    Dict = readUngenerateFile(fd)
2851    fd.close()
2852    return Dict
2853
2854def readUngenerateFile(fd):
2855    """
2856    import a file, ofile, with the format
2857    [poly]
2858    poly format:
2859    First line:  <# of polynomial> <x centroid> <y centroid>
2860    Following lines: <x> <y>
2861    last line:  "END"
2862    """
2863    END_DELIMITER = 'END'
2864   
2865    points = []
2866    segments = []
2867   
2868    isEnd = False
2869    line = fd.readline() #not used <# of polynomial> <x> <y>
2870    while not isEnd:
2871        line = fd.readline()
2872        fragments = line.split()
2873        vert = [float(fragments.pop(0)),float(fragments.pop(0))]
2874        points.append(vert)
2875        PreviousVertIndex = len(points)-1
2876        firstVertIndex = PreviousVertIndex
2877       
2878        line = fd.readline() #Read the next line
2879        while not line.startswith(END_DELIMITER): 
2880            #print "line >" + line + "<"
2881            fragments = line.split()
2882            vert = [float(fragments.pop(0)),float(fragments.pop(0))]
2883            points.append(vert)
2884            thisVertIndex = len(points)-1
2885            segment = [PreviousVertIndex,thisVertIndex]
2886            segments.append(segment)
2887            PreviousVertIndex = thisVertIndex
2888            line = fd.readline() #Read the next line
2889            i =+ 1
2890        # If the last and first segments are the same,
2891        # Remove the last segment and the last vertex
2892        # then add a segment from the second last vert to the 1st vert
2893        thisVertIndex = len(points)-1
2894        firstVert = points[firstVertIndex]
2895        thisVert = points[thisVertIndex]
2896        #print "firstVert",firstVert
2897        #print "thisVert",thisVert
2898        if (firstVert[0] == thisVert[0] and firstVert[1] == thisVert[1]):
2899            points.pop()
2900            segments.pop()
2901            thisVertIndex = len(points)-1
2902            segments.append([thisVertIndex, firstVertIndex])
2903       
2904        line = fd.readline() # read <# of polynomial> <x> <y> OR END
2905        #print "line >>" + line + "<<"
2906        if line.startswith(END_DELIMITER):
2907            isEnd = True
2908   
2909    #print "points", points       
2910    #print "segments", segments
2911    ungenerated_dict = {}
2912    ungenerated_dict['points'] = points
2913    ungenerated_dict['segments'] = segments
2914    return ungenerated_dict
2915
2916def importMeshFromFile(ofile):
2917    """returns a mesh object, made from a .xya/.pts or .tsh/.msh file
2918    Often raises IOError,RuntimeError
2919    """
2920    newmesh = None
2921    if (ofile[-4:]== ".xya" or ofile[-4:]== ".pts"):
2922        dict = load_mesh.loadASCII.import_points_file(ofile)
2923        dict['points'] = dict['pointlist']
2924        dict['outline_segments'] = []
2925        dict['outline_segment_tags'] = []
2926        dict['regions'] = []
2927        dict['region_tags'] = []
2928        dict['region_max_areas'] = []
2929        dict['holes'] = [] 
2930        newmesh= Mesh(geo_reference = dict['geo_reference'])
2931        newmesh.IOOutline2Mesh(dict)
2932        counter = newmesh.removeDuplicatedUserVertices()
2933        if (counter >0):
2934            print "%i duplicate vertices removed from dataset" % (counter)
2935    elif (ofile[-4:]== ".tsh" or ofile[-4:]== ".msh"):
2936        dict = load_mesh.loadASCII.import_mesh_file(ofile)
2937        #print "********"
2938        #print "zq mesh.dict",dict
2939        #print "********"
2940        newmesh= Mesh()
2941        newmesh.IOOutline2Mesh(dict)
2942        newmesh.IOTriangulation2Mesh(dict)
2943    else:
2944        raise RuntimeError
2945   
2946    if dict.has_key('geo_reference') and not dict['geo_reference'] == None:
2947        newmesh.geo_reference = dict['geo_reference']
2948    return newmesh
2949
2950def loadPickle(currentName):
2951    fd = open(currentName)
2952    mesh = pickle.load(fd)
2953    fd.close()
2954    return mesh
2955   
2956def square_outline(side_length = 1,up = "top", left = "left", right = "right",
2957                   down = "bottom", regions = False):
2958   
2959        a = Vertex (0,0)
2960        b = Vertex (0,side_length)
2961        c = Vertex (side_length,0)
2962        d = Vertex (side_length,side_length)
2963     
2964        s2 = Segment(b,d, tag = up)
2965        s3 = Segment(b,a, tag = left)
2966        s4 = Segment(d,c, tag = right)
2967        s5 = Segment(a,c, tag = down)
2968
2969        if regions:
2970            e =  Vertex (side_length/2,side_length/2)
2971            s6 = Segment(a,e, tag = down + left)
2972            s7 = Segment(b,e, tag = up + left)
2973            s8 = Segment(c,e, tag = down + right)
2974            s9 = Segment(d,e, tag = up + right)
2975            r1 = Region(side_length/2,3.*side_length/4, tag = up)
2976            r2 = Region(1.*side_length/4,side_length/2, tag = left)
2977            r3 = Region(3.*side_length/4,side_length/2, tag = right)
2978            r4 = Region(side_length/2,1.*side_length/4, tag = down)
2979            mesh = Mesh(userVertices=[a,b,c,d,e],
2980                        userSegments=[s2,s3,s4,s5,s6,s7,s8,s9],
2981                        regions = [r1,r2,r3,r4])
2982        else:
2983            mesh = Mesh(userVertices=[a,b,c,d],
2984                 userSegments=[s2,s3,s4,s5])
2985     
2986        return mesh
2987
2988   
2989
2990def region_strings2ints(region_list):
2991    """Given a list of (x_int,y_int,tag_string) lists it returns a list of
2992    (x_int,y_int,tag_int) and a list to convert the tag_int's back to
2993    the tag_strings 
2994    """
2995    # Make sure "" has an index of 0
2996    region_list.reverse()
2997    region_list.append((1.0,2.0,""))
2998    region_list.reverse()
2999    convertint2string = []
3000    for i in xrange(len(region_list)):
3001        convertint2string.append(region_list[i][2])
3002        if len(region_list[i]) == 4: # there's an area value
3003            region_list[i] = (region_list[i][0],
3004                              region_list[i][1],i,region_list[i][3])
3005        elif len(region_list[i]) == 3: # no area value
3006            region_list[i] = (region_list[i][0],region_list[i][1],i)
3007        else:
3008            print "The region list has a bad size"
3009            # raise an error ..
3010            raise Error
3011
3012    #remove "" from the region_list
3013    region_list.pop(0)
3014   
3015    return [region_list, convertint2string]
3016
3017def region_ints2strings(region_list,convertint2string):
3018    """Reverses the transformation of region_strings2ints
3019    """
3020    if region_list[0] != []:
3021        for i in xrange(len(region_list)):
3022            region_list[i] = [convertint2string[int(region_list[i][0])]]
3023    return region_list
3024
3025def segment_ints2strings(intlist, convertint2string):
3026    """Reverses the transformation of segment_strings2ints """
3027    stringlist = []
3028    for x in intlist:
3029        stringlist.append(convertint2string[x])
3030    return stringlist
3031
3032def segment_strings2ints(stringlist, preset):
3033    """Given a list of strings return a list of 0 to n ints which represent
3034    the strings and a converting list of the strings, indexed by 0 to n ints.
3035    Also, input an initial converting list of the strings
3036    Note, the converting list starts off with
3037    ["internal boundary", "external boundary", "internal boundary"]
3038    example input and output
3039    input ["a","b","a","c"],["c"]
3040    output [[2, 1, 2, 0], ['c', 'b', 'a']]
3041
3042    the first element in the converting list will be
3043    overwritten with "".
3044    ?This will become the third item in the converting list?
3045   
3046    # note, order the initial converting list is important,
3047     since the index = the int tag
3048    """
3049    nodups = unique(stringlist)
3050    # note, order is important, the index = the int tag
3051    #preset = ["internal boundary", "external boundary"]
3052    #Remove the preset tags from the list with no duplicates
3053    nodups = [x for x in nodups if x not in preset]
3054
3055    try:
3056        nodups.remove("") # this has to go to zero
3057    except ValueError:
3058        pass
3059   
3060    # Add the preset list at the beginning of no duplicates
3061    preset.reverse()
3062    nodups.extend(preset)
3063    nodups.reverse()
3064
3065       
3066    convertstring2int = {}
3067    convertint2string = []
3068    index = 0
3069    for x in nodups:
3070        convertstring2int[x] = index
3071        convertint2string.append(x)
3072        index += 1
3073    convertstring2int[""] = 0
3074
3075    intlist = []
3076    for x in stringlist:
3077        intlist.append(convertstring2int[x])
3078    return [intlist, convertint2string]
3079
3080
3081
3082
3083def unique(s):
3084    """Return a list of the elements in s, but without duplicates.
3085
3086    For example, unique([1,2,3,1,2,3]) is some permutation of [1,2,3],
3087    unique("abcabc") some permutation of ["a", "b", "c"], and
3088    unique(([1, 2], [2, 3], [1, 2])) some permutation of
3089    [[2, 3], [1, 2]].
3090
3091    For best speed, all sequence elements should be hashable.  Then
3092    unique() will usually work in linear time.
3093
3094    If not possible, the sequence elements should enjoy a total
3095    ordering, and if list(s).sort() doesn't raise TypeError it's
3096    assumed that they do enjoy a total ordering.  Then unique() will
3097    usually work in O(N*log2(N)) time.
3098
3099    If that's not possible either, the sequence elements must support
3100    equality-testing.  Then unique() will usually work in quadratic
3101    time.
3102    """
3103
3104    n = len(s)
3105    if n == 0:
3106        return []
3107
3108    # Try using a dict first, as that's the fastest and will usually
3109    # work.  If it doesn't work, it will usually fail quickly, so it
3110    # usually doesn't cost much to *try* it.  It requires that all the
3111    # sequence elements be hashable, and support equality comparison.
3112    u = {}
3113    try:
3114        for x in s:
3115            u[x] = 1
3116    except TypeError:
3117        del u  # move on to the next method
3118    else:
3119        return u.keys()
3120
3121    # We can't hash all the elements.  Second fastest is to sort,
3122    # which brings the equal elements together; then duplicates are
3123    # easy to weed out in a single pass.
3124    # NOTE:  Python's list.sort() was designed to be efficient in the
3125    # presence of many duplicate elements.  This isn't true of all
3126    # sort functions in all languages or libraries, so this approach
3127    # is more effective in Python than it may be elsewhere.
3128    try:
3129        t = list(s)
3130        t.sort()
3131    except TypeError:
3132        del t  # move on to the next method
3133    else:
3134        assert n > 0
3135        last = t[0]
3136        lasti = i = 1
3137        while i < n:
3138            if t[i] != last:
3139                t[lasti] = last = t[i]
3140                lasti += 1
3141            i += 1
3142        return t[:lasti]
3143
3144    # Brute force is all that's left.
3145    u = []
3146    for x in s:
3147        if x not in u:
3148            u.append(x)
3149    return u
3150
3151"""Refines triangles
3152
3153   Implements the #triangular bisection?# algorithm.
3154 
3155   
3156"""
3157
3158def Refine(mesh, triangles):
3159    """
3160    Given a general_mesh, and a triangle number, split
3161    that triangle in the mesh in half. Then to prevent
3162    vertices and edges from meeting, keep refining
3163    neighbouring triangles until the mesh is clean.
3164    """
3165    state = BisectionState(mesh)
3166    for triangle in triangles:
3167        if not state.refined_triangles.has_key(triangle):
3168            triangle.rotate_longest_side()
3169            state.start(triangle)
3170            Refine_mesh(mesh, state)
3171
3172def Refine_mesh(mesh, state):
3173    """
3174    """
3175    state.getState(mesh)
3176    refine_triangle(mesh,state)
3177    state.evolve()
3178    if not state.end:
3179        Refine_mesh(mesh,state)
3180
3181def refine_triangle(mesh,state):
3182    split(mesh,state.current_triangle,state.new_point)
3183    if state.case == 'one':
3184        state.r[3]=state.current_triangle#triangle 2
3185
3186        new_triangle_id = len(mesh.meshTriangles)-1
3187        new_triangle = mesh.meshTriangles[new_triangle_id]
3188
3189        split(mesh,new_triangle,state.old_point)
3190        state.r[2]=new_triangle#triangle 1.2
3191        state.r[4]=mesh.meshTriangles[len(mesh.meshTriangles)-1]#triangle 1.1
3192        r = state.r
3193        state.repairCaseOne()
3194
3195    if state.case == 'two':
3196        state.r[2]=mesh.meshTriangles[len(mesh.meshTriangles)-1]#triangle 1
3197
3198        new_triangle = state.current_triangle
3199
3200        split(mesh,new_triangle,state.old_point)
3201
3202        state.r[3]=mesh.meshTriangles[len(mesh.meshTriangles)-1]#triangle 2.1
3203        state.r[4]=new_triangle#triangle 2.2
3204        r = state.r
3205
3206        state.repairCaseTwo()
3207
3208    if state.case == 'vertex':
3209        state.r[2]=state.current_triangle#triangle 2
3210        state.r[3]=mesh.meshTriangles[len(mesh.meshTriangles)-1]#triangle 1
3211        r = state.r
3212        state.repairCaseVertex()
3213       
3214    if state.case == 'start':
3215        state.r[2]=mesh.meshTriangles[len(mesh.meshTriangles)-1]#triangle 1
3216        state.r[3]=state.current_triangle#triangle 2
3217
3218    if state.next_case == 'boundary':
3219        state.repairCaseBoundary()
3220
3221
3222def split(mesh, triangle, new_point):
3223        """
3224        Given a mesh, triangle_id and a new point,
3225        split the corrosponding triangle into two
3226        new triangles and update the mesh.
3227        """
3228
3229        new_triangle1 = Triangle(new_point,triangle.vertices[0],triangle.vertices[1],attribute = triangle.attribute, neighbors = None)
3230        new_triangle2 = Triangle(new_point,triangle.vertices[2],triangle.vertices[0],attribute = triangle.attribute, neighbors = None)
3231
3232        new_triangle1.setNeighbors(triangle.neighbors[2],None,new_triangle2)
3233        new_triangle2.setNeighbors(triangle.neighbors[1],new_triangle1,None)
3234
3235        mesh.meshTriangles.append(new_triangle1)
3236
3237        triangle.vertices = new_triangle2.vertices
3238        triangle.neighbors = new_triangle2.neighbors
3239
3240
3241class State:
3242
3243    def __init__(self):
3244        pass
3245
3246class BisectionState(State):
3247
3248
3249    def __init__(self,mesh):
3250        self.len = len(mesh.meshTriangles)
3251        self.refined_triangles = {}
3252        self.mesh = mesh
3253        self.current_triangle = None
3254        self.case = 'start'
3255        self.end = False
3256        self.r = [None,None,None,None,None]
3257
3258    def start(self, triangle):
3259        self.current_triangle = triangle
3260        self.case = 'start'
3261        self.end = False
3262        self.r = [None,None,None,None,None]
3263
3264    def getState(self,mesh):
3265        if not self.case == 'vertex':
3266            self.new_point=self.getNewVertex(mesh, self.current_triangle)
3267            #self.neighbour=self.getNeighbour(mesh, self.current_triangle)
3268            self.neighbour = self.current_triangle.neighbors[0]
3269            if not self.neighbour is None:
3270                self.neighbour.rotate_longest_side()
3271            self.next_case = self.get_next_case(mesh,self.neighbour,self.current_triangle)
3272        if self.case == 'vertex':
3273            self.new_point=self.old_point
3274
3275
3276    def evolve(self):
3277        if self.case == 'vertex':
3278            self.end = True
3279
3280        self.last_case = self.case
3281        self.case = self.next_case
3282
3283        self.old_point = self.new_point
3284        self.current_triangle = self.neighbour
3285
3286        if self.case == 'boundary':
3287            self.end = True
3288        self.refined_triangles[self.r[2]]=1
3289        self.refined_triangles[self.r[3]]=1
3290        if not self.r[4] is None:
3291            self.refined_triangles[self.r[4]]=1
3292        self.r[0]=self.r[2]
3293        self.r[1]=self.r[3]
3294
3295
3296    def getNewVertex(self,mesh,triangle):
3297        coordinate1 = triangle.vertices[1]
3298        coordinate2 = triangle.vertices[2]
3299        a = ([coordinate1.x*1.,coordinate1.y*1.])
3300        b = ([coordinate2.x*1.,coordinate2.y*1.])
3301        attributes = []
3302        for i in range(len(coordinate1.attributes)):
3303            att = (coordinate1.attributes[i]+coordinate2.attributes[i])/2
3304            attributes.append(att)
3305        new_coordinate = [((a[0]-b[0])/2+b[0]),((a[1]-b[1])/2+b[1])]
3306        newVertex = Vertex(new_coordinate[0],new_coordinate[1], attributes = attributes)
3307        mesh.maxVertexIndex+=1
3308        newVertex.index = mesh.maxVertexIndex
3309        mesh.meshVertices.append(newVertex)
3310        return newVertex
3311
3312    def get_next_case(self, mesh,neighbour,triangle):
3313        """
3314        Given the locations of two neighbouring triangles,
3315        examine the interior indices of their vertices (i.e.
3316        0,1 or 2) to determine what how the neighbour needs
3317        to be refined.
3318        """
3319        if (neighbour is None):
3320                next_case = 'boundary'
3321        else:
3322                if triangle.vertices[1].x==neighbour.vertices[2].x:
3323                    if triangle.vertices[1].y==neighbour.vertices[2].y:
3324                        next_case = 'vertex'
3325                if triangle.vertices[1].x==neighbour.vertices[0].x:
3326                    if triangle.vertices[1].y==neighbour.vertices[0].y:
3327                        next_case = 'two'
3328                if triangle.vertices[1].x==neighbour.vertices[1].x:
3329                    if triangle.vertices[1].y==neighbour.vertices[1].y:
3330                        next_case = 'one'
3331        return next_case
3332
3333
3334
3335    def repairCaseVertex(self):
3336
3337        r = self.r
3338
3339
3340        self.link(r[0],r[2])
3341        self.repair(r[0])
3342
3343        self.link(r[1],r[3])
3344        self.repair(r[1])
3345
3346        self.repair(r[2])
3347
3348        self.repair(r[3])
3349
3350
3351    def repairCaseOne(self):
3352        r = self.rkey
3353
3354
3355        self.link(r[0],r[2])
3356        self.repair(r[0])
3357
3358        self.link(r[1],r[4])
3359        self.repair(r[1])
3360
3361        self.repair(r[4])
3362
3363    def repairCaseTwo(self):
3364        r = self.r
3365
3366        self.link(r[0],r[4])
3367        self.repair(r[0])
3368
3369        self.link(r[1],r[3])
3370        self.repair(r[1])
3371
3372        self.repair(r[4])
3373
3374    def repairCaseBoundary(self):
3375        r = self.r
3376        self.repair(r[2])
3377        self.repair(r[3])
3378
3379
3380
3381    def repair(self,triangle):
3382        """
3383        Given a triangle that knows its neighbours, this will
3384        force the neighbours to comply.
3385
3386        However, it needs to compare the vertices of triangles
3387        for this implementation
3388
3389        But it doesn't work for invalid neighbour structures
3390        """
3391        n=triangle.neighbors
3392        for i in (0,1,2):
3393            if not n[i] is None:
3394                for j in (0,1,2):#to find which side of the list is broken
3395                    if not (n[i].vertices[j] in triangle.vertices):
3396                    #ie if j is the side of n that needs fixing
3397                        k = j
3398                n[i].neighbors[k]=triangle
3399
3400
3401
3402    def link(self,triangle1,triangle2):
3403        """
3404        make triangle1 neighbors point to t
3405                #count = 0riangle2
3406        """
3407        count = 0
3408        for i in (0,1,2):#to find which side of the list is broken
3409            if not (triangle1.vertices[i] in triangle2.vertices):
3410                j = i
3411                count+=1
3412        assert count == 1
3413        triangle1.neighbors[j]=triangle2
3414
3415class Discretised_Tuple_Set:
3416    """
3417    if a={(0.0):[(0.01),(0.02)],(0.2):[(0.17)]}
3418    a[(0.01)]=a[(0.0)]=[(0.01),(0.02)]
3419    a[(10000)]=[] #NOT KEYERROR
3420
3421    a.append[(0.01)]
3422    => {0.0:[(0.01),(0.02),(0.01)],0.2:[(0.17)]}
3423
3424    #NOT IMPLIMENTED
3425    a.remove[(0.01)]
3426    => {(0.0):[(0.02),(0.01)],0.2:[(0.17)]}
3427
3428    a.remove[(0.17)]
3429    => {(0.0):[(0.02),(0.01)],0.2:[]}
3430    #NOT IMPLIMENTED
3431    at a.precision = 2:
3432        a.round_up_rel[0.0]=
3433        a.round_flat[0.0]=
3434        a.round_down_rel[0.0]=
3435
3436        a.up((0.1,2.04))=
3437
3438    If t_rel = 0, nothing gets rounded into
3439    two bins. If t_rel = 0.5, everything does.
3440
3441    Ideally, precision can be set high, so that multiple
3442    entries are rarely in the same bin. And t_rel should
3443    be low (<0.1 for 1 dimension!,<(0.1/n) for small n!!)
3444    so that it is rare to put items in mutiple bins.
3445
3446    Ex bins per entry = product(a,b,c...,n)
3447    a = 1 or 2 s.t. Ex(a) = 1+2*t_rel
3448    b = 1 or 2 ...
3449
3450    BUT!!! to avoid missing the right bin:
3451    (-10)**(precision+1)*t_rel must be greater than the
3452    greatest possible variation that an identical element
3453    can display.
3454
3455
3456    Note that if tol = 0.5 (the max allowed) 0.6 will round to .7 and .5
3457    but not .6 - this looks wrong, but note that *everything* will round,
3458    so .6 wont be missed as everything close to it will check in .7 and .5.
3459    """
3460    def __init__(self,p_rel = 6,t_rel = 0.01):
3461        self.__p_rel__ = p_rel
3462        self.__t_rel__ = t_rel
3463
3464        self.__p_abs__ = p_rel+1
3465        self.__t_abs__ = t_rel
3466
3467        assert t_rel <= 0.5
3468        self.__items__ = {}
3469        from math import frexp
3470        self.frexp = frexp
3471        roundings = [self.round_up_rel,\
3472        self.round_down_rel,self.round_flat_rel,\
3473        self.round_down_abs,self.round_up_abs,\
3474        self.round_flat_abs]#
3475
3476        self.roundings = roundings
3477
3478    def __repr__(self):
3479        return '%s'%self.__items__
3480
3481    def rounded_keys(self,key):
3482        key = tuple(key)
3483        keys = [key]
3484        keys = self.__rounded_keys__(key)
3485        return (keys)
3486
3487    def __rounded_keys__(self,key):
3488        keys = []
3489        rounded_key=list(key)
3490        rounded_values=list(key)
3491
3492        roundings = list(self.roundings)
3493
3494        #initialise rounded_values
3495        round = roundings.pop(0)
3496        for i in range(len(rounded_values)):
3497            rounded_key[i]=round(key[i])
3498            rounded_values[i]={}
3499            rounded_values[i][rounded_key[i]]=None
3500        keys.append(tuple(rounded_key))
3501       
3502        for round in roundings:
3503            for i in range(len(rounded_key)):
3504                rounded_value=round(key[i])
3505                if not rounded_values[i].has_key(rounded_value):
3506                    #ie unless round_up_rel = round_down_rel
3507                    #so the keys stay unique
3508                    for j in range(len(keys)):
3509                        rounded_key = list(keys[j])
3510                        rounded_key[i]=rounded_value
3511                        keys.append(tuple(rounded_key))
3512        return keys
3513
3514    def append(self,item):
3515        keys = self.rounded_keys(item)
3516        for key in keys:
3517            if self.__items__.has_key(key):
3518                self.__items__[key].append(item)
3519            else:
3520                self.__items__[key]=[item]
3521
3522    def __getitem__(self,key):
3523        answer = []
3524        keys = self.rounded_keys(key)
3525        for key in keys:
3526            if self.__items__.has_key(key):
3527                answer.extend(self.__items__[key])
3528        #if len(answer)==0:
3529        #    raise KeyError#FIXME or return KeyError
3530        #                  #FIXME or just return []?
3531        else:
3532            return answer #FIXME or unique(answer)?
3533
3534    def __delete__(self,item):
3535        keys = self.rounded_keys(item)
3536        answer = False
3537        #if any of the possible keys contains
3538        #a list, return true
3539        for key in keys:       
3540            if self.__items__.has_key(key):
3541                if item in self.__items__[key]:
3542                    self.__items__[key].remove(item)
3543
3544    def remove(self,item):
3545        self.__delete__(item)
3546
3547    def __contains__(self,item):
3548
3549        keys = self.rounded_keys(item)
3550        answer = False
3551        #if any of the possible keys contains
3552        #a list, return true
3553        for key in keys:       
3554            if self.__items__.has_key(key):
3555                if item in self.__items__[key]:
3556                    answer = True
3557        return answer
3558
3559
3560    def has_item(self,item):
3561        return self.__contains__(item)
3562
3563    def round_up_rel2(self,value):
3564         t_rel=self.__t_rel__
3565         #Rounding up the value
3566         m,e = self.frexp(value)
3567         m = m/2
3568         e = e + 1
3569         #m is the mantissa, e the exponent
3570         # 0.5 < |m| < 1.0
3571         m = m+t_rel*(10**-(self.__p_rel__))
3572         #bump m up
3573         m = round(m,self.__p_rel__)
3574         return m*(2.**e)
3575
3576    def round_down_rel2(self,value):
3577         t_rel=self.__t_rel__
3578         #Rounding down the value
3579         m,e = self.frexp(value)
3580         m = m/2
3581         e = e + 1
3582         #m is the mantissa, e the exponent
3583         # 0.5 < m < 1.0
3584         m = m-t_rel*(10**-(self.__p_rel__))
3585         #bump the |m| down, by 5% or whatever
3586         #self.p_rel dictates
3587         m = round(m,self.__p_rel__)
3588         return m*(2.**e)
3589
3590    def round_flat_rel2(self,value):
3591    #redundant
3592         m,e = self.frexp(value)
3593         m = m/2
3594         e = e + 1
3595         m = round(m,self.__p_rel__)
3596         return m*(2.**e)
3597
3598    def round_up_rel(self,value):
3599         t_rel=self.__t_rel__
3600         #Rounding up the value
3601         m,e = self.frexp(value)
3602         #m is the mantissa, e the exponent
3603         # 0.5 < |m| < 1.0
3604         m = m+t_rel*(10**-(self.__p_rel__))
3605         #bump m up
3606         m = round(m,self.__p_rel__)
3607         return m*(2.**e)
3608
3609    def round_down_rel(self,value):
3610         t_rel=self.__t_rel__
3611         #Rounding down the value
3612         m,e = self.frexp(value)
3613         #m is the mantissa, e the exponent
3614         # 0.5 < m < 1.0
3615         m = m-t_rel*(10**-(self.__p_rel__))
3616         #bump the |m| down, by 5% or whatever
3617         #self.p_rel dictates
3618         m = round(m,self.__p_rel__)
3619         return m*(2.**e)
3620
3621    def round_flat_rel(self,value):
3622    #redundant
3623         m,e = self.frexp(value)
3624         m = round(m,self.__p_rel__)
3625         return m*(2.**e)
3626
3627    def round_up_abs(self,value):
3628         t_abs=self.__t_abs__
3629         #Rounding up the value
3630         m = value+t_abs*(10**-(self.__p_abs__))
3631         #bump m up
3632         m = round(m,self.__p_abs__)
3633         return m
3634
3635    def round_down_abs(self,value):
3636         t_abs=self.__t_abs__
3637         #Rounding down the value
3638         m = value-t_abs*(10**-(self.__p_abs__))
3639         #bump the |m| down, by 5% or whatever
3640         #self.p_rel dictates
3641         m = round(m,self.__p_abs__)
3642         return m
3643
3644    def round_flat_abs(self,value):
3645    #redundant?
3646         m = round(value,self.__p_abs__)
3647         return m
3648
3649    def keys(self):
3650        return self.__items__.keys()
3651
3652
3653class Mapped_Discretised_Tuple_Set(Discretised_Tuple_Set):
3654    """
3655    This is a discretised tuple set, but
3656    based on a mapping. The mapping MUST
3657    return a sequence.
3658
3659    example:
3660    def weight(animal):
3661        return [animal.weight]
3662   
3663    a = Mapped_Discretised_Tuple_Set(weight)
3664    a.append[cow]
3665    a.append[fox]
3666    a.append[horse]
3667
3668    a[horse]    -> [cow,horse]
3669    a[dog]      -> [fox]
3670    a[elephant] -> []
3671    """
3672    def __init__(self,mapping,p_rel = 6, t_rel=0.01):
3673        Discretised_Tuple_Set.__init__\
3674         (self,p_rel,t_rel = t_rel)
3675        self.mapping = mapping
3676
3677    def rounded_keys(self,key):
3678        mapped_key = tuple(self.mapping(key))
3679        keys = self.__rounded_keys__(mapped_key)
3680        return keys
3681
3682class Affine_Linespace(Mapped_Discretised_Tuple_Set):
3683    """
3684    The affine linespace creates a way to record and compare lines.
3685    Precision is a bit of a hack, but it creates a way to avoid
3686    misses caused by round offs (between two lines of different
3687    lenghts, the short one gets rounded off more).
3688    I am starting to think that a quadratic search would be faster.
3689    Nearly.
3690    """
3691    def __init__(self,p_rel=4,t_rel=0.2):
3692        Mapped_Discretised_Tuple_Set.__init__\
3693            (self,self.affine_line,\
3694            p_rel=p_rel,t_rel=t_rel)
3695
3696        roundings = \
3697        [self.round_down_rel,self.round_up_rel,self.round_flat_rel]
3698        self.roundings = roundings
3699        #roundings = \
3700        #[self.round_down_abs,self.round_up_abs,self.round_flat_abs]
3701        #self.roundings = roundings
3702
3703    def affine_line(self,line):
3704        point_1 = line[0]
3705        point_2 = line[1]
3706        #returns the equation of a line
3707        #between two points, in the from
3708        #(a,b,-c), as in ax+by-c=0
3709        #or line *dot* (x,y,1) = (0,0,0)
3710
3711        #Note that it normalises the line
3712        #(a,b,-c) so Line*Line = 1.
3713        #This does not change the mathematical
3714        #properties, but it makes comparism
3715        #easier.
3716
3717        #There are probably better algorithms.
3718        x1 = point_1[0]
3719        x2 = point_2[0]
3720        y1 = point_1[1]
3721        y2 = point_2[1]
3722        dif_x = x1-x2
3723        dif_y = y1-y2
3724
3725        if dif_x == dif_y == 0:
3726            msg = 'points are the same'
3727            raise msg
3728        elif abs(dif_x)>=abs(dif_y):
3729            alpha = (-dif_y)/dif_x
3730            #a = alpha * b
3731            b = -1.
3732            c = (x1*alpha+x2*alpha+y1+y2)/2.
3733            a = alpha*b
3734        else:
3735            beta = dif_x/(-dif_y)
3736            #b = beta * a
3737            a = 1.
3738            c = (x1+x2+y1*beta+y2*beta)/2.
3739            b = beta*a
3740        mag = abs(a)+abs(b)
3741        #This does not change the mathematical
3742        #properties, but it makes comparism possible.
3743
3744        #note that the gradient is b/a, or (a/b)**-1.
3745        #so
3746
3747        #if a == 0:
3748        #    sign_a = 1.
3749        #else:
3750        #    sign_a = a/((a**2)**0.5)
3751        #if b == 0:
3752        #    sign_b = 1.
3753        #else:
3754        #    sign_b = b/((b**2)**0.5)
3755        #if c == 0:
3756        #    sign_c = 1.
3757        #else:
3758        #    sign_c = c/((c**2)**0.5)
3759        #a = a/mag*sign_a
3760        #b = b/mag*sign_b
3761        #c = c/mag*sign_c
3762        a = a/mag
3763        b = b/mag
3764        c = c/mag
3765        return a,b,c
3766
3767
3768# FIXME (DSG-DSG)
3769# To do-
3770# Create a clear interface. eg
3771# have the interface methods more at the top of this file and add comments
3772# for the interface functions/methods, use function_name (not functionName),
3773
3774#Currently
3775#function_name methods assume absolute values.  Geo-refs can be passed in.
3776#
3777
3778# instead of functionName
3779if __name__ == "__main__":
3780    #from mesh import *
3781    # THIS CAN BE DELETED
3782    m = Mesh()
3783    dict = importUngenerateFile("ungen_test.txt")
3784    m.addVertsSegs(dict)
3785    print m3
Note: See TracBrowser for help on using the repository browser.