source: inundation/pmesh/mesh.py @ 2778

Last change on this file since 2778 was 2778, checked in by steve, 18 years ago

On our 64 bit machine, ran into problem in pmesh/mesh.py where seg[0], seg[1], triangle,
and neighbor where not seen as integers. Had to wrap them in int(..)

File size: 124.0 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 triangle.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        """
1311        #Clear the current generated mesh values
1312        self.meshTriangles=[]
1313        self.attributeTitles=[]
1314        self.meshSegments=[]
1315        self.meshVertices=[]
1316
1317        #print "mesh.setTriangulation@#@#@#"
1318        #print genDict
1319        #print "@#@#@#"
1320
1321        self.maxVertexIndex = 0
1322        for point in genDict['generatedpointlist']:
1323            v=Vertex(point[0], point[1])
1324            v.index =  self.maxVertexIndex
1325            self.maxVertexIndex +=1
1326            self.meshVertices.append(v)
1327
1328        self.attributeTitles = genDict['generatedpointattributetitlelist']
1329
1330        index = 0
1331        for seg,marker in map(None,genDict['generatedsegmentlist'],
1332                              genDict['generatedsegmentmarkerlist']):
1333            segObject = Segment( self.meshVertices[int(seg[0])],
1334                           self.meshVertices[int(seg[1])], tag = marker )
1335            segObject.index = index
1336            index +=1
1337            self.meshSegments.append(segObject)
1338
1339        index = 0
1340        for triangle in genDict['generatedtrianglelist']:
1341            tObject =Triangle( self.meshVertices[int(triangle[0])],
1342                        self.meshVertices[int(triangle[1])],
1343                        self.meshVertices[int(triangle[2])] )
1344            tObject.index = index
1345            index +=1
1346            self.meshTriangles.append(tObject)
1347
1348        index = 0
1349        for att in genDict['generatedtriangleattributelist']:
1350            if att == []:
1351                self.meshTriangles[index].setAttribute("")
1352            else:
1353                self.meshTriangles[index].setAttribute(att[0])
1354            index += 1
1355           
1356        index = 0
1357        for att in genDict['generatedpointattributelist']:
1358            if att == None:
1359                self.meshVertices[index].setAttributes([])
1360            else:
1361                self.meshVertices[index].setAttributes(att)
1362            index += 1
1363   
1364        index = 0
1365        for triangle in genDict['generatedtriangleneighborlist']:
1366            # Build a list of triangle object neighbors
1367            ObjectNeighbor = []
1368            for neighbor in triangle:
1369                if ( neighbor != -1):
1370                    ObjectNeighbor.append(self.meshTriangles[int(neighbor)])
1371                else:
1372                    ObjectNeighbor.append(None)
1373            self.meshTriangles[index].setNeighbors(ObjectNeighbor[0],
1374                                                   ObjectNeighbor[1],
1375                                                   ObjectNeighbor[2])
1376            index += 1
1377
1378
1379    def setMesh(self, genDict):
1380        """
1381        Set the user Mesh attributes given a dictionary of the lists
1382        point list: [(x1,y1),(x2,y2),...] (Tuples of doubles) 
1383        point attribute list:[(P1att1,P1attt2, ...),(P2att1,P2attt2,...),...]
1384        segment list: [(point1,point2),(p3,p4),...] (Tuples of integers)
1385        segment tag list: [S1Tag, S2Tag, ...] (list of ints)
1386        region list: [(x1,y1),(x2,y2),...] (Tuples of doubles)
1387        region attribute list: ["","reservoir",""] list of strings
1388        region max area list:[real, None, Real,...] list of None and reals
1389       
1390        mesh is an instance of a mesh object
1391        """
1392        #Clear the current user mesh values
1393        self.clearUserSegments()
1394        self.userVertices=[]
1395        self.Holes=[]
1396        self.Regions=[]
1397
1398        #print "mesh.setMesh@#@#@#"
1399        #print genDict
1400        #print "@#@#@#"
1401       
1402        #index = 0
1403        for point in genDict['pointlist']:
1404            v=Vertex(point[0], point[1])
1405            #v.index = index
1406            #index +=1
1407            self.userVertices.append(v)
1408
1409        #index = 0
1410        for seg,tag in map(None,genDict['segmentlist'],
1411                           genDict['segmenttaglist']):
1412            segObject = Segment( self.userVertices[int(seg[0])],
1413                           self.userVertices[int(seg[1])], tag = tag )
1414            #segObject.index = index
1415            #index +=1
1416            self.userSegments.append(segObject)
1417
1418# Remove the loading of attribute info.
1419# Have attribute info added using least_squares in pyvolution
1420#         index = 0
1421#         for att in genDict['pointattributelist']:
1422#             if att == None:
1423#                 self.userVertices[index].setAttributes([])
1424#             else:
1425#                 self.userVertices[index].setAttributes(att)
1426#            index += 1
1427       
1428        #index = 0
1429        for point in genDict['holelist']:
1430            h=Hole(point[0], point[1])
1431            #h.index = index
1432            #index +=1
1433            self.holes.append(h)
1434
1435        #index = 0
1436        for reg,att,maxArea in map(None,
1437                                   genDict['regionlist'],
1438                                   genDict['regionattributelist'],
1439                                   genDict['regionmaxarealist']):
1440            Object = Region( reg[0],
1441                             reg[1],
1442                             tag = att,
1443                             maxArea = maxArea)
1444            #Object.index = index
1445            #index +=1
1446            self.regions.append(Object)
1447           
1448    def TestautoSegment(self):
1449        newsegs = []
1450        s1 = Segment(self.userVertices[0],
1451                               self.userVertices[1])
1452        s2 = Segment(self.userVertices[0],
1453                               self.userVertices[2])
1454        s3 = Segment(self.userVertices[2],
1455                               self.userVertices[1])
1456        if self.isUserSegmentNew(s1.vertices[0],s1.vertices[1]):
1457            newsegs.append(s1)
1458        if self.isUserSegmentNew(s2.vertices[0],s2.vertices[1]):
1459            newsegs.append(s2)
1460        if self.isUserSegmentNew(s3.vertices[0],s3.vertices[1]):
1461            newsegs.append(s3)
1462        #DSG!!!
1463        self.userSegments.extend(newsegs)
1464        return newsegs
1465
1466   
1467    def savePickle(self, currentName):
1468        fd = open(currentName, 'w')
1469        pickle.dump(self,fd)
1470        fd.close()
1471
1472    def autoSegmentHull(self):
1473        """
1474        initially work by running an executable
1475        Later compile the c code with a python wrapper.
1476
1477        Precon: There must be 3 or more vertices in the userVertices structure
1478        """
1479        newsegs = []
1480        inputfile = 'hull_in.txt'
1481        outputfile = inputfile + '-alf'
1482        #write vertices to file
1483        fd = open(inputfile,'w')
1484        for v in self.userVertices:
1485            fd.write(str(v.x))
1486            fd.write(' ')
1487            fd.write(str(v.y))
1488            fd.write('\n')
1489        fd.close()
1490       
1491        #run hull executable
1492        #warning need to compile hull for the current operating system
1493        command = 'hull.exe -A -i ' + inputfile
1494        os.system(command)
1495       
1496        #read results into this object
1497        fd = open(outputfile)
1498        lines = fd.readlines()
1499        fd.close()
1500        #print "(*(*(*("
1501        #print lines
1502        #print "(*(*(*("
1503        lines.pop(0) #remove the first (title) line
1504        for line in lines:
1505            vertindexs = line.split()
1506            #print 'int(vertindexs[0])', int(vertindexs[0])
1507            #print 'int(vertindexs[1])', int(vertindexs[1])
1508            #print 'self.userVertices[int(vertindexs[0])]' ,self.userVertices[int(vertindexs[0])]
1509            #print 'self.userVertices[int(vertindexs[1])]' ,self.userVertices[int(vertindexs[1])]
1510            v1 = self.userVertices[int(vertindexs[0])]
1511            v2 = self.userVertices[int(vertindexs[1])]
1512           
1513            if self.isUserSegmentNew(v1,v2):
1514                newseg = Segment(v1, v2)
1515                newsegs.append(newseg)
1516            #DSG!!!
1517        self.userSegments.extend(newsegs)
1518        return newsegs
1519    def autoSegmentFilter(self,raw_boundary=True,
1520                    remove_holes=False,
1521                    smooth_indents=False,
1522                    expand_pinch=False):
1523        """
1524        Precon: There is a self.shape
1525        """
1526        #FIXME remove the precon.  Internally check
1527        return self._boundary2mesh(raw_boundary=raw_boundary,
1528                    remove_holes=remove_holes,
1529                    smooth_indents=smooth_indents,
1530                    expand_pinch=expand_pinch)
1531   
1532       
1533   
1534    def autoSegment(self, alpha = None,
1535                    raw_boundary=True,
1536                    remove_holes=False,
1537                    smooth_indents=False,
1538                    expand_pinch=False): 
1539        """
1540        Precon: There must be 3 or more vertices in the userVertices structure
1541        """
1542        self._createBoundary(alpha=alpha)
1543        return self._boundary2mesh(raw_boundary=raw_boundary,
1544                    remove_holes=remove_holes,
1545                    smooth_indents=smooth_indents,
1546                    expand_pinch=expand_pinch)
1547
1548    def _createBoundary(self,alpha=None):
1549        """
1550        """
1551        points=[]
1552        for vertex in self.getUserVertices():
1553            points.append((vertex.x,vertex.y))
1554        self.shape = alpha_shape.alpha_shape.Alpha_Shape(points, alpha = alpha)
1555
1556
1557    def _boundary2mesh(self, raw_boundary=True,
1558                    remove_holes=False,
1559                    smooth_indents=False,
1560                    expand_pinch=False):
1561        """
1562        Precon there must be a shape object.
1563        """
1564        self.shape.set_boundary_type(raw_boundary=raw_boundary,
1565                                 remove_holes=remove_holes,
1566                                 smooth_indents=smooth_indents,
1567                                 expand_pinch=expand_pinch)
1568        boundary_segs = self.shape.get_boundary()
1569        #print "boundary_segs",boundary_segs
1570        segs2delete = self.alphaUserSegments
1571        #FIXME(DSG-DSG) this algorithm needs comments
1572        new_segs = {}
1573        #alpha_segs = []
1574        #user_segs = []
1575        for seg in boundary_segs:
1576            v1 = self.userVertices[int(seg[0])]
1577            v2 = self.userVertices[int(seg[1])]
1578            boundary_seg = Segment(v1, v2)
1579            new_segs[(v1,v2)] = boundary_seg
1580
1581        for user_seg in self.userSegments:
1582            if new_segs.has_key((user_seg.vertices[0],
1583                                user_seg.vertices[1])):
1584                del new_segs[user_seg.vertices[0],
1585                                user_seg.vertices[1]]
1586            elif new_segs.has_key((user_seg.vertices[1],
1587                                user_seg.vertices[0])):
1588                del new_segs[user_seg.vertices[1],
1589                                user_seg.vertices[0]]
1590               
1591        optimum_alpha = self.shape.get_alpha()
1592        alpha_segs_no_user_segs  = new_segs.values()
1593        self.alphaUserSegments = alpha_segs_no_user_segs
1594        return alpha_segs_no_user_segs, segs2delete, optimum_alpha
1595   
1596    def _boundary2mesh_old(self, raw_boundary=True,
1597                    remove_holes=False,
1598                    smooth_indents=False,
1599                    expand_pinch=False):
1600        """
1601        Precon there must be a shape object.
1602        """
1603        self.shape.set_boundary_type(raw_boundary=raw_boundary,
1604                                 remove_holes=remove_holes,
1605                                 smooth_indents=smooth_indents,
1606                                 expand_pinch=expand_pinch)
1607        boundary_segs = self.shape.get_boundary()
1608        #print "boundary_segs",boundary_segs
1609        segs2delete = self.alphaUserSegments
1610
1611        #FIXME(DSG-DSG) this algorithm needs comments
1612        #FIXME(DSG-DSG) can it be sped up?  It's slow
1613        new_segs = []
1614        alpha_segs = []
1615        user_segs = []
1616        for seg in boundary_segs:
1617            v1 = self.userVertices[int(seg[0])]
1618            v2 = self.userVertices[int(seg[1])]
1619            alpha_seg = self.representedAlphaUserSegment(v1, v2)
1620            user_seg = self.representedUserSegment(v1, v2)
1621            #DSG!!!
1622            assert not(not (alpha_seg == None) and not (user_seg == None))
1623            if not alpha_seg == None:
1624                alpha_segs.append(alpha_seg)
1625            elif not user_seg  == None:
1626                user_segs.append(user_seg)
1627            else:
1628                unique_seg = Segment(v1, v2)
1629                new_segs.append(unique_seg)
1630               
1631            for seg in alpha_segs:
1632                try:
1633                    segs2delete.remove(seg)
1634                except:
1635                    pass
1636       
1637        self.alphaUserSegments = []
1638        self.alphaUserSegments.extend(new_segs)
1639        self.alphaUserSegments.extend(alpha_segs)
1640
1641        optimum_alpha = self.shape.get_alpha()
1642        # need to draw newsegs
1643        return new_segs, segs2delete, optimum_alpha
1644   
1645    def representedAlphaUserSegment(self, v1,v2):
1646        identicalSegs= [x for x in self.alphaUserSegments \
1647                        if (x.vertices[0] == v1 and x.vertices[1] == v2)
1648        or (x.vertices[0] == v2 and x.vertices[1] == v1) ]
1649
1650        if identicalSegs == []:
1651            return None
1652        else:
1653            # Only return the first one.
1654            return identicalSegs[0]
1655   
1656    def representedUserSegment(self, v1,v2):
1657        identicalSegs= [x for x in self.userSegments \
1658                        if (x.vertices[0] == v1 and x.vertices[1] == v2)
1659        or (x.vertices[0] == v2 and x.vertices[1] == v1) ]
1660
1661        if identicalSegs == []:
1662            return None
1663        else:
1664            # Only return the first one.
1665            return identicalSegs[0]
1666       
1667    def joinVertices(self):
1668        """
1669        Return list of segments connecting all userVertices
1670        in the order they were given
1671       
1672        Precon: There must be 3 or more vertices in the userVertices structure
1673        """
1674
1675        newsegs = []
1676       
1677        v1 = self.userVertices[0]
1678        for v2 in self.userVertices[1:]:
1679            if self.isUserSegmentNew(v1,v2):           
1680                newseg = Segment(v1, v2)
1681                newsegs.append(newseg)
1682            v1 = v2
1683
1684        #Connect last point to the first
1685        v2 = self.userVertices[0]       
1686        if self.isUserSegmentNew(v1,v2):           
1687                newseg = Segment(v1, v2)
1688                newsegs.append(newseg)
1689           
1690
1691        #Update list of user segments       
1692        #DSG!!!
1693        self.userSegments.extend(newsegs)               
1694        return newsegs
1695   
1696    def normaliseMesh(self,scale, offset, height_scale):
1697        [xmin, ymin, xmax, ymax] = self.boxsize()
1698        [attmin0, attmax0] = self.maxMinVertAtt(0)
1699        #print "[attmin0, attmax0]" ,[attmin0, attmax0]
1700        [attmin1, attmax1] = self.maxMinVertAtt(1)
1701        #print [xmin, ymin, xmax, ymax]
1702        xrange = xmax - xmin
1703        yrange = ymax - ymin
1704        if xrange > yrange:
1705            min,max = xmin, xmax
1706        else:
1707            min,max = ymin, ymax
1708           
1709        for obj in self.getUserVertices():
1710            obj.x = (obj.x - xmin)/(max- min)*scale + offset
1711            obj.y = (obj.y - ymin)/(max- min)*scale + offset
1712            if len(obj.attributes)  > 0 and attmin0 != attmax0:
1713                obj.attributes[0] = (obj.attributes[0]-attmin0)/ \
1714                                    (attmax0-attmin0)*height_scale
1715            if len(obj.attributes) > 1 and attmin1 != attmax1:
1716                obj.attributes[1] = (obj.attributes[1]-attmin1)/ \
1717                                    (attmax1-attmin1)*height_scale
1718           
1719        for obj in self.getMeshVertices():
1720            obj.x = (obj.x - xmin)/(max- min)*scale + offset
1721            obj.y = (obj.y - ymin)/(max- min)*scale + offset
1722            if len(obj.attributes)  > 0 and attmin0 != attmax0:
1723                obj.attributes[0] = (obj.attributes[0]-attmin0)/ \
1724                                    (attmax0-attmin0)*height_scale
1725            if len(obj.attributes) > 1 and attmin1 != attmax1:
1726                obj.attributes[1] = (obj.attributes[1]-attmin1)/ \
1727                                    (attmax1-attmin1)*height_scale
1728               
1729        for obj in self.getHoles():
1730            obj.x = (obj.x - xmin)/(max- min)*scale + offset
1731            obj.y = (obj.y - ymin)/(max- min)*scale + offset
1732        for obj in self.getRegions():
1733            obj.x = (obj.x - xmin)/(max- min)*scale + offset
1734            obj.y = (obj.y - ymin)/(max- min)*scale + offset
1735        [xmin, ymin, xmax, ymax] = self.boxsize()
1736        #print [xmin, ymin, xmax, ymax]
1737     
1738    def boxsizeVerts(self):
1739        """
1740        Returns a list of verts denoting a box or triangle that contains
1741        verts on the xmin, ymin, xmax and ymax axis.
1742        Structure: list of verts
1743        """
1744        # FIXME dsg!!! large is a hack
1745        #You want the kinds package, part of Numeric:
1746        #In [2]: import kinds
1747       
1748        #In [3]: kinds.default_float_kind.M
1749        #kinds.default_float_kind.MAX         kinds.default_float_kind.MIN
1750    #kinds.default_float_kind.MAX_10_EXP  kinds.default_float_kind.MIN_10_EXP
1751        #kinds.default_float_kind.MAX_EXP     kinds.default_float_kind.MIN_EXP
1752
1753        #In [3]: kinds.default_float_kind.MIN
1754        #Out[3]: 2.2250738585072014e-308
1755
1756        large = 1e100
1757        xmin= large
1758        xmax=-large
1759        ymin= large
1760        ymax=-large
1761        for vertex in self.userVertices:
1762            if vertex.x < xmin:
1763                xmin = vertex.x
1764                xminVert = vertex
1765            if vertex.x > xmax:
1766                xmax = vertex.x
1767                xmaxVert = vertex
1768               
1769            if vertex.y < ymin:
1770                ymin = vertex.y
1771                yminVert = vertex
1772            if vertex.y > ymax:
1773                ymax = vertex.y
1774                ymaxVert = vertex
1775        verts, count = self.removeDuplicatedVertices([xminVert,
1776                                                      xmaxVert,
1777                                                      yminVert,
1778                                                      ymaxVert])
1779         
1780        return verts
1781   
1782    def boxsize(self):
1783        """
1784        Returns a list denoting a box that contains the entire
1785        structure of vertices
1786        Structure: [xmin, ymin, xmax, ymax]
1787        """
1788        # FIXME dsg!!! large is a hack
1789        #You want the kinds package, part of Numeric:
1790        #In [2]: import kinds
1791       
1792        #In [3]: kinds.default_float_kind.M
1793        #kinds.default_float_kind.MAX         kinds.default_float_kind.MIN
1794    #kinds.default_float_kind.MAX_10_EXP  kinds.default_fltesting oat_kind.MIN_10_EXP
1795        #kinds.default_float_kind.MAX_EXP     kinds.default_float_kind.MIN_EXP
1796
1797        #In [3]: kinds.default_float_kind.MIN
1798        #Out[3]: 2.2250738585072014e-308
1799
1800        large = 1e100
1801        xmin= large
1802        xmax=-large
1803        ymin= large
1804        ymax=-large
1805        for vertex in self.userVertices:
1806            if vertex.x < xmin:
1807                xmin = vertex.x
1808            if vertex.x > xmax:
1809                xmax = vertex.x
1810               
1811            if vertex.y < ymin:
1812                ymin = vertex.y
1813            if vertex.y > ymax:
1814                ymax = vertex.y
1815        return [xmin, ymin, xmax, ymax]
1816 
1817    def maxMinVertAtt(self, iatt):
1818        """
1819        Returns a list denoting a box that contains the entire structure
1820        of vertices
1821        Structure: [xmin, ymin, xmax, ymax]
1822        """
1823        # FIXME dsg!!! large is a hacktesting
1824        #You want the kinds package, part of Numeric:
1825        #In [2]: import kinds
1826       
1827        #In [3]: kinds.default_float_kind.M
1828        #kinds.default_float_kind.MAX         kinds.default_float_kind.MIN
1829    #kinds.default_float_kind.MAX_10_EXP  kinds.default_float_kind.MIN_10_EXP
1830        #kinds.default_float_kind.MAX_EXP     kinds.default_float_kind.MIN_EXP
1831
1832        #In [3]: kinds.default_float_kind.MIN
1833        #Out[3]: 2.2250738585072014e-308
1834
1835        large = 1e100
1836        min= large
1837        max=-large
1838        for vertex in self.userVertices:
1839            if len(vertex.attributes) > iatt:
1840                if vertex.attributes[iatt] < min:
1841                    min = vertex.attributes[iatt]
1842                if vertex.attributes[iatt] > max:
1843                    max = vertex.attributes[iatt] 
1844        for vertex in self.meshVertices:
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        return [min, max]
1851   
1852    def scaleoffset(self, WIDTH, HEIGHT):
1853        """
1854        Returns a list denoting the scale and offset terms that need to be
1855        applied when converting  mesh co-ordinates onto grid co-ordinates
1856        Structure: [scale, xoffset, yoffset]
1857        """   
1858        OFFSET = 0.05*min([WIDTH, HEIGHT])
1859        [xmin, ymin, xmax, ymax] = self.boxsize()
1860        SCALE = min([0.9*WIDTH, 0.9*HEIGHT])/max([xmax-xmin, ymax-ymin])
1861       
1862        if SCALE*xmin < OFFSET:
1863            xoffset = abs(SCALE*xmin) + OFFSET
1864        if SCALE*xmax > WIDTH - OFFSET:
1865            xoffset= -(SCALE*xmax - WIDTH + OFFSET)
1866        if SCALE*ymin < OFFSET:
1867            b = abs(SCALE*ymin)+OFFSET
1868        if SCALE*ymax > HEIGHT-OFFSET:
1869            b = -(SCALE*ymax - HEIGHT + OFFSET)
1870        yoffset = HEIGHT - b
1871        return [SCALE, xoffset, yoffset]
1872
1873         
1874    def exportASCIIobj(self,ofile):
1875        """
1876        export a file, ofile, with the format
1877         lines:  v <x> <y> <first attribute>
1878        f <vertex #>  <vertex #> <vertex #> (of the triangles)
1879        """
1880        fd = open(ofile,'w')
1881        self.writeASCIIobj(fd)   
1882        fd.close()
1883
1884
1885    def writeASCIIobj(self,fd):
1886        fd.write(" # Triangulation as an obj file\n")
1887        numVert = str(len(self.meshVertices))
1888       
1889        index1 = 1 
1890        for vert in self.meshVertices:
1891            vert.index1 = index1
1892            index1 += 1
1893           
1894            fd.write("v "
1895                     + str(vert.x) + " "
1896                     + str(vert.y) + " "
1897                     + str(vert.attributes[0]) + "\n")
1898           
1899        for tri in self.meshTriangles:
1900            fd.write("f "
1901                     + str(tri.vertices[0].index1) + " " 
1902                     + str(tri.vertices[1].index1) + " " 
1903                     + str(tri.vertices[2].index1) + "\n")
1904           
1905    def exportASCIIsegmentoutlinefile(self,ofile):
1906        """Write the boundary user mesh info, eg
1907         vertices that are connected to segments,
1908         segments
1909        """
1910       
1911        verts = {}
1912        for seg in self.getUserSegments():
1913            verts[seg.vertices[0]] = seg.vertices[0]
1914            verts[seg.vertices[1]] = seg.vertices[1]
1915        meshDict = self.Mesh2IOOutlineDict(userVertices=verts.values())
1916        load_mesh.loadASCII.export_mesh_file(ofile,meshDict)
1917       
1918        # exportASCIImeshfile   - this function is used
1919    def export_mesh_file(self,ofile):
1920        """
1921        export a file, ofile, with the format
1922        """
1923       
1924        dict = self.Mesh2IODict()
1925        load_mesh.loadASCII.export_mesh_file(ofile,dict)
1926
1927    def exportPointsFile(self,ofile):
1928        """
1929        export a points (.xya or .pts)  file, ofile.
1930       
1931        """
1932       
1933        mesh_dict = self.Mesh2IODict()
1934        point_dict = {}
1935        point_dict['attributelist'] = {} #this will need to be expanded..
1936                                         # if attributes are brought back in.
1937        point_dict['geo_reference'] = self.geo_reference
1938        if mesh_dict['vertices'] == []:
1939            point_dict['pointlist'] = mesh_dict['points']
1940        else:
1941            point_dict['pointlist'] = mesh_dict['vertices']
1942
1943        load_mesh.loadASCII.export_points_file(ofile,point_dict)
1944
1945       
1946########### IO CONVERTERS ##################
1947        """
1948        The dict fromat for IO with .tsh files is;
1949        (the triangulation)
1950        vertices: [[x1,y1],[x2,y2],...] (lists of doubles)
1951        vertex_attributes: [[a11,a12,...],[a21,a22],...] (lists of doubles)
1952        vertex_attribute_titles:[A1Title, A2Title ...] (A list of strings)
1953        segments: [[v1,v2],[v3,v4],...] (lists of integers)
1954        segment_tags : [tag,tag,...] list of strings
1955        triangles : [(v1,v2,v3), (v4,v5,v6),....] lists of points
1956        triangle tags: [s1,s2,...] A list of strings
1957        triangle neighbors: [[t1,t2,t3], [t4,t5,t6],..] lists of triangles
1958       
1959        (the outline)   
1960        points: [[x1,y1],[x2,y2],...] (lists of doubles)
1961        point_attributes: [[a11,a12,...],[a21,a22],...] (lists of doubles)
1962        outline_segments: [[point1,point2],[p3,p4],...] (lists of integers)
1963        outline_segment_tags : [tag1,tag2,...] list of strings
1964        holes : [[x1,y1],...](List of doubles, one inside each hole region)
1965        regions : [ [x1,y1],...] (List of 4 doubles)
1966        region_tags : [tag1,tag2,...] (list of strings)
1967        region_max_areas: [ma1,ma2,...] (A list of doubles)
1968        {Convension: A -ve max area means no max area}
1969       
1970        """
1971     
1972
1973                               
1974    def Mesh2IODict(self):
1975        """
1976        Convert the triangulation and outline info of a mesh to a dictionary
1977        structure
1978        """
1979        dict = self.Mesh2IOTriangulationDict()
1980        dict_mesh = self.Mesh2IOOutlineDict()
1981        for element in dict_mesh.keys():
1982            dict[element] = dict_mesh[element]
1983
1984        # add the geo reference
1985        dict['geo_reference'] = self.geo_reference
1986        return dict
1987   
1988    def Mesh2IOTriangulationDict(self):
1989        """
1990        Convert the Mesh to a dictionary of lists describing the
1991        triangulation variables;
1992       
1993        Used to produce .tsh file
1994        """
1995       
1996        meshDict = {}       
1997        vertices=[]
1998        vertex_attributes=[]
1999
2000        self.maxVertexIndex=0
2001        for vertex in self.meshVertices:
2002            vertex.index = self.maxVertexIndex
2003            vertices.append([vertex.x,vertex.y])
2004            vertex_attributes.append(vertex.attributes)           
2005            self.maxVertexIndex += 1
2006
2007        meshDict['vertices'] = vertices
2008        meshDict['vertex_attributes'] = vertex_attributes
2009        meshDict['vertex_attribute_titles'] = self.attributeTitles
2010        #segments
2011        segments=[]
2012        segment_tags=[]
2013        for seg in self.meshSegments:
2014            segments.append([seg.vertices[0].index,seg.vertices[1].index])
2015            segment_tags.append(seg.tag)
2016        meshDict['segments'] =segments
2017        meshDict['segment_tags'] =segment_tags
2018
2019        # Make sure that the indexation is correct
2020        index = 0
2021        for tri in self.meshTriangles:
2022            tri.index = index           
2023            index += 1
2024
2025        triangles = []
2026        triangle_tags = []
2027        triangle_neighbors = []
2028        for tri in self.meshTriangles: 
2029            triangles.append([tri.vertices[0].index,
2030                              tri.vertices[1].index,
2031                              tri.vertices[2].index]) 
2032            triangle_tags.append(tri.attribute)
2033            neighborlist = [-1,-1,-1]
2034            for neighbor,index in map(None,tri.neighbors,
2035                                      range(len(tri.neighbors))):
2036                if neighbor:   
2037                    neighborlist[index] = neighbor.index
2038            triangle_neighbors.append(neighborlist)
2039       
2040        meshDict['triangles'] = triangles
2041        meshDict['triangle_tags'] = triangle_tags
2042        meshDict['triangle_neighbors'] = triangle_neighbors
2043       
2044        #print "mesh.Mesh2IOTriangulationDict*)*)"
2045        #print meshDict
2046        #print "mesh.Mesh2IOTriangulationDict*)*)"
2047
2048        return meshDict
2049
2050                                                     
2051    def Mesh2IOOutlineDict(self, userVertices=None,
2052                        userSegments=None,
2053                        holes=None,
2054                        regions=None):
2055        """
2056        Convert the mesh outline to a dictionary of the lists needed for the
2057        triang module;
2058       
2059        Note, this adds an index attribute to the user Vertex objects.
2060
2061        Used to produce .tsh file and output to triangle
2062        """
2063        if userVertices is None:
2064            userVertices = self.getUserVertices()
2065        if userSegments is None:
2066            userSegments = self.getUserSegments()
2067        if holes is None:
2068            holes = self.getHoles()
2069        if regions is None:
2070            regions = self.getRegions()
2071           
2072        meshDict = {}
2073        #print "userVertices",userVertices
2074        #print "userSegments",userSegments
2075        pointlist=[]
2076        pointattributelist=[]
2077        index = 0
2078        for vertex in userVertices:
2079            vertex.index = index
2080            pointlist.append([vertex.x,vertex.y])
2081            pointattributelist.append(vertex.attributes)
2082           
2083            index += 1
2084        meshDict['points'] = pointlist
2085        meshDict['point_attributes'] = pointattributelist
2086
2087        segmentlist=[]
2088        segmenttaglist=[]
2089        for seg in userSegments:
2090            segmentlist.append([seg.vertices[0].index,seg.vertices[1].index])
2091            segmenttaglist.append(seg.tag)
2092        meshDict['outline_segments'] =segmentlist
2093        meshDict['outline_segment_tags'] =segmenttaglist
2094       
2095        holelist=[]
2096        for hole in holes:
2097            holelist.append([hole.x,hole.y]) 
2098        meshDict['holes'] = holelist
2099       
2100        regionlist=[]
2101        regiontaglist = []
2102        regionmaxarealist = []
2103        for region in regions:
2104            regionlist.append([region.x,region.y])
2105            regiontaglist.append(region.getTag())
2106           
2107            if (region.getMaxArea() != None): 
2108                regionmaxarealist.append(region.getMaxArea())
2109            else:
2110                regionmaxarealist.append( load_mesh.loadASCII.NOMAXAREA)
2111        meshDict['regions'] = regionlist
2112        meshDict['region_tags'] = regiontaglist
2113        meshDict['region_max_areas'] = regionmaxarealist
2114        #print "*(*("
2115        #print meshDict
2116        #print meshDict['regionlist']
2117        #print "*(*("
2118        return meshDict
2119
2120    def IOTriangulation2Mesh(self, genDict):
2121        """
2122        Set the mesh attributes given an tsh IO dictionary
2123        """
2124        #Clear the current generated mesh values
2125        self.meshTriangles=[]
2126        self.attributeTitles=[]
2127        self.meshSegments=[]
2128        self.meshVertices=[]
2129
2130        #print "mesh.setTriangulation@#@#@#"
2131        #print genDict
2132        #print "@#@#@#"
2133
2134        self.maxVertexIndex = 0
2135        for point in genDict['vertices']:
2136            v=Vertex(point[0], point[1])
2137            v.index =  self.maxVertexIndex
2138            self.maxVertexIndex +=1
2139            self.meshVertices.append(v)
2140
2141        self.attributeTitles = genDict['vertex_attribute_titles']
2142
2143        index = 0
2144        for seg,tag in map(None,genDict['segments'],genDict['segment_tags']):
2145            segObject = Segment( self.meshVertices[int(seg[0])],
2146                           self.meshVertices[int(seg[1])], tag = tag )
2147            segObject.index = index
2148            index +=1
2149            self.meshSegments.append(segObject)
2150
2151        index = 0
2152        for triangle in genDict['triangles']:
2153            tObject =Triangle( self.meshVertices[int(triangle[0])],
2154                        self.meshVertices[int(triangle[1])],
2155                        self.meshVertices[int(triangle[2])] )
2156            tObject.index = index
2157            index +=1
2158            self.meshTriangles.append(tObject)
2159
2160        index = 0
2161        for att in genDict['triangle_tags']:
2162            if att == []:
2163                self.meshTriangles[index].setAttribute("")
2164            else:
2165                self.meshTriangles[index].setAttribute(att)
2166            index += 1
2167           
2168        index = 0
2169        for att in genDict['vertex_attributes']:
2170            if att == None:
2171                self.meshVertices[index].setAttributes([])
2172            else:
2173                self.meshVertices[index].setAttributes(att)
2174            index += 1
2175   
2176        index = 0
2177        for triangle in genDict['triangle_neighbors']:
2178            # Build a list of triangle object neighbors
2179            ObjectNeighbor = []
2180            for neighbor in triangle:
2181                if ( neighbor != -1):
2182                    ObjectNeighbor.append(self.meshTriangles[int(neighbor)])
2183                else:
2184                    ObjectNeighbor.append(None)
2185            self.meshTriangles[index].setNeighbors(ObjectNeighbor[0],
2186                                                   ObjectNeighbor[1],
2187                                                   ObjectNeighbor[2])
2188            index += 1
2189
2190
2191    def IOOutline2Mesh(self, genDict):
2192        """
2193        Set the outline (user Mesh attributes) given a IO tsh dictionary
2194       
2195        mesh is an instance of a mesh object
2196        """
2197        #Clear the current user mesh values
2198        self.clearUserSegments()
2199        self.userVertices=[]
2200        self.Holes=[]
2201        self.Regions=[]
2202
2203        #print "mesh.IOOutline2Mesh@#@#@#"
2204        #print "genDict",genDict
2205        #print "@#@#@#"
2206       
2207        #index = 0
2208        for point in genDict['points']:
2209            v=Vertex(point[0], point[1])
2210            #v.index = index
2211            #index +=1
2212            self.userVertices.append(v)
2213
2214        #index = 0
2215        for seg,tag in map(None,genDict['outline_segments'],
2216                           genDict['outline_segment_tags']):
2217
2218            segObject = Segment( self.userVertices[int(seg[0])],
2219                           self.userVertices[int(seg[1])], tag = tag )
2220            #segObject.index = index
2221            #index +=1
2222            self.userSegments.append(segObject)
2223
2224# Remove the loading of attribute info.
2225# Have attribute info added using least_squares in pyvolution
2226#         index = 0
2227#         for att in genDict['point_attributes']:
2228#             if att == None:
2229#                 self.userVertices[index].setAttributes([])
2230#             else:
2231#                 self.userVertices[index].setAttributes(att)
2232#            index += 1
2233       
2234        #index = 0
2235        for point in genDict['holes']:
2236            h=Hole(point[0], point[1])
2237            #h.index = index
2238            #index +=1
2239            self.holes.append(h)
2240
2241        #index = 0
2242        for reg,att,maxArea in map(None,
2243                                   genDict['regions'],
2244                                   genDict['region_tags'],
2245                                   genDict['region_max_areas']):
2246            if maxArea > 0:  # maybe I should ref NOMAXAREA? Prob' not though
2247                Object = Region( reg[0],
2248                                 reg[1],
2249                                 tag = att,
2250                                 maxArea = maxArea)
2251            else:
2252                Object = Region( reg[0],
2253                                 reg[1],
2254                                 tag = att)
2255               
2256            #Object.index = index
2257            #index +=1
2258            self.regions.append(Object)
2259 
2260############################################
2261
2262       
2263    def refineSet(self,setName):
2264        Triangles = self.sets[self.setID[setName]]
2265        Refine(self,Triangles)
2266
2267    def selectAllTriangles(self):
2268        A=[]
2269        A.extend(self.meshTriangles)
2270        if not('All' in self.setID.keys()):
2271            self.setID['All']=len(self.sets)
2272            self.sets.append(A)
2273        else:
2274            self.sets[self.setID['All']]=A
2275        return 'All'
2276        # and objectIDs
2277
2278
2279    def clearSelection(self):
2280        A = []
2281        if not('None' in self.setID.keys()):
2282            self.setID['None']=len(self.sets)
2283            self.sets.append(A)
2284        return 'None'
2285
2286    def drawSet(self,canvas,setName,SCALE,colour=SET_COLOUR):
2287    #FIXME Draws over previous triangles - may bloat canvas
2288        Triangles = self.sets[self.setID[setName]]
2289        for triangle in Triangles:
2290            triangle.draw(canvas,1,
2291                          scale = SCALE,
2292                          colour = colour)
2293           
2294    def undrawSet(self,canvas,setName,SCALE,colour='green'):
2295    #FIXME Draws over previous lines - may bloat canvas
2296        Triangles = self.sets[self.setID[setName]]
2297        for triangle in Triangles:
2298            triangle.draw(canvas,1,
2299                          scale = SCALE,
2300                          colour = colour)
2301
2302    def weed(self,Vertices,Segments):
2303        #Depreciated
2304        #weed out existing duplicates
2305        print 'len(self.getUserSegments())'
2306        print len(self.getUserSegments())
2307        print 'len(self.getUserVertices())'
2308        print len(self.getUserVertices())
2309
2310        point_keys = {}
2311        for vertex in Vertices:
2312            point = (vertex.x,vertex.y)
2313            point_keys[point]=vertex
2314        #inlined would looks very ugly
2315
2316        line_keys = {}
2317        for segment in Segments:
2318            vertex1 = segment.vertices[0]
2319            vertex2 = segment.vertices[1]
2320            point1 = (vertex1.x,vertex1.y)
2321            point2 = (vertex2.x,vertex2.y)
2322            segment.vertices[0]=point_keys[point1]
2323            segment.vertices[1]=point_keys[point2]
2324            vertex1 = segment.vertices[0]
2325            vertex2 = segment.vertices[1]
2326            point1 = (vertex1.x,vertex1.y)
2327            point2 = (vertex2.x,vertex2.y)
2328            line1 = (point1,point2)
2329            line2 = (point2,point1)
2330            if not (line_keys.has_key(line1) \
2331                 or line_keys.has_key(line2)):
2332                 line_keys[line1]=segment
2333        Vertices=point_keys.values()
2334        Segments=line_keys.values()
2335        return Vertices,Segments
2336
2337    def segs_to_dict(self,segments):
2338        dict={}
2339        for segment in segments:
2340            vertex1 = segment.vertices[0]
2341            vertex2 = segment.vertices[1]
2342            point1 = (vertex1.x,vertex1.y)
2343            point2 = (vertex2.x,vertex2.y)
2344            line = (point1,point2)
2345            dict[line]=segment
2346        return dict
2347
2348    def seg2line(self,s):
2349        return ((s.vertices[0].x,s.vertices[0].y,)\
2350                (s.vertices[1].x,s.vertices[1].y))
2351
2352    def line2seg(self,line,tag=None):
2353        point0 = self.point2ver(line[0])
2354        point1 = self.point2ver(line[1])
2355        return Segment(point0,point1,tag=tag)
2356
2357    def ver2point(self,vertex):
2358        return (vertex.x,vertex.y)
2359
2360    def point2ver(self,point):
2361        return Vertex(point[0],point[1])
2362
2363    def smooth_polySet(self,min_radius=0.05):
2364        #for all pairs of connecting segments:
2365        #    propose a new segment that replaces the 2
2366
2367        #    If the difference between the new segment
2368        #    and the old lines is small: replace the
2369        #    old lines.
2370
2371        seg2line = self.seg2line
2372        ver2point= self.ver2point
2373        line2seg = self.line2seg
2374        point2ver= self.point2ver
2375
2376        #create dictionaries of lines -> segments
2377        userSegments = self.segs_to_dict(self.userSegments)
2378        alphaSegments = self.segs_to_dict(self.alphaUserSegments)
2379
2380        #lump user and alpha segments
2381        for key in alphaSegments.keys():
2382            userSegments[key]=alphaSegments[key]
2383
2384        #point_keys = tuple -> vertex
2385        #userVertices = vertex -> [line,line] - lines from that node
2386        point_keys = {}
2387        userVertices={}
2388        for vertex in self.getUserVertices():
2389            point = ver2point(vertex)
2390            if not point_keys.has_key(point):
2391                point_keys[point]=vertex
2392                userVertices[vertex]=[]
2393        for key in userSegments.keys():
2394            line = key
2395            point_0 = key[0]
2396            point_1 = key[1]
2397            userVertices[point_keys[point_0]].append(line)
2398            userVertices[point_keys[point_1]].append(line)
2399
2400        for point in point_keys.keys():
2401            try:
2402            #removed keys can cause keyerrors
2403                vertex = point_keys[point]
2404                lines = userVertices[vertex]
2405   
2406                #if there are 2 lines on the node
2407                if len(lines)==2:
2408                    line_0 = lines[0]
2409                    line_1 = lines[1]
2410   
2411                    #if the tags are the the same on the 2 lines
2412                    if userSegments[line_0].tag == userSegments[line_1].tag:
2413                        tag = userSegments[line_0].tag
2414   
2415                        #point_a is one of the next nodes, point_b is the other
2416                        if point==line_0[0]:
2417                            point_a = line_0[1]
2418                        if point==line_0[1]:
2419                            point_a = line_0[0]
2420                        if point==line_1[0]:
2421                            point_b = line_1[1]
2422                        if point==line_1[1]:
2423                            point_b = line_1[0]
2424   
2425   
2426                        #line_2 is proposed
2427                        line_2 = (point_a,point_b)
2428
2429                        #calculate the area of the triangle between
2430                        #the two existing segments and the proposed
2431                        #new segment
2432                        ax = point_a[0]
2433                        ay = point_a[1]
2434                        bx = point_b[0]
2435                        by = point_b[1]
2436                        cx = point[0]
2437                        cy = point[1]
2438                        area=abs((bx*ay-ax*by)+(cx*by-bx*cy)+(ax*cy-cx*ay))/2
2439
2440                        #calculate the perimeter
2441                        len_a =  ((cx-bx)**2+(cy-by)**2)**0.5 
2442                        len_b =  ((ax-cx)**2+(ay-cy)**2)**0.5 
2443                        len_c =  ((bx-ax)**2+(by-ay)**2)**0.5 
2444                        perimeter = len_a+len_b+len_c
2445
2446                        #calculate the radius
2447                        r = area/(2*perimeter)
2448
2449                        #if the radius is small: then replace the existing
2450                        #segments with the new one
2451                        if r < min_radius:
2452                            if len_c < min_radius: append = False
2453                            else: append = True
2454                            #if the new seg is also time, don't add it
2455                            if append:
2456                                segment = self.line2seg(line_2,tag=tag)
2457
2458                            list_a=userVertices[point_keys[point_a]]
2459                            list_b=userVertices[point_keys[point_b]]
2460
2461                            if line_0 in list_a:
2462                                list_a.remove(line_0)
2463                            else:
2464                                list_a.remove(line_1)
2465
2466                            if line_0 in list_b:
2467                                list_b.remove(line_0)
2468                            else:
2469                                list_b.remove(line_1)
2470
2471                            if append:
2472                                list_a.append(line_2)
2473                                list_b.append(line_2)
2474                            else:
2475                                if len(list_a)==0:
2476                                    userVertices.pop(point_keys[point_a])
2477                                    point_keys.pop(point_a)
2478                                if len(list_b)==0:
2479                                    userVertices.pop(point_keys[point_b])
2480                                    point_keys.pop(point_b)
2481
2482                            userVertices.pop(point_keys[point])
2483                            point_keys.pop(point)
2484                            userSegments.pop(line_0)
2485                            userSegments.pop(line_1)
2486
2487                            if append:
2488                                userSegments[line_2]=segment
2489            except:
2490                pass
2491
2492        #self.userVerticies = userVertices.keys()
2493        #self.userSegments = []
2494        #for key in userSegments.keys():
2495        #    self.userSegments.append(userSegments[key])
2496        #self.alphaUserSegments = []
2497
2498        self.userVerticies = []
2499        self.userSegments = []
2500        self.alphaUserSegments = []
2501
2502        return userVertices,userSegments,alphaSegments
2503
2504    def triangles_to_polySet(self,setName):
2505        #self.smooth_polySet()
2506
2507        seg2line = self.seg2line
2508        ver2point= self.ver2point
2509        line2seg = self.line2seg
2510        point2ver= self.point2ver
2511
2512        from Numeric import array,allclose
2513        #turn the triangles into a set
2514        Triangles = self.sets[self.setID[setName]]
2515        Triangles_dict = {}
2516        for triangle in Triangles:
2517            Triangles_dict[triangle]=None
2518 
2519
2520        #create a dict of points to vertexes (tuple -> object)
2521        #also create a set of vertexes (object -> True)
2522        point_keys = {}
2523        userVertices={}
2524        for vertex in self.getUserVertices():
2525            point = ver2point(vertex)
2526            if not point_keys.has_key(point):
2527                point_keys[point]=vertex
2528                userVertices[vertex]=True
2529
2530        #create a dict of lines to segments (tuple -> object)
2531        userSegments = self.segs_to_dict(self.userSegments)
2532        #append the userlines in an affine linespace
2533        affine_lines = Affine_Linespace()
2534        for line in userSegments.keys():
2535            affine_lines.append(line)
2536        alphaSegments = self.segs_to_dict(self.alphaUserSegments)
2537        for line in alphaSegments.keys():
2538            affine_lines.append(line)
2539       
2540        for triangle in Triangles:
2541            for i in (0,1,2):
2542                #for every triangles neighbour:
2543                if not Triangles_dict.has_key(triangle.neighbors[i]):
2544                #if the neighbour is not in the set:
2545                    a = triangle.vertices[i-1]
2546                    b = triangle.vertices[i-2]
2547                    #Get possible matches:
2548                    point_a = ver2point(a)
2549                    point_b = ver2point(b)
2550                    midpoint = ((a.x+b.x)/2,(a.y+b.y)/2)
2551                    line = (point_a,point_b)
2552                    tag = None
2553
2554
2555                    #this bit checks for matching lines
2556                    possible_lines = affine_lines[line] 
2557                    possible_lines = unique(possible_lines)
2558                    found = 0                           
2559                    for user_line in possible_lines:
2560                        if self.point_on_line(midpoint,user_line):
2561                            found+=1
2562                            assert found<2
2563                            if userSegments.has_key(user_line):
2564                                parent_segment = userSegments.pop(user_line)
2565                            if alphaSegments.has_key(user_line):
2566                                parent_segment = alphaSegments.pop(user_line)
2567                            tag = parent_segment.tag
2568                            offspring = [line]
2569                            offspring.extend(self.subtract_line(user_line,line))
2570                            affine_lines.remove(user_line)
2571                            for newline in offspring:
2572                                line_vertices = []
2573                                for point in newline:
2574                                    if point_keys.has_key(point):
2575                                        vert = point_keys[point]
2576                                    else:
2577                                        vert = Vertex(point[0],point[1])
2578                                        userVertices[vert]=True
2579                                        point_keys[point]=vert
2580                                    line_vertices.append(vert)
2581                                segment = Segment(line_vertices[0],line_vertices[1],tag)
2582                                userSegments[newline]=segment
2583                                affine_lines.append(newline)
2584                            #break
2585                    assert found<2
2586
2587
2588
2589                    #if no matching lines
2590                    if not found:
2591                        line_vertices = []
2592                        for point in line:
2593                            if point_keys.has_key(point):
2594                                vert = point_keys[point]
2595                            else:
2596                                vert = Vertex(point[0],point[1])
2597                                userVertices[vert]=True
2598                                point_keys[point]=vert
2599                            line_vertices.append(vert)
2600                        segment = Segment(line_vertices[0],line_vertices[1],tag)
2601                        userSegments[line]=segment
2602                        affine_lines.append(line)
2603       
2604        self.userVerticies = []
2605        self.userSegments = []
2606        self.alphaUserSegments = []
2607
2608        return userVertices,userSegments,alphaSegments
2609
2610    def subtract_line(self,parent,child):
2611        #Subtracts child from parent
2612        #Requires that the child is a
2613        #subline of parent to work.
2614
2615        from Numeric import allclose,dot,array
2616        A= parent[0]
2617        B= parent[1]
2618        a = child[0]
2619        b = child[1]
2620
2621        A_array = array(parent[0])
2622        B_array = array(parent[1])
2623        a_array   = array(child[0])
2624        b_array   = array(child[1])
2625
2626        assert not A == B
2627        assert not a == b
2628
2629        answer = []
2630
2631        #if the new line does not share a
2632        #vertex with the old one
2633        if not (allclose(A_array,a_array)\
2634             or allclose(B_array,b_array)\
2635             or allclose(A_array,b_array)\
2636             or allclose(a_array,B_array)):
2637            if dot(A_array-a_array,A_array-a_array) \
2638            < dot(A_array-b_array,A_array-b_array):
2639                sibling1 = (A,a)
2640                sibling2 = (B,b)
2641                return [sibling1,sibling2]
2642            else:
2643                sibling1 = (A,b)
2644                sibling2 = (B,a)
2645                return [sibling1,sibling2]
2646
2647        elif allclose(A_array,a_array):
2648            if allclose(B_array,b_array):
2649                return []
2650            else:
2651                sibling = (b,B)
2652                return [sibling]
2653        elif allclose(B_array,b_array):
2654            sibling = (a,A)
2655            return [sibling]
2656
2657        elif allclose(A_array,b_array):
2658            if allclose(B,a):
2659                return []
2660            else:
2661                sibling = (a,B)
2662                return [sibling]
2663        elif allclose(a_array,B_array):
2664            sibling = (b,A)
2665            return [sibling]
2666
2667    def point_on_line(self,point,line):
2668        #returns true within a tolerance of 3 degrees
2669        x=point[0]
2670        y=point[1]
2671        x0=line[0][0]
2672        x1=line[1][0]
2673        y0=line[0][1]
2674        y1=line[1][1]
2675        from Numeric import array, dot, allclose
2676        from math import sqrt
2677        tol = 3. #DEGREES
2678        tol = tol*3.1415/180
2679
2680        a = array([x - x0, y - y0]) 
2681        a_normal = array([a[1], -a[0]])     
2682        len_a_normal = sqrt(sum(a_normal**2)) 
2683
2684        b = array([x1 - x0, y1 - y0])         
2685        len_b = sqrt(sum(b**2)) 
2686   
2687        if abs(dot(a_normal, b)/(len_b*len_a_normal))< tol:
2688            #Point is somewhere on the infinite extension of the line
2689
2690            len_a = sqrt(sum(a**2))                                     
2691            if dot(a, b) >= 0 and len_a <= len_b:
2692               return True
2693            else:   
2694               return False
2695        else:
2696          return False
2697
2698    def line_length(self,line):
2699        x0=line[0][0]
2700        x1=line[1][0]
2701        y0=line[0][1]
2702        y1=line[1][1]
2703        return ((x1-x0)**2-(y1-y0)**2)**0.5     
2704
2705    def threshold(self,setName,min=None,max=None,attribute_name = 'elevation'):
2706        """
2707        threshold using  d
2708        """
2709        triangles = self.sets[self.setID[setName]]
2710        A = []
2711
2712        if attribute_name in self.attributeTitles:
2713            i = self.attributeTitles.index(attribute_name)
2714        else: i = -1#no attribute
2715        if not max == None:
2716            for t in triangles:
2717                if (min<self.av_att(t,i)<max):
2718                    A.append(t)
2719        else:
2720            for t in triangles:
2721                if (min<self.av_att(t,i)):
2722                    A.append(t)
2723        self.sets[self.setID[setName]] = A
2724
2725    def general_threshold(self,setName,min=None,max=None\
2726              ,attribute_name = 'elevation',function=None):
2727        """
2728        Thresholds the triangles
2729        """
2730        from visual.graph import arange,ghistogram,color as colour
2731        triangles = self.sets[self.setID[setName]]
2732        A = []
2733        data=[]
2734        #data is for the graph
2735
2736        if attribute_name in self.attributeTitles:
2737            i = self.attributeTitles.index(attribute_name)
2738        else: i = -1
2739        if not max == None:
2740            for t in triangles:
2741                value=function(t,i)
2742                if (min<value<max):
2743                    A.append(t)
2744                data.append(value)
2745        else:
2746            for t in triangles:
2747                value=function(t,i)
2748                if (min<value):
2749                    A.append(t)
2750                data.append(value)
2751        self.sets[self.setID[setName]] = A
2752
2753        if self.visualise_graph:
2754            if len(data)>0:
2755                max=data[0]
2756                min=data[0]
2757                for value in data:
2758                    if value > max:
2759                        max = value
2760                    if value < min:
2761                        min = value
2762
2763                inc = (max-min)/100
2764
2765                histogram = ghistogram(bins=arange(min,max,inc),\
2766                             color = colour.red)
2767                histogram.plot(data=data)
2768       
2769    def av_att(self,triangle,i):
2770        if i==-1: return 1
2771        else:
2772            #evaluates the average attribute of the vertices of a triangle.
2773            V = triangle.getVertices()
2774            a0 = (V[0].attributes[i])
2775            a1 = (V[1].attributes[i])
2776            a2 = (V[2].attributes[i])
2777            return (a0+a1+a2)/3
2778
2779    def Courant_ratio(self,triangle,index):
2780        """
2781        Uses the courant threshold
2782        """
2783        e = self.av_att(triangle,index)
2784        A = triangle.calcArea()
2785        P = triangle.calcP()
2786        r = A/(2*P)
2787        e = max(0.1,abs(e))
2788        return r/e**0.5
2789
2790    def Gradient(self,triangle,index):
2791        V = triangle.vertices
2792        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]
2793        grad_x,grad_y = gradient(x0, y0, x1, y1, x2, y2, q0, q1, q2)
2794        if ((grad_x**2)+(grad_y**2))**(0.5)<0:
2795            print ((grad_x**2)+(grad_y**2))**(0.5)
2796        return ((grad_x**2)+(grad_y**2))**(0.5)
2797   
2798
2799    def append_triangle(self,triangle):
2800        self.meshTriangles.append(triangle)
2801
2802    def replace_triangle(self,triangle,replacement):
2803        i = self.meshTriangles.index(triangle)
2804        self.meshTriangles[i]=replacement
2805        assert replacement in self.meshTriangles
2806
2807def importUngenerateFile(ofile):
2808    """
2809    import a file, ofile, with the format
2810    [poly]
2811    poly format:
2812    First line:  <# of vertices> <x centroid> <y centroid>
2813    Following lines: <x> <y>
2814    last line:  "END"
2815
2816    Note: These are clockwise.
2817    """
2818    fd = open(ofile,'r')
2819    Dict = readUngenerateFile(fd)
2820    fd.close()
2821    return Dict
2822
2823def readUngenerateFile(fd):
2824    """
2825    import a file, ofile, with the format
2826    [poly]
2827    poly format:
2828    First line:  <# of polynomial> <x centroid> <y centroid>
2829    Following lines: <x> <y>
2830    last line:  "END"
2831    """
2832    END_DELIMITER = 'END\n'
2833   
2834    points = []
2835    segments = []
2836   
2837    isEnd = False
2838    line = fd.readline() #not used <# of polynomial> <x> <y>
2839    while not isEnd:
2840        line = fd.readline()
2841        fragments = line.split()
2842        vert = [float(fragments.pop(0)),float(fragments.pop(0))]
2843        points.append(vert)
2844        PreviousVertIndex = len(points)-1
2845        firstVertIndex = PreviousVertIndex
2846       
2847        line = fd.readline() #Read the next line
2848        while line <> END_DELIMITER: 
2849            #print "line >" + line + "<"
2850            fragments = line.split()
2851            vert = [float(fragments.pop(0)),float(fragments.pop(0))]
2852            points.append(vert)
2853            thisVertIndex = len(points)-1
2854            segment = [PreviousVertIndex,thisVertIndex]
2855            segments.append(segment)
2856            PreviousVertIndex = thisVertIndex
2857            line = fd.readline() #Read the next line
2858            i =+ 1
2859        # If the last and first segments are the same,
2860        # Remove the last segment and the last vertex
2861        # then add a segment from the second last vert to the 1st vert
2862        thisVertIndex = len(points)-1
2863        firstVert = points[firstVertIndex]
2864        thisVert = points[thisVertIndex]
2865        #print "firstVert",firstVert
2866        #print "thisVert",thisVert
2867        if (firstVert[0] == thisVert[0] and firstVert[1] == thisVert[1]):
2868            points.pop()
2869            segments.pop()
2870            thisVertIndex = len(points)-1
2871            segments.append([thisVertIndex, firstVertIndex])
2872       
2873        line = fd.readline() # read <# of polynomial> <x> <y> OR END
2874        #print "line >>" + line + "<<"
2875        if line == END_DELIMITER:
2876            isEnd = True
2877   
2878    #print "points", points       
2879    #print "segments", segments
2880    ungenerated_dict = {}
2881    ungenerated_dict['points'] = points
2882    ungenerated_dict['segments'] = segments
2883    return ungenerated_dict
2884
2885def importMeshFromFile(ofile):
2886    """returns a mesh object, made from a .xya/.pts or .tsh/.msh file
2887    Often raises IOError,RuntimeError
2888    """
2889    newmesh = None
2890    if (ofile[-4:]== ".xya" or ofile[-4:]== ".pts"):
2891        dict = load_mesh.loadASCII.import_points_file(ofile)
2892        dict['points'] = dict['pointlist']
2893        dict['outline_segments'] = []
2894        dict['outline_segment_tags'] = []
2895        dict['regions'] = []
2896        dict['region_tags'] = []
2897        dict['region_max_areas'] = []
2898        dict['holes'] = [] 
2899        newmesh= Mesh(geo_reference = dict['geo_reference'])
2900        newmesh.IOOutline2Mesh(dict)
2901        counter = newmesh.removeDuplicatedUserVertices()
2902        if (counter >0):
2903            print "%i duplicate vertices removed from dataset" % (counter)
2904    elif (ofile[-4:]== ".tsh" or ofile[-4:]== ".msh"):
2905        dict = load_mesh.loadASCII.import_mesh_file(ofile)
2906        #print "********"
2907        #print "zq mesh.dict",dict
2908        #print "********"
2909        newmesh= Mesh()
2910        newmesh.IOOutline2Mesh(dict)
2911        newmesh.IOTriangulation2Mesh(dict)
2912    else:
2913        raise RuntimeError
2914   
2915    if dict.has_key('geo_reference') and not dict['geo_reference'] == None:
2916        newmesh.geo_reference = dict['geo_reference']
2917    return newmesh
2918
2919def loadPickle(currentName):
2920    fd = open(currentName)
2921    mesh = pickle.load(fd)
2922    fd.close()
2923    return mesh
2924   
2925def square_outline(side_length = 1,up = "top", left = "left", right = "right",
2926                   down = "bottom", regions = False):
2927   
2928        a = Vertex (0,0)
2929        b = Vertex (0,side_length)
2930        c = Vertex (side_length,0)
2931        d = Vertex (side_length,side_length)
2932     
2933        s2 = Segment(b,d, tag = up)
2934        s3 = Segment(b,a, tag = left)
2935        s4 = Segment(d,c, tag = right)
2936        s5 = Segment(a,c, tag = down)
2937
2938        if regions:
2939            e =  Vertex (side_length/2,side_length/2)
2940            s6 = Segment(a,e, tag = down + left)
2941            s7 = Segment(b,e, tag = up + left)
2942            s8 = Segment(c,e, tag = down + right)
2943            s9 = Segment(d,e, tag = up + right)
2944            r1 = Region(side_length/2,3.*side_length/4, tag = up)
2945            r2 = Region(1.*side_length/4,side_length/2, tag = left)
2946            r3 = Region(3.*side_length/4,side_length/2, tag = right)
2947            r4 = Region(side_length/2,1.*side_length/4, tag = down)
2948            mesh = Mesh(userVertices=[a,b,c,d,e],
2949                        userSegments=[s2,s3,s4,s5,s6,s7,s8,s9],
2950                        regions = [r1,r2,r3,r4])
2951        else:
2952            mesh = Mesh(userVertices=[a,b,c,d],
2953                 userSegments=[s2,s3,s4,s5])
2954     
2955        return mesh
2956
2957   
2958
2959def region_strings2ints(region_list):
2960    """Given a list of (x_int,y_int,tag_string) lists it returns a list of
2961    (x_int,y_int,tag_int) and a list to convert the tag_int's back to
2962    the tag_strings 
2963    """
2964    # Make sure "" has an index of 0
2965    region_list.reverse()
2966    region_list.append((1.0,2.0,""))
2967    region_list.reverse()
2968    convertint2string = []
2969    for i in xrange(len(region_list)):
2970        convertint2string.append(region_list[i][2])
2971        if len(region_list[i]) == 4: # there's an area value
2972            region_list[i] = (region_list[i][0],
2973                              region_list[i][1],i,region_list[i][3])
2974        elif len(region_list[i]) == 3: # no area value
2975            region_list[i] = (region_list[i][0],region_list[i][1],i)
2976        else:
2977            print "The region list has a bad size"
2978            # raise an error ..
2979            raise Error
2980
2981    #remove "" from the region_list
2982    region_list.pop(0)
2983   
2984    return [region_list, convertint2string]
2985
2986def region_ints2strings(region_list,convertint2string):
2987    """Reverses the transformation of region_strings2ints
2988    """
2989    if region_list[0] != []:
2990        for i in xrange(len(region_list)):
2991            region_list[i] = [convertint2string[int(region_list[i][0])]]
2992    return region_list
2993
2994def segment_ints2strings(intlist, convertint2string):
2995    """Reverses the transformation of segment_strings2ints """
2996    stringlist = []
2997    for x in intlist:
2998        stringlist.append(convertint2string[x])
2999    return stringlist
3000
3001def segment_strings2ints(stringlist, preset):
3002    """Given a list of strings return a list of 0 to n ints which represent
3003    the strings and a converting list of the strings, indexed by 0 to n ints.
3004    Also, input an initial converting list of the strings
3005    Note, the converting list starts off with
3006    ["internal boundary", "external boundary", "internal boundary"]
3007    example input and output
3008    input ["a","b","a","c"],["c"]
3009    output [[2, 1, 2, 0], ['c', 'b', 'a']]
3010
3011    the first element in the converting list will be
3012    overwritten with "".
3013    ?This will become the third item in the converting list?
3014   
3015    # note, order the initial converting list is important,
3016     since the index = the int tag
3017    """
3018    nodups = unique(stringlist)
3019    # note, order is important, the index = the int tag
3020    #preset = ["internal boundary", "external boundary"]
3021    #Remove the preset tags from the list with no duplicates
3022    nodups = [x for x in nodups if x not in preset]
3023
3024    try:
3025        nodups.remove("") # this has to go to zero
3026    except ValueError:
3027        pass
3028   
3029    # Add the preset list at the beginning of no duplicates
3030    preset.reverse()
3031    nodups.extend(preset)
3032    nodups.reverse()
3033
3034       
3035    convertstring2int = {}
3036    convertint2string = []
3037    index = 0
3038    for x in nodups:
3039        convertstring2int[x] = index
3040        convertint2string.append(x)
3041        index += 1
3042    convertstring2int[""] = 0
3043
3044    intlist = []
3045    for x in stringlist:
3046        intlist.append(convertstring2int[x])
3047    return [intlist, convertint2string]
3048
3049
3050
3051
3052def unique(s):
3053    """Return a list of the elements in s, but without duplicates.
3054
3055    For example, unique([1,2,3,1,2,3]) is some permutation of [1,2,3],
3056    unique("abcabc") some permutation of ["a", "b", "c"], and
3057    unique(([1, 2], [2, 3], [1, 2])) some permutation of
3058    [[2, 3], [1, 2]].
3059
3060    For best speed, all sequence elements should be hashable.  Then
3061    unique() will usually work in linear time.
3062
3063    If not possible, the sequence elements should enjoy a total
3064    ordering, and if list(s).sort() doesn't raise TypeError it's
3065    assumed that they do enjoy a total ordering.  Then unique() will
3066    usually work in O(N*log2(N)) time.
3067
3068    If that's not possible either, the sequence elements must support
3069    equality-testing.  Then unique() will usually work in quadratic
3070    time.
3071    """
3072
3073    n = len(s)
3074    if n == 0:
3075        return []
3076
3077    # Try using a dict first, as that's the fastest and will usually
3078    # work.  If it doesn't work, it will usually fail quickly, so it
3079    # usually doesn't cost much to *try* it.  It requires that all the
3080    # sequence elements be hashable, and support equality comparison.
3081    u = {}
3082    try:
3083        for x in s:
3084            u[x] = 1
3085    except TypeError:
3086        del u  # move on to the next method
3087    else:
3088        return u.keys()
3089
3090    # We can't hash all the elements.  Second fastest is to sort,
3091    # which brings the equal elements together; then duplicates are
3092    # easy to weed out in a single pass.
3093    # NOTE:  Python's list.sort() was designed to be efficient in the
3094    # presence of many duplicate elements.  This isn't true of all
3095    # sort functions in all languages or libraries, so this approach
3096    # is more effective in Python than it may be elsewhere.
3097    try:
3098        t = list(s)
3099        t.sort()
3100    except TypeError:
3101        del t  # move on to the next method
3102    else:
3103        assert n > 0
3104        last = t[0]
3105        lasti = i = 1
3106        while i < n:
3107            if t[i] != last:
3108                t[lasti] = last = t[i]
3109                lasti += 1
3110            i += 1
3111        return t[:lasti]
3112
3113    # Brute force is all that's left.
3114    u = []
3115    for x in s:
3116        if x not in u:
3117            u.append(x)
3118    return u
3119
3120"""Refines triangles
3121
3122   Implements the #triangular bisection?# algorithm.
3123 
3124   
3125"""
3126
3127def Refine(mesh, triangles):
3128    """
3129    Given a general_mesh, and a triangle number, split
3130    that triangle in the mesh in half. Then to prevent
3131    vertices and edges from meeting, keep refining
3132    neighbouring triangles until the mesh is clean.
3133    """
3134    state = BisectionState(mesh)
3135    for triangle in triangles:
3136        if not state.refined_triangles.has_key(triangle):
3137            triangle.rotate_longest_side()
3138            state.start(triangle)
3139            Refine_mesh(mesh, state)
3140
3141def Refine_mesh(mesh, state):
3142    """
3143    """
3144    state.getState(mesh)
3145    refine_triangle(mesh,state)
3146    state.evolve()
3147    if not state.end:
3148        Refine_mesh(mesh,state)
3149
3150def refine_triangle(mesh,state):
3151    split(mesh,state.current_triangle,state.new_point)
3152    if state.case == 'one':
3153        state.r[3]=state.current_triangle#triangle 2
3154
3155        new_triangle_id = len(mesh.meshTriangles)-1
3156        new_triangle = mesh.meshTriangles[new_triangle_id]
3157
3158        split(mesh,new_triangle,state.old_point)
3159        state.r[2]=new_triangle#triangle 1.2
3160        state.r[4]=mesh.meshTriangles[len(mesh.meshTriangles)-1]#triangle 1.1
3161        r = state.r
3162        state.repairCaseOne()
3163
3164    if state.case == 'two':
3165        state.r[2]=mesh.meshTriangles[len(mesh.meshTriangles)-1]#triangle 1
3166
3167        new_triangle = state.current_triangle
3168
3169        split(mesh,new_triangle,state.old_point)
3170
3171        state.r[3]=mesh.meshTriangles[len(mesh.meshTriangles)-1]#triangle 2.1
3172        state.r[4]=new_triangle#triangle 2.2
3173        r = state.r
3174
3175        state.repairCaseTwo()
3176
3177    if state.case == 'vertex':
3178        state.r[2]=state.current_triangle#triangle 2
3179        state.r[3]=mesh.meshTriangles[len(mesh.meshTriangles)-1]#triangle 1
3180        r = state.r
3181        state.repairCaseVertex()
3182       
3183    if state.case == 'start':
3184        state.r[2]=mesh.meshTriangles[len(mesh.meshTriangles)-1]#triangle 1
3185        state.r[3]=state.current_triangle#triangle 2
3186
3187    if state.next_case == 'boundary':
3188        state.repairCaseBoundary()
3189
3190
3191def split(mesh, triangle, new_point):
3192        """
3193        Given a mesh, triangle_id and a new point,
3194        split the corrosponding triangle into two
3195        new triangles and update the mesh.
3196        """
3197
3198        new_triangle1 = Triangle(new_point,triangle.vertices[0],triangle.vertices[1],attribute = triangle.attribute, neighbors = None)
3199        new_triangle2 = Triangle(new_point,triangle.vertices[2],triangle.vertices[0],attribute = triangle.attribute, neighbors = None)
3200
3201        new_triangle1.setNeighbors(triangle.neighbors[2],None,new_triangle2)
3202        new_triangle2.setNeighbors(triangle.neighbors[1],new_triangle1,None)
3203
3204        mesh.meshTriangles.append(new_triangle1)
3205
3206        triangle.vertices = new_triangle2.vertices
3207        triangle.neighbors = new_triangle2.neighbors
3208
3209
3210class State:
3211
3212    def __init__(self):
3213        pass
3214
3215class BisectionState(State):
3216
3217
3218    def __init__(self,mesh):
3219        self.len = len(mesh.meshTriangles)
3220        self.refined_triangles = {}
3221        self.mesh = mesh
3222        self.current_triangle = None
3223        self.case = 'start'
3224        self.end = False
3225        self.r = [None,None,None,None,None]
3226
3227    def start(self, triangle):
3228        self.current_triangle = triangle
3229        self.case = 'start'
3230        self.end = False
3231        self.r = [None,None,None,None,None]
3232
3233    def getState(self,mesh):
3234        if not self.case == 'vertex':
3235            self.new_point=self.getNewVertex(mesh, self.current_triangle)
3236            #self.neighbour=self.getNeighbour(mesh, self.current_triangle)
3237            self.neighbour = self.current_triangle.neighbors[0]
3238            if not self.neighbour is None:
3239                self.neighbour.rotate_longest_side()
3240            self.next_case = self.get_next_case(mesh,self.neighbour,self.current_triangle)
3241        if self.case == 'vertex':
3242            self.new_point=self.old_point
3243
3244
3245    def evolve(self):
3246        if self.case == 'vertex':
3247            self.end = True
3248
3249        self.last_case = self.case
3250        self.case = self.next_case
3251
3252        self.old_point = self.new_point
3253        self.current_triangle = self.neighbour
3254
3255        if self.case == 'boundary':
3256            self.end = True
3257        self.refined_triangles[self.r[2]]=1
3258        self.refined_triangles[self.r[3]]=1
3259        if not self.r[4] is None:
3260            self.refined_triangles[self.r[4]]=1
3261        self.r[0]=self.r[2]
3262        self.r[1]=self.r[3]
3263
3264
3265    def getNewVertex(self,mesh,triangle):
3266        coordinate1 = triangle.vertices[1]
3267        coordinate2 = triangle.vertices[2]
3268        a = ([coordinate1.x*1.,coordinate1.y*1.])
3269        b = ([coordinate2.x*1.,coordinate2.y*1.])
3270        attributes = []
3271        for i in range(len(coordinate1.attributes)):
3272            att = (coordinate1.attributes[i]+coordinate2.attributes[i])/2
3273            attributes.append(att)
3274        new_coordinate = [((a[0]-b[0])/2+b[0]),((a[1]-b[1])/2+b[1])]
3275        newVertex = Vertex(new_coordinate[0],new_coordinate[1], attributes = attributes)
3276        mesh.maxVertexIndex+=1
3277        newVertex.index = mesh.maxVertexIndex
3278        mesh.meshVertices.append(newVertex)
3279        return newVertex
3280
3281    def get_next_case(self, mesh,neighbour,triangle):
3282        """
3283        Given the locations of two neighbouring triangles,
3284        examine the interior indices of their vertices (i.e.
3285        0,1 or 2) to determine what how the neighbour needs
3286        to be refined.
3287        """
3288        if (neighbour is None):
3289                next_case = 'boundary'
3290        else:
3291                if triangle.vertices[1].x==neighbour.vertices[2].x:
3292                    if triangle.vertices[1].y==neighbour.vertices[2].y:
3293                        next_case = 'vertex'
3294                if triangle.vertices[1].x==neighbour.vertices[0].x:
3295                    if triangle.vertices[1].y==neighbour.vertices[0].y:
3296                        next_case = 'two'
3297                if triangle.vertices[1].x==neighbour.vertices[1].x:
3298                    if triangle.vertices[1].y==neighbour.vertices[1].y:
3299                        next_case = 'one'
3300        return next_case
3301
3302
3303
3304    def repairCaseVertex(self):
3305
3306        r = self.r
3307
3308
3309        self.link(r[0],r[2])
3310        self.repair(r[0])
3311
3312        self.link(r[1],r[3])
3313        self.repair(r[1])
3314
3315        self.repair(r[2])
3316
3317        self.repair(r[3])
3318
3319
3320    def repairCaseOne(self):
3321        r = self.rkey
3322
3323
3324        self.link(r[0],r[2])
3325        self.repair(r[0])
3326
3327        self.link(r[1],r[4])
3328        self.repair(r[1])
3329
3330        self.repair(r[4])
3331
3332    def repairCaseTwo(self):
3333        r = self.r
3334
3335        self.link(r[0],r[4])
3336        self.repair(r[0])
3337
3338        self.link(r[1],r[3])
3339        self.repair(r[1])
3340
3341        self.repair(r[4])
3342
3343    def repairCaseBoundary(self):
3344        r = self.r
3345        self.repair(r[2])
3346        self.repair(r[3])
3347
3348
3349
3350    def repair(self,triangle):
3351        """
3352        Given a triangle that knows its neighbours, this will
3353        force the neighbours to comply.
3354
3355        However, it needs to compare the vertices of triangles
3356        for this implementation
3357
3358        But it doesn't work for invalid neighbour structures
3359        """
3360        n=triangle.neighbors
3361        for i in (0,1,2):
3362            if not n[i] is None:
3363                for j in (0,1,2):#to find which side of the list is broken
3364                    if not (n[i].vertices[j] in triangle.vertices):
3365                    #ie if j is the side of n that needs fixing
3366                        k = j
3367                n[i].neighbors[k]=triangle
3368
3369
3370
3371    def link(self,triangle1,triangle2):
3372        """
3373        make triangle1 neighbors point to t
3374                #count = 0riangle2
3375        """
3376        count = 0
3377        for i in (0,1,2):#to find which side of the list is broken
3378            if not (triangle1.vertices[i] in triangle2.vertices):
3379                j = i
3380                count+=1
3381        assert count == 1
3382        triangle1.neighbors[j]=triangle2
3383
3384class Discretised_Tuple_Set:
3385    """
3386    if a={(0.0):[(0.01),(0.02)],(0.2):[(0.17)]}
3387    a[(0.01)]=a[(0.0)]=[(0.01),(0.02)]
3388    a[(10000)]=[] #NOT KEYERROR
3389
3390    a.append[(0.01)]
3391    => {0.0:[(0.01),(0.02),(0.01)],0.2:[(0.17)]}
3392
3393    #NOT IMPLIMENTED
3394    a.remove[(0.01)]
3395    => {(0.0):[(0.02),(0.01)],0.2:[(0.17)]}
3396
3397    a.remove[(0.17)]
3398    => {(0.0):[(0.02),(0.01)],0.2:[]}
3399    #NOT IMPLIMENTED
3400    at a.precision = 2:
3401        a.round_up_rel[0.0]=
3402        a.round_flat[0.0]=
3403        a.round_down_rel[0.0]=
3404
3405        a.up((0.1,2.04))=
3406
3407    If t_rel = 0, nothing gets rounded into
3408    two bins. If t_rel = 0.5, everything does.
3409
3410    Ideally, precision can be set high, so that multiple
3411    entries are rarely in the same bin. And t_rel should
3412    be low (<0.1 for 1 dimension!,<(0.1/n) for small n!!)
3413    so that it is rare to put items in mutiple bins.
3414
3415    Ex bins per entry = product(a,b,c...,n)
3416    a = 1 or 2 s.t. Ex(a) = 1+2*t_rel
3417    b = 1 or 2 ...
3418
3419    BUT!!! to avoid missing the right bin:
3420    (-10)**(precision+1)*t_rel must be greater than the
3421    greatest possible variation that an identical element
3422    can display.
3423
3424
3425    Note that if tol = 0.5 (the max allowed) 0.6 will round to .7 and .5
3426    but not .6 - this looks wrong, but note that *everything* will round,
3427    so .6 wont be missed as everything close to it will check in .7 and .5.
3428    """
3429    def __init__(self,p_rel = 6,t_rel = 0.01):
3430        self.__p_rel__ = p_rel
3431        self.__t_rel__ = t_rel
3432
3433        self.__p_abs__ = p_rel+1
3434        self.__t_abs__ = t_rel
3435
3436        assert t_rel <= 0.5
3437        self.__items__ = {}
3438        from math import frexp
3439        self.frexp = frexp
3440        roundings = [self.round_up_rel,\
3441        self.round_down_rel,self.round_flat_rel,\
3442        self.round_down_abs,self.round_up_abs,\
3443        self.round_flat_abs]#
3444
3445        self.roundings = roundings
3446
3447    def __repr__(self):
3448        return '%s'%self.__items__
3449
3450    def rounded_keys(self,key):
3451        key = tuple(key)
3452        keys = [key]
3453        keys = self.__rounded_keys__(key)
3454        return (keys)
3455
3456    def __rounded_keys__(self,key):
3457        keys = []
3458        rounded_key=list(key)
3459        rounded_values=list(key)
3460
3461        roundings = list(self.roundings)
3462
3463        #initialise rounded_values
3464        round = roundings.pop(0)
3465        for i in range(len(rounded_values)):
3466            rounded_key[i]=round(key[i])
3467            rounded_values[i]={}
3468            rounded_values[i][rounded_key[i]]=None
3469        keys.append(tuple(rounded_key))
3470       
3471        for round in roundings:
3472            for i in range(len(rounded_key)):
3473                rounded_value=round(key[i])
3474                if not rounded_values[i].has_key(rounded_value):
3475                    #ie unless round_up_rel = round_down_rel
3476                    #so the keys stay unique
3477                    for j in range(len(keys)):
3478                        rounded_key = list(keys[j])
3479                        rounded_key[i]=rounded_value
3480                        keys.append(tuple(rounded_key))
3481        return keys
3482
3483    def append(self,item):
3484        keys = self.rounded_keys(item)
3485        for key in keys:
3486            if self.__items__.has_key(key):
3487                self.__items__[key].append(item)
3488            else:
3489                self.__items__[key]=[item]
3490
3491    def __getitem__(self,key):
3492        answer = []
3493        keys = self.rounded_keys(key)
3494        for key in keys:
3495            if self.__items__.has_key(key):
3496                answer.extend(self.__items__[key])
3497        #if len(answer)==0:
3498        #    raise KeyError#FIXME or return KeyError
3499        #                  #FIXME or just return []?
3500        else:
3501            return answer #FIXME or unique(answer)?
3502
3503    def __delete__(self,item):
3504        keys = self.rounded_keys(item)
3505        answer = False
3506        #if any of the possible keys contains
3507        #a list, return true
3508        for key in keys:       
3509            if self.__items__.has_key(key):
3510                if item in self.__items__[key]:
3511                    self.__items__[key].remove(item)
3512
3513    def remove(self,item):
3514        self.__delete__(item)
3515
3516    def __contains__(self,item):
3517
3518        keys = self.rounded_keys(item)
3519        answer = False
3520        #if any of the possible keys contains
3521        #a list, return true
3522        for key in keys:       
3523            if self.__items__.has_key(key):
3524                if item in self.__items__[key]:
3525                    answer = True
3526        return answer
3527
3528
3529    def has_item(self,item):
3530        return self.__contains__(item)
3531
3532    def round_up_rel2(self,value):
3533         t_rel=self.__t_rel__
3534         #Rounding up the value
3535         m,e = self.frexp(value)
3536         m = m/2
3537         e = e + 1
3538         #m is the mantissa, e the exponent
3539         # 0.5 < |m| < 1.0
3540         m = m+t_rel*(10**-(self.__p_rel__))
3541         #bump m up
3542         m = round(m,self.__p_rel__)
3543         return m*(2.**e)
3544
3545    def round_down_rel2(self,value):
3546         t_rel=self.__t_rel__
3547         #Rounding down the value
3548         m,e = self.frexp(value)
3549         m = m/2
3550         e = e + 1
3551         #m is the mantissa, e the exponent
3552         # 0.5 < m < 1.0
3553         m = m-t_rel*(10**-(self.__p_rel__))
3554         #bump the |m| down, by 5% or whatever
3555         #self.p_rel dictates
3556         m = round(m,self.__p_rel__)
3557         return m*(2.**e)
3558
3559    def round_flat_rel2(self,value):
3560    #redundant
3561         m,e = self.frexp(value)
3562         m = m/2
3563         e = e + 1
3564         m = round(m,self.__p_rel__)
3565         return m*(2.**e)
3566
3567    def round_up_rel(self,value):
3568         t_rel=self.__t_rel__
3569         #Rounding up the value
3570         m,e = self.frexp(value)
3571         #m is the mantissa, e the exponent
3572         # 0.5 < |m| < 1.0
3573         m = m+t_rel*(10**-(self.__p_rel__))
3574         #bump m up
3575         m = round(m,self.__p_rel__)
3576         return m*(2.**e)
3577
3578    def round_down_rel(self,value):
3579         t_rel=self.__t_rel__
3580         #Rounding down the value
3581         m,e = self.frexp(value)
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_rel(self,value):
3591    #redundant
3592         m,e = self.frexp(value)
3593         m = round(m,self.__p_rel__)
3594         return m*(2.**e)
3595
3596    def round_up_abs(self,value):
3597         t_abs=self.__t_abs__
3598         #Rounding up the value
3599         m = value+t_abs*(10**-(self.__p_abs__))
3600         #bump m up
3601         m = round(m,self.__p_abs__)
3602         return m
3603
3604    def round_down_abs(self,value):
3605         t_abs=self.__t_abs__
3606         #Rounding down the value
3607         m = value-t_abs*(10**-(self.__p_abs__))
3608         #bump the |m| down, by 5% or whatever
3609         #self.p_rel dictates
3610         m = round(m,self.__p_abs__)
3611         return m
3612
3613    def round_flat_abs(self,value):
3614    #redundant?
3615         m = round(value,self.__p_abs__)
3616         return m
3617
3618    def keys(self):
3619        return self.__items__.keys()
3620
3621
3622class Mapped_Discretised_Tuple_Set(Discretised_Tuple_Set):
3623    """
3624    This is a discretised tuple set, but
3625    based on a mapping. The mapping MUST
3626    return a sequence.
3627
3628    example:
3629    def weight(animal):
3630        return [animal.weight]
3631   
3632    a = Mapped_Discretised_Tuple_Set(weight)
3633    a.append[cow]
3634    a.append[fox]
3635    a.append[horse]
3636
3637    a[horse]    -> [cow,horse]
3638    a[dog]      -> [fox]
3639    a[elephant] -> []
3640    """
3641    def __init__(self,mapping,p_rel = 6, t_rel=0.01):
3642        Discretised_Tuple_Set.__init__\
3643         (self,p_rel,t_rel = t_rel)
3644        self.mapping = mapping
3645
3646    def rounded_keys(self,key):
3647        mapped_key = tuple(self.mapping(key))
3648        keys = self.__rounded_keys__(mapped_key)
3649        return keys
3650
3651class Affine_Linespace(Mapped_Discretised_Tuple_Set):
3652    """
3653    The affine linespace creates a way to record and compare lines.
3654    Precision is a bit of a hack, but it creates a way to avoid
3655    misses caused by round offs (between two lines of different
3656    lenghts, the short one gets rounded off more).
3657    I am starting to think that a quadratic search would be faster.
3658    Nearly.
3659    """
3660    def __init__(self,p_rel=4,t_rel=0.2):
3661        Mapped_Discretised_Tuple_Set.__init__\
3662            (self,self.affine_line,\
3663            p_rel=p_rel,t_rel=t_rel)
3664
3665        roundings = \
3666        [self.round_down_rel,self.round_up_rel,self.round_flat_rel]
3667        self.roundings = roundings
3668        #roundings = \
3669        #[self.round_down_abs,self.round_up_abs,self.round_flat_abs]
3670        #self.roundings = roundings
3671
3672    def affine_line(self,line):
3673        point_1 = line[0]
3674        point_2 = line[1]
3675        #returns the equation of a line
3676        #between two points, in the from
3677        #(a,b,-c), as in ax+by-c=0
3678        #or line *dot* (x,y,1) = (0,0,0)
3679
3680        #Note that it normalises the line
3681        #(a,b,-c) so Line*Line = 1.
3682        #This does not change the mathematical
3683        #properties, but it makes comparism
3684        #easier.
3685
3686        #There are probably better algorithms.
3687        x1 = point_1[0]
3688        x2 = point_2[0]
3689        y1 = point_1[1]
3690        y2 = point_2[1]
3691        dif_x = x1-x2
3692        dif_y = y1-y2
3693
3694        if dif_x == dif_y == 0:
3695            msg = 'points are the same'
3696            raise msg
3697        elif abs(dif_x)>=abs(dif_y):
3698            alpha = (-dif_y)/dif_x
3699            #a = alpha * b
3700            b = -1.
3701            c = (x1*alpha+x2*alpha+y1+y2)/2.
3702            a = alpha*b
3703        else:
3704            beta = dif_x/(-dif_y)
3705            #b = beta * a
3706            a = 1.
3707            c = (x1+x2+y1*beta+y2*beta)/2.
3708            b = beta*a
3709        mag = abs(a)+abs(b)
3710        #This does not change the mathematical
3711        #properties, but it makes comparism possible.
3712
3713        #note that the gradient is b/a, or (a/b)**-1.
3714        #so
3715
3716        #if a == 0:
3717        #    sign_a = 1.
3718        #else:
3719        #    sign_a = a/((a**2)**0.5)
3720        #if b == 0:
3721        #    sign_b = 1.
3722        #else:
3723        #    sign_b = b/((b**2)**0.5)
3724        #if c == 0:
3725        #    sign_c = 1.
3726        #else:
3727        #    sign_c = c/((c**2)**0.5)
3728        #a = a/mag*sign_a
3729        #b = b/mag*sign_b
3730        #c = c/mag*sign_c
3731        a = a/mag
3732        b = b/mag
3733        c = c/mag
3734        return a,b,c
3735
3736
3737# FIXME (DSG-DSG)
3738# To do-
3739# Create a clear interface. eg
3740# have the interface methods more at the top of this file and add comments
3741# for the interface functions/methods, use function_name (not functionName),
3742
3743#Currently
3744#function_name methods assume absolute values.  Geo-refs can be passed in.
3745#
3746
3747# instead of functionName
3748if __name__ == "__main__":
3749    #from mesh import *
3750    # THIS CAN BE DELETED
3751    m = Mesh()
3752    dict = importUngenerateFile("ungen_test.txt")
3753    m.addVertsSegs(dict)
3754    print m3
Note: See TracBrowser for help on using the repository browser.