# source:inundation/pyvolution/util.py@1859

Last change on this file since 1859 was 1859, checked in by ole, 18 years ago

File size: 26.0 KB
Line
1"""This module contains various auxiliary function used by pyvolution.
2
3It is also a clearing house for functions that may later earn a module
4of their own.
5"""
6
7#FIXME: Move polygon stuff out
8
9def angle(v):
10    """Compute angle between e1 (the unit vector in the x-direction)
11    and the specified vector
12    """
13
14    from math import acos, pi, sqrt
15    from Numeric import sum, array
16
17    l = sqrt( sum (array(v)**2))
18    v1 = v[0]/l
19    v2 = v[1]/l
20
21
22    theta = acos(v1)
23
24    #try:
25    #   theta = acos(v1)
26    #except ValueError, e:
27    #    print 'WARNING (util.py): Angle acos(%s) failed: %s' %(str(v1), e)
28    #
29    #    #FIXME, hack to avoid acos(1.0) Value error
30    #    # why is it happening?
31    #    # is it catching something we should avoid?
32    #    #FIXME (Ole): When did this happen? We need a unit test to
33    #    #reveal this condition or otherwise remove the hack
34    #
35    #    s = 1e-6
36    #    if (v1+s >  1.0) and (v1-s < 1.0) :
37    #        theta = 0.0
38    #    elif (v1+s >  -1.0) and (v1-s < -1.0):
39    #        theta = 3.1415926535897931
40    #    print 'WARNING (util.py): angle v1 is %f, setting acos to %f '\
41    #          %(v1, theta)
42
43    if v2 < 0:
44        #Quadrant 3 or 4
45        theta = 2*pi-theta
46
47    return theta
48
49
50def anglediff(v0, v1):
51    """Compute difference between angle of vector x0, y0 and x1, y1.
52    This is used for determining the ordering of vertices,
53    e.g. for checking if they are counter clockwise.
54
55    Always return a positive value
56    """
57
58    from math import pi
59
60    a0 = angle(v0)
61    a1 = angle(v1)
62
63    #Ensure that difference will be positive
64    if a0 < a1:
65        a0 += 2*pi
66
67    return a0-a1
68
69
70def mean(x):
71    from Numeric import sum
72    return sum(x)/len(x)
73
74
75def point_on_line(x, y, x0, y0, x1, y1):
76    """Determine whether a point is on a line segment
77
78    Input: x, y, x0, x0, x1, y1: where
79        point is given by x, y
80        line is given by (x0, y0) and (x1, y1)
81
82    """
83
84    from Numeric import array, dot, allclose
85    from math import sqrt
86
87    a = array([x - x0, y - y0])
88    a_normal = array([a[1], -a[0]])
89
90    b = array([x1 - x0, y1 - y0])
91
92    if dot(a_normal, b) == 0:
93        #Point is somewhere on the infinite extension of the line
94
95        len_a = sqrt(sum(a**2))
96        len_b = sqrt(sum(b**2))
97        if dot(a, b) >= 0 and len_a <= len_b:
98           return True
99        else:
100           return False
101    else:
102      return False
103
104def ensure_numeric(A, typecode = None):
105    """Ensure that sequence is a Numeric array.
106    Inputs:
107        A: Sequance. If A is already a Numeric array it will be returned
108                     unaltered
109                     If not, an attempt is made to convert it to a Numeric
110                     array
111        typecode: Numeric type. If specified, use this in the conversion.
112                                If not, let Numeric decide
113
114    This function is necessary as array(A) can cause memory overflow.
115    """
116
117    from Numeric import ArrayType, array
118
119    if typecode is None:
120        if type(A) == ArrayType:
121            return A
122        else:
123            return array(A)
124    else:
125        if type(A) == ArrayType:
126            if A.typecode == typecode:
127                return array(A)  #FIXME: Shouldn't this be just return A?
128            else:
129                return A.astype(typecode)
130        else:
131            return array(A).astype(typecode)
132
133
134
135def file_function(filename,
136                  domain = None,
137                  quantities = None,
138                  interpolation_points = None, verbose = False):
139    """Read time history of spatial data from NetCDF file and return
140    a callable object.
141
142    If the file has extension 'sww' then it is assumed to be spatio-temporal
143    or temporal and the callable object will have the form f(t,x,y) or f(t)
144    depending on whether the file contains spatial data
145
146    If the file has extension 'tms' then it is assumed to be temporal only
147    and the callable object will have the form f(t)
148
149    Either form will return interpolated values based on the input file
150    using the underlying interpolation_function.
151
152    If domain is specified, model time (domain.starttime)
153    will be checked and possibly modified.
154
155    All times are assumed to be in UTC
156
157    All spatial information is assumed to be in UTM coordinates.
158
159    See Interpolation function for further documentation
160    """
161
162
163    #FIXME (OLE): Should check origin of domain against that of file
164    #In fact, this is where origin should be converted to that of domain
165    #Also, check that file covers domain fully.
166    #If we use the suggested Point_set class for interpolation points
167    #here it would be easier
168
169    #Take into account:
170    #- domain's georef
171    #- sww file's georef
172    #- interpolation points georef
173
174
175    assert type(filename) == type(''),\
176               'First argument to File_function must be a string'
177
178    try:
179        fid = open(filename)
180    except Exception, e:
181        msg = 'File "%s" could not be opened: Error="%s"'\
182                  %(filename, e)
183        raise msg
184
185    line = fid.readline()
186    fid.close()
187
188    if quantities is None:
189        if domain is not None:
190            quantities = domain.conserved_quantities
191
192
193
194    if line[:3] == 'CDF':
195        return get_netcdf_file_function(filename, domain, quantities,
196                                        interpolation_points, verbose = verbose)
197    else:
198        raise 'Must be a NetCDF File'
199
200
201
202def get_netcdf_file_function(filename, domain=None, quantity_names=None,
203                             interpolation_points=None, verbose = False):
204    """Read time history of spatial data from NetCDF sww file and
205    return a callable object f(t,x,y)
206    which will return interpolated values based on the input file.
207
208    If domain is specified, model time (domain.starttime)
209    will be checked and possibly modified
210
211    All times are assumed to be in UTC
212
213    See Interpolation function for further documetation
214
215
216    """
217
218    #verbose = True
219
220    #FIXME: Check that model origin is the same as file's origin
221    #(both in UTM coordinates)
222    #If not - modify those from file to match domain
223    #Take this code from e.g. dem2pts in data_manager.py
224    #FIXME: Use geo_reference to read and write xllcorner...
225
226
227    import time, calendar, types
228    from config import time_format
229    from Scientific.IO.NetCDF import NetCDFFile
230    from Numeric import array, zeros, Float, alltrue, concatenate, reshape
231
232    #Open NetCDF file
233    if verbose: print 'Reading', filename
234    fid = NetCDFFile(filename, 'r')
235
236    if type(quantity_names) == types.StringType:
237        quantity_names = [quantity_names]
238
239    if quantity_names is None or len(quantity_names) < 1:
240        #If no quantities are specified get quantities from file
241        #x, y, time are assumed as the independent variables so
242        #they are excluded as potentiol quantities
243        quantity_names = []
244        for name in fid.variables.keys():
245            if name not in ['x', 'y', 'time']:
246                quantity_names.append(name)
247
248    if len(quantity_names) < 1:
249        msg = 'ERROR: At least one independent value must be specified'
250        raise msg
251
252
253    #Now assert that requested quantitites (and the independent ones)
254    #are present in file
255    missing = []
256    for quantity in ['time'] + quantity_names:
257        if not fid.variables.has_key(quantity):
258            missing.append(quantity)
259
260    if len(missing) > 0:
261        msg = 'Quantities %s could not be found in file %s'\
262              %(str(missing), filename)
263        fid.close()
264        raise msg
265
266    #Decide whether this data has a spatial dimension
267    spatial = True
268    for quantity in ['x', 'y']:
269        if not fid.variables.has_key(quantity):
270            spatial = False
271
272    if filename[-3:] == 'tms' and spatial is True:
273        msg = 'Files of type tms must not contain spatial information'
274        raise msg
275
276    if filename[-3:] == 'sww' and spatial is False:
277        msg = 'Files of type sww must contain spatial information'
278        raise msg
279
280    #Get first timestep
281    try:
282        starttime = fid.starttime[0]
283    except ValueError:
284        msg = 'Could not read starttime from file %s' %filename
285        raise msg
286
287    #Get variables
288    if verbose: print 'Get variables'
289    time = fid.variables['time'][:]
290
291    if spatial:
292        #Get origin
293        xllcorner = fid.xllcorner[0]
294        yllcorner = fid.yllcorner[0]
295        zone = fid.zone[0]
296
297        x = fid.variables['x'][:]
298        y = fid.variables['y'][:]
299        triangles = fid.variables['volumes'][:]
300
301        x = reshape(x, (len(x),1))
302        y = reshape(y, (len(y),1))
303        vertex_coordinates = concatenate((x,y), axis=1) #m x 2 array
304
305
306    if domain is not None:
307        #Update domain.startime if it is *earlier* than starttime
308        if starttime > domain.starttime:
309            msg = 'WARNING: Start time as specified in domain (%f)'\
310                  %domain.starttime
311            msg += ' is earlier than the starttime of file %s (%f).'\
312                     %(filename, starttime)
313            msg += ' Modifying domain starttime accordingly.'
314
315            if verbose: print msg
316            domain.starttime = starttime #Modifying model time
317            if verbose: print 'Domain starttime is now set to %f'\
318               %domain.starttime
319
320
321        #If domain.startime is *later* than starttime,
322        #move time back - relative to domain's time
323        if domain.starttime > starttime:
324            time = time - domain.starttime + starttime
325
326        #FIXME Use method in geo to reconcile
327        #if spatial:
328        #assert domain.geo_reference.xllcorner == xllcorner
329        #assert domain.geo_reference.yllcorner == yllcorner
330        #assert domain.geo_reference.zone == zone
331
332    if verbose:
333        print 'File_function data obtained from: %s' %filename
334        print '  References:'
335        #print '    Datum ....' #FIXME
336        if spatial:
337            print '    Lower left corner: [%f, %f]'\
338                  %(xllcorner, yllcorner)
339        print '    Start time:   %f' %starttime
340
341
342    #Produce values for desired data points at
343    #each timestep
344
345    quantities = {}
346    for i, name in enumerate(quantity_names):
347        quantities[name] = fid.variables[name][:]
348    fid.close()
349
350
351    from least_squares import Interpolation_function
352
353    if not spatial:
354        vertex_coordinates = triangles = interpolation_points = None
355
356    return Interpolation_function(time,
357                                  quantities,
358                                  quantity_names,
359                                  vertex_coordinates,
360                                  triangles,
361                                  interpolation_points,
362                                  verbose = verbose)
363
364
365
366
367
368def read_xya(filename, format = 'netcdf'):
369    """Read simple xya file, possibly with a header in the first line, just
370    x y [attributes]
371    separated by whitespace
372
373    Format can be either ASCII or NetCDF
374
375    Return list of points, list of attributes and
376    attribute names if present otherwise None
377    """
378    #FIXME: Probably obsoleted by similar function in load_ASCII
379    #FIXME: Name read_data_points (or see 'load' in pylab)
380
381    from Scientific.IO.NetCDF import NetCDFFile
382
383    if format.lower() == 'netcdf':
384        #Get NetCDF
385
386        fid = NetCDFFile(filename, 'r')
387
388        # Get the variables
389        points = fid.variables['points'][:]
390        keys = fid.variables.keys()
391        attributes = {}
392        for key in keys:
393            if key != 'points':
394                attributes[key] = fid.variables[key][:]
395
396        fid.close()
397    else:
398        #Read as ASCII file assuming that it is separated by whitespace
399        fid = open(filename)
400        lines = fid.readlines()
401        fid.close()
402
403        #Check if there is a header line
404        fields = lines[0].strip().split()
405        try:
406            float(fields[0])
407        except:
408            #This must be a header line
409            attribute_names = fields
410            lines = lines[1:]
411        else:
412            attribute_names = ['elevation']  #HACK
413
414        attributes = {}
415        for key in attribute_names:
416            attributes[key] = []
417
418        points = []
419        for line in lines:
420            fields = line.strip().split()
421            points.append( (float(fields[0]), float(fields[1])) )
422            for i, key in enumerate(attribute_names):
423                attributes[key].append( float(fields[2+i]) )
424
425    return points, attributes
426
427
428#####################################
429#POLYGON STUFF
430#
431#FIXME: All these should be put into new module polygon.py
432
433
434
435#These have been redefined to use separate_by_polygon which will
436#put all inside indices in the first part of the indices array and
437#outside indices in the last
438
439def inside_polygon(points, polygon, closed = True, verbose = False):
440    """See separate_points_by_polygon for documentation
441    """
442
443    from Numeric import array, Float, reshape
444
445    if verbose: print 'Checking input to inside_polygon'
446    try:
447        points = ensure_numeric(points, Float)
448    except:
449        msg = 'Points could not be converted to Numeric array'
450        raise msg
451
452    try:
453        polygon = ensure_numeric(polygon, Float)
454    except:
455        msg = 'Polygon could not be converted to Numeric array'
456        raise msg
457
458
459
460    if len(points.shape) == 1:
461        one_point = True
462        points = reshape(points, (1,2))
463    else:
464        one_point = False
465
466    indices, count = separate_points_by_polygon(points, polygon,
467                                                closed, verbose)
468
469    if one_point:
470        return count == 1
471    else:
472        return indices[:count]
473
474def outside_polygon(points, polygon, closed = True, verbose = False):
475    """See separate_points_by_polygon for documentation
476    """
477
478    from Numeric import array, Float, reshape
479
480    if verbose: print 'Checking input to outside_polygon'
481    try:
482        points = ensure_numeric(points, Float)
483    except:
484        msg = 'Points could not be converted to Numeric array'
485        raise msg
486
487    try:
488        polygon = ensure_numeric(polygon, Float)
489    except:
490        msg = 'Polygon could not be converted to Numeric array'
491        raise msg
492
493
494
495    if len(points.shape) == 1:
496        one_point = True
497        points = reshape(points, (1,2))
498    else:
499        one_point = False
500
501    indices, count = separate_points_by_polygon(points, polygon,
502                                                closed, verbose)
503
504    if one_point:
505        return count != 1
506    else:
507        return indices[count:][::-1]  #return reversed
508
509
510def separate_points_by_polygon(points, polygon,
511                               closed = True, verbose = False):
512    """Determine whether points are inside or outside a polygon
513
514    Input:
515       points - Tuple of (x, y) coordinates, or list of tuples
516       polygon - list of vertices of polygon
517       closed - (optional) determine whether points on boundary should be
518       regarded as belonging to the polygon (closed = True)
519       or not (closed = False)
520
521    Outputs:
522       indices: array of same length as points with indices of points falling
523       inside the polygon listed from the beginning and indices of points
524       falling outside listed from the end.
525
526       count: count of points falling inside the polygon
527
528       The indices of points inside are obtained as indices[:count]
529       The indices of points outside are obtained as indices[count:]
530
531
532    Examples:
533       U = [[0,0], [1,0], [1,1], [0,1]] #Unit square
534
535       separate_points_by_polygon( [[0.5, 0.5], [1, -0.5], [0.3, 0.2]], U)
536       will return the indices [0, 2, 1] and count == 2 as only the first
537       and the last point are inside the unit square
538
539    Remarks:
540       The vertices may be listed clockwise or counterclockwise and
541       the first point may optionally be repeated.
542       Polygons do not need to be convex.
543       Polygons can have holes in them and points inside a hole is
544       regarded as being outside the polygon.
545
546    Algorithm is based on work by Darel Finley,
547    http://www.alienryderflex.com/polygon/
548
549    """
550
551    from Numeric import array, Float, reshape, Int, zeros
552
553
554    #Input checks
555    try:
556        points = ensure_numeric(points, Float)
557    except:
558        msg = 'Points could not be converted to Numeric array'
559        raise msg
560
561    try:
562        polygon = ensure_numeric(polygon, Float)
563    except:
564        msg = 'Polygon could not be converted to Numeric array'
565        raise msg
566
567    assert len(polygon.shape) == 2,\
568       'Polygon array must be a 2d array of vertices'
569
570    assert polygon.shape[1] == 2,\
571       'Polygon array must have two columns'
572
573    assert len(points.shape) == 2,\
574       'Points array must be a 2d array'
575
576    assert points.shape[1] == 2,\
577       'Points array must have two columns'
578
579    N = polygon.shape[0] #Number of vertices in polygon
580    M = points.shape[0]  #Number of points
581
582    px = polygon[:,0]
583    py = polygon[:,1]
584
585    #Used for an optimisation when points are far away from polygon
586    minpx = min(px); maxpx = max(px)
587    minpy = min(py); maxpy = max(py)
588
589
590    #Begin main loop
591    indices = zeros(M, Int)
592
593    inside_index = 0    #Keep track of points inside
594    outside_index = M-1 #Keep track of points outside (starting from end)
595
596    for k in range(M):
597
598        if verbose:
599            if k %((M+10)/10)==0: print 'Doing %d of %d' %(k, M)
600
601        #for each point
602        x = points[k, 0]
603        y = points[k, 1]
604
605        inside = False
606
607        if not x > maxpx or x < minpx or y > maxpy or y < minpy:
608            #Check polygon
609            for i in range(N):
610                j = (i+1)%N
611
612                #Check for case where point is contained in line segment
613                if point_on_line(x, y, px[i], py[i], px[j], py[j]):
614                    if closed:
615                        inside = True
616                    else:
617                        inside = False
618                    break
619                else:
620                    #Check if truly inside polygon
621                    if py[i] < y and py[j] >= y or\
622                       py[j] < y and py[i] >= y:
623                        if px[i] + (y-py[i])/(py[j]-py[i])*(px[j]-px[i]) < x:
624                            inside = not inside
625
626        if inside:
627            indices[inside_index] = k
628            inside_index += 1
629        else:
630            indices[outside_index] = k
631            outside_index -= 1
632
633    return indices, inside_index
634
635
636def separate_points_by_polygon_c(points, polygon,
637                                 closed = True, verbose = False):
638    """Determine whether points are inside or outside a polygon
639
640    C-wrapper
641    """
642
643    from Numeric import array, Float, reshape, zeros, Int
644
645
646    if verbose: print 'Checking input to separate_points_by_polygon'
647    #Input checks
648    try:
649        points = ensure_numeric(points, Float)
650    except:
651        msg = 'Points could not be converted to Numeric array'
652        raise msg
653
654    #if verbose: print 'Checking input to separate_points_by_polygon 2'
655    try:
656        polygon = ensure_numeric(polygon, Float)
657    except:
658        msg = 'Polygon could not be converted to Numeric array'
659        raise msg
660
661    if verbose: print 'check'
662
663    assert len(polygon.shape) == 2,\
664       'Polygon array must be a 2d array of vertices'
665
666    assert polygon.shape[1] == 2,\
667       'Polygon array must have two columns'
668
669    assert len(points.shape) == 2,\
670       'Points array must be a 2d array'
671
672    assert points.shape[1] == 2,\
673       'Points array must have two columns'
674
675    N = polygon.shape[0] #Number of vertices in polygon
676    M = points.shape[0]  #Number of points
677
678    from util_ext import separate_points_by_polygon
679
680    if verbose: print 'Allocating array for indices'
681
682    indices = zeros( M, Int )
683
684    if verbose: print 'Calling C-version of inside poly'
685    count = separate_points_by_polygon(points, polygon, indices,
686                                       int(closed), int(verbose))
687
688    return indices, count
689
690
691
692class Polygon_function:
693    """Create callable object f: x,y -> z, where a,y,z are vectors and
694    where f will return different values depending on whether x,y belongs
695    to specified polygons.
696
697    To instantiate:
698
699       Polygon_function(polygons)
700
701    where polygons is a list of tuples of the form
702
703      [ (P0, v0), (P1, v1), ...]
704
705      with Pi being lists of vertices defining polygons and vi either
706      constants or functions of x,y to be applied to points with the polygon.
707
708    The function takes an optional argument, default which is the value
709    (or function) to used for points not belonging to any polygon.
710    For example:
711
712       Polygon_function(polygons, default = 0.03)
713
714    If omitted the default value will be 0.0
715
716    Note: If two polygons overlap, the one last in the list takes precedence
717
718    FIXME? : Currently, coordinates specified here are assumed to be relative to the origin (georeference) used by domain. This function is more general than domain so perhaps it'd be an idea to allow the optional argument georeference???
719
720    """
721
722    def __init__(self, regions, default = 0.0):
723
724        try:
725            len(regions)
726        except:
727            msg = 'Polygon_function takes a list of pairs (polygon, value). Got %s' %polygons
728            raise msg
729
730
731        T = regions[0]
732        try:
733            a = len(T)
734        except:
735            msg = 'Polygon_function takes a list of pairs (polygon, value). Got %s' %polygons
736            raise msg
737
738        assert a == 2, 'Must have two component each: %s' %T
739
740        self.regions = regions
741        self.default = default
742
743
744    def __call__(self, x, y):
745        from util import inside_polygon
746        from Numeric import ones, Float, concatenate, array, reshape, choose
747
748        x = array(x).astype(Float)
749        y = array(y).astype(Float)
750
751        N = len(x)
752        assert len(y) == N
753
754        points = concatenate( (reshape(x, (N, 1)),
755                               reshape(y, (N, 1))), axis=1 )
756
757        if callable(self.default):
758            z = self.default(x,y)
759        else:
760            z = ones(N, Float) * self.default
761
762        for polygon, value in self.regions:
763            indices = inside_polygon(points, polygon)
764
765            #FIXME: This needs to be vectorised
766            if callable(value):
767                for i in indices:
768                    xx = array([x[i]])
769                    yy = array([y[i]])
770                    z[i] = value(xx, yy)[0]
771            else:
772                for i in indices:
773                    z[i] = value
774
775        return z
776
778    """Read points assumed to form a polygon
779       There must be exactly two numbers in each line
780    """
781
782    #Get polygon
783    fid = open(filename)
784    lines = fid.readlines()
785    fid.close()
786    polygon = []
787    for line in lines:
788        fields = line.split(',')
789        polygon.append( [float(fields[0]), float(fields[1])] )
790
791    return polygon
792
793def populate_polygon(polygon, number_of_points, seed = None, exclude = None):
794    """Populate given polygon with uniformly distributed points.
795
796    Input:
797       polygon - list of vertices of polygon
798       number_of_points - (optional) number of points
799       seed - seed for random number generator (default=None)
800       exclude - list of polygons (inside main polygon) from where points should be excluded
801
802    Output:
803       points - list of points inside polygon
804
805    Examples:
806       populate_polygon( [[0,0], [1,0], [1,1], [0,1]], 5 )
807       will return five randomly selected points inside the unit square
808    """
809
810    from random import uniform, seed
811
812    seed(seed)
813
814    points = []
815
816    #Find outer extent of polygon
817    max_x = min_x = polygon[0][0]
818    max_y = min_y = polygon[0][1]
819    for point in polygon[1:]:
820        x = point[0]
821        if x > max_x: max_x = x
822        if x < min_x: min_x = x
823        y = point[1]
824        if y > max_y: max_y = y
825        if y < min_y: min_y = y
826
827
828    while len(points) < number_of_points:
829        x = uniform(min_x, max_x)
830        y = uniform(min_y, max_y)
831
832        append = False
833        if inside_polygon( [x,y], polygon ):
834
835            append = True
836
837            #Check exclusions
838            if exclude is not None:
839                for ex_poly in exclude:
840                    if inside_polygon( [x,y], ex_poly ):
841                        append = False
842
843
844        if append is True:
845            points.append([x,y])
846
847    return points
848
849
850
851def tsh2sww(filename, verbose=True):
852    """
853    to check if a tsh/msh file 'looks' good.
854    """
855
856    #FIXME Move to data_manager
857    from shallow_water import Domain
858    from pmesh2domain import pmesh_to_domain_instance
859    import time, os
860    from data_manager import get_dataobject
861    from os import sep, path
862    #from util import mean
863
864    if verbose == True:print 'Creating domain from', filename
865    domain = pmesh_to_domain_instance(filename, Domain)
866    if verbose == True:print "Number of triangles = ", len(domain)
867
868    domain.smooth = True
869    domain.format = 'sww'   #Native netcdf visualisation format
870    file_path, filename = path.split(filename)
871    filename, ext = path.splitext(filename)
872    domain.filename = filename
873    domain.reduction = mean
874    if verbose == True:print "file_path",file_path
875    if file_path == "":file_path = "."
877
878    if verbose == True:
879        print "Output written to " + domain.get_datadir() + sep + \
880              domain.filename + "." + domain.format
881    sww = get_dataobject(domain)
882    sww.store_connectivity()
883    sww.store_timestep('stage')
884
885####################################################################
886#Python versions of function that are also implemented in util_gateway.c
887#
888
889def gradient_python(x0, y0, x1, y1, x2, y2, q0, q1, q2):
890    """
891    """
892
893    det = (y2-y0)*(x1-x0) - (y1-y0)*(x2-x0)
894    a = (y2-y0)*(q1-q0) - (y1-y0)*(q2-q0)
895    a /= det
896
897    b = (x1-x0)*(q2-q0) - (x2-x0)*(q1-q0)
898    b /= det
899
900    return a, b
901
902
903def gradient2_python(x0, y0, x1, y1, q0, q1):
904    """Compute radient based on two points and enforce zero gradient
905    in the direction orthogonal to (x1-x0), (y1-y0)
906    """
907
908    #Old code
909    #det = x0*y1 - x1*y0
910    #if det != 0.0:
911    #    a = (y1*q0 - y0*q1)/det
912    #    b = (x0*q1 - x1*q0)/det
913
914    #Correct code (ON)
915    det = (x1-x0)**2 + (y1-y0)**2
916    if det != 0.0:
917        a = (x1-x0)*(q1-q0)/det
918        b = (y1-y0)*(q1-q0)/det
919
920    return a, b
921
922
923##############################################
924#Initialise module
925
926import compile
927if compile.can_use_C_extension('util_ext.c'):
928    from util_ext import gradient, gradient2, point_on_line
929    separate_points_by_polygon = separate_points_by_polygon_c
930else: