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

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

Made interpolation_points in file_function absolute UTM coordinates and wrote test

File size: 26.4 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/l
19    v2 = v/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:
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, -a])
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 absolute 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 as absolute UTM coordinates
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
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    from util import ensure_numeric
232
233    #Open NetCDF file
234    if verbose: print 'Reading', filename
235    fid = NetCDFFile(filename, 'r')
236
237    if type(quantity_names) == types.StringType:
238        quantity_names = [quantity_names]
239
240    if quantity_names is None or len(quantity_names) < 1:
241        #If no quantities are specified get quantities from file
242        #x, y, time are assumed as the independent variables so
243        #they are excluded as potentiol quantities
244        quantity_names = []
245        for name in fid.variables.keys():
246            if name not in ['x', 'y', 'time']:
247                quantity_names.append(name)
248
249    if len(quantity_names) < 1:
250        msg = 'ERROR: At least one independent value must be specified'
251        raise msg
252
253
254    if interpolation_points is not None:
255        interpolation_points = ensure_numeric(interpolation_points)
256
257
258
259    #Now assert that requested quantitites (and the independent ones)
260    #are present in file
261    missing = []
262    for quantity in ['time'] + quantity_names:
263        if not fid.variables.has_key(quantity):
264            missing.append(quantity)
265
266    if len(missing) > 0:
267        msg = 'Quantities %s could not be found in file %s'\
268              %(str(missing), filename)
269        fid.close()
270        raise msg
271
272    #Decide whether this data has a spatial dimension
273    spatial = True
274    for quantity in ['x', 'y']:
275        if not fid.variables.has_key(quantity):
276            spatial = False
277
278    if filename[-3:] == 'tms' and spatial is True:
279        msg = 'Files of type tms must not contain spatial information'
280        raise msg
281
282    if filename[-3:] == 'sww' and spatial is False:
283        msg = 'Files of type sww must contain spatial information'
284        raise msg
285
286    #Get first timestep
287    try:
288        starttime = fid.starttime
289    except ValueError:
290        msg = 'Could not read starttime from file %s' %filename
291        raise msg
292
293    #Get variables
294    if verbose: print 'Get variables'
295    time = fid.variables['time'][:]
296
297    if spatial:
298        #Get origin
299        xllcorner = fid.xllcorner
300        yllcorner = fid.yllcorner
301        zone = fid.zone
302
303        x = fid.variables['x'][:]
304        y = fid.variables['y'][:]
305        triangles = fid.variables['volumes'][:]
306
307        x = reshape(x, (len(x),1))
308        y = reshape(y, (len(y),1))
309        vertex_coordinates = concatenate((x,y), axis=1) #m x 2 array
310
311        if interpolation_points is not None:
313            interpolation_points[:,0] -= xllcorner
314            interpolation_points[:,1] -= yllcorner
315
316
317
318
319    if domain is not None:
320        #Update domain.startime if it is *earlier* than starttime
321        if starttime > domain.starttime:
322            msg = 'WARNING: Start time as specified in domain (%f)'\
323                  %domain.starttime
324            msg += ' is earlier than the starttime of file %s (%f).'\
325                     %(filename, starttime)
326            msg += ' Modifying domain starttime accordingly.'
327
328            if verbose: print msg
329            domain.starttime = starttime #Modifying model time
330            if verbose: print 'Domain starttime is now set to %f'\
331               %domain.starttime
332
333
334        #If domain.startime is *later* than starttime,
335        #move time back - relative to domain's time
336        if domain.starttime > starttime:
337            time = time - domain.starttime + starttime
338
339        #FIXME Use method in geo to reconcile
340        #if spatial:
341        #assert domain.geo_reference.xllcorner == xllcorner
342        #assert domain.geo_reference.yllcorner == yllcorner
343        #assert domain.geo_reference.zone == zone
344
345    if verbose:
346        print 'File_function data obtained from: %s' %filename
347        print '  References:'
348        #print '    Datum ....' #FIXME
349        if spatial:
350            print '    Lower left corner: [%f, %f]'\
351                  %(xllcorner, yllcorner)
352        print '    Start time:   %f' %starttime
353
354
355    #Produce values for desired data points at
356    #each timestep
357
358    quantities = {}
359    for i, name in enumerate(quantity_names):
360        quantities[name] = fid.variables[name][:]
361    fid.close()
362
363
364    from least_squares import Interpolation_function
365
366    if not spatial:
367        vertex_coordinates = triangles = interpolation_points = None
368
369    return Interpolation_function(time,
370                                  quantities,
371                                  quantity_names,
372                                  vertex_coordinates,
373                                  triangles,
374                                  interpolation_points,
375                                  verbose = verbose)
376
377
378
379
380
382    """Read simple xya file, possibly with a header in the first line, just
383    x y [attributes]
384    separated by whitespace
385
386    Format can be either ASCII or NetCDF
387
388    Return list of points, list of attributes and
389    attribute names if present otherwise None
390    """
391    #FIXME: Probably obsoleted by similar function in load_ASCII
393
394    from Scientific.IO.NetCDF import NetCDFFile
395
396    if format.lower() == 'netcdf':
397        #Get NetCDF
398
399        fid = NetCDFFile(filename, 'r')
400
401        # Get the variables
402        points = fid.variables['points'][:]
403        keys = fid.variables.keys()
404        attributes = {}
405        for key in keys:
406            if key != 'points':
407                attributes[key] = fid.variables[key][:]
408
409        fid.close()
410    else:
411        #Read as ASCII file assuming that it is separated by whitespace
412        fid = open(filename)
414        fid.close()
415
416        #Check if there is a header line
417        fields = lines.strip().split()
418        try:
419            float(fields)
420        except:
421            #This must be a header line
422            attribute_names = fields
423            lines = lines[1:]
424        else:
425            attribute_names = ['elevation']  #HACK
426
427        attributes = {}
428        for key in attribute_names:
429            attributes[key] = []
430
431        points = []
432        for line in lines:
433            fields = line.strip().split()
434            points.append( (float(fields), float(fields)) )
435            for i, key in enumerate(attribute_names):
436                attributes[key].append( float(fields[2+i]) )
437
438    return points, attributes
439
440
441#####################################
442#POLYGON STUFF
443#
444#FIXME: All these should be put into new module polygon.py
445
446
447
448#These have been redefined to use separate_by_polygon which will
449#put all inside indices in the first part of the indices array and
450#outside indices in the last
451
452def inside_polygon(points, polygon, closed = True, verbose = False):
453    """See separate_points_by_polygon for documentation
454    """
455
456    from Numeric import array, Float, reshape
457
458    if verbose: print 'Checking input to inside_polygon'
459    try:
460        points = ensure_numeric(points, Float)
461    except:
462        msg = 'Points could not be converted to Numeric array'
463        raise msg
464
465    try:
466        polygon = ensure_numeric(polygon, Float)
467    except:
468        msg = 'Polygon could not be converted to Numeric array'
469        raise msg
470
471
472
473    if len(points.shape) == 1:
474        one_point = True
475        points = reshape(points, (1,2))
476    else:
477        one_point = False
478
479    indices, count = separate_points_by_polygon(points, polygon,
480                                                closed, verbose)
481
482    if one_point:
483        return count == 1
484    else:
485        return indices[:count]
486
487def outside_polygon(points, polygon, closed = True, verbose = False):
488    """See separate_points_by_polygon for documentation
489    """
490
491    from Numeric import array, Float, reshape
492
493    if verbose: print 'Checking input to outside_polygon'
494    try:
495        points = ensure_numeric(points, Float)
496    except:
497        msg = 'Points could not be converted to Numeric array'
498        raise msg
499
500    try:
501        polygon = ensure_numeric(polygon, Float)
502    except:
503        msg = 'Polygon could not be converted to Numeric array'
504        raise msg
505
506
507
508    if len(points.shape) == 1:
509        one_point = True
510        points = reshape(points, (1,2))
511    else:
512        one_point = False
513
514    indices, count = separate_points_by_polygon(points, polygon,
515                                                closed, verbose)
516
517    if one_point:
518        return count != 1
519    else:
520        return indices[count:][::-1]  #return reversed
521
522
523def separate_points_by_polygon(points, polygon,
524                               closed = True, verbose = False):
525    """Determine whether points are inside or outside a polygon
526
527    Input:
528       points - Tuple of (x, y) coordinates, or list of tuples
529       polygon - list of vertices of polygon
530       closed - (optional) determine whether points on boundary should be
531       regarded as belonging to the polygon (closed = True)
532       or not (closed = False)
533
534    Outputs:
535       indices: array of same length as points with indices of points falling
536       inside the polygon listed from the beginning and indices of points
537       falling outside listed from the end.
538
539       count: count of points falling inside the polygon
540
541       The indices of points inside are obtained as indices[:count]
542       The indices of points outside are obtained as indices[count:]
543
544
545    Examples:
546       U = [[0,0], [1,0], [1,1], [0,1]] #Unit square
547
548       separate_points_by_polygon( [[0.5, 0.5], [1, -0.5], [0.3, 0.2]], U)
549       will return the indices [0, 2, 1] and count == 2 as only the first
550       and the last point are inside the unit square
551
552    Remarks:
553       The vertices may be listed clockwise or counterclockwise and
554       the first point may optionally be repeated.
555       Polygons do not need to be convex.
556       Polygons can have holes in them and points inside a hole is
557       regarded as being outside the polygon.
558
559    Algorithm is based on work by Darel Finley,
560    http://www.alienryderflex.com/polygon/
561
562    """
563
564    from Numeric import array, Float, reshape, Int, zeros
565
566
567    #Input checks
568    try:
569        points = ensure_numeric(points, Float)
570    except:
571        msg = 'Points could not be converted to Numeric array'
572        raise msg
573
574    try:
575        polygon = ensure_numeric(polygon, Float)
576    except:
577        msg = 'Polygon could not be converted to Numeric array'
578        raise msg
579
580    assert len(polygon.shape) == 2,\
581       'Polygon array must be a 2d array of vertices'
582
583    assert polygon.shape == 2,\
584       'Polygon array must have two columns'
585
586    assert len(points.shape) == 2,\
587       'Points array must be a 2d array'
588
589    assert points.shape == 2,\
590       'Points array must have two columns'
591
592    N = polygon.shape #Number of vertices in polygon
593    M = points.shape  #Number of points
594
595    px = polygon[:,0]
596    py = polygon[:,1]
597
598    #Used for an optimisation when points are far away from polygon
599    minpx = min(px); maxpx = max(px)
600    minpy = min(py); maxpy = max(py)
601
602
603    #Begin main loop
604    indices = zeros(M, Int)
605
606    inside_index = 0    #Keep track of points inside
607    outside_index = M-1 #Keep track of points outside (starting from end)
608
609    for k in range(M):
610
611        if verbose:
612            if k %((M+10)/10)==0: print 'Doing %d of %d' %(k, M)
613
614        #for each point
615        x = points[k, 0]
616        y = points[k, 1]
617
618        inside = False
619
620        if not x > maxpx or x < minpx or y > maxpy or y < minpy:
621            #Check polygon
622            for i in range(N):
623                j = (i+1)%N
624
625                #Check for case where point is contained in line segment
626                if point_on_line(x, y, px[i], py[i], px[j], py[j]):
627                    if closed:
628                        inside = True
629                    else:
630                        inside = False
631                    break
632                else:
633                    #Check if truly inside polygon
634                    if py[i] < y and py[j] >= y or\
635                       py[j] < y and py[i] >= y:
636                        if px[i] + (y-py[i])/(py[j]-py[i])*(px[j]-px[i]) < x:
637                            inside = not inside
638
639        if inside:
640            indices[inside_index] = k
641            inside_index += 1
642        else:
643            indices[outside_index] = k
644            outside_index -= 1
645
646    return indices, inside_index
647
648
649def separate_points_by_polygon_c(points, polygon,
650                                 closed = True, verbose = False):
651    """Determine whether points are inside or outside a polygon
652
653    C-wrapper
654    """
655
656    from Numeric import array, Float, reshape, zeros, Int
657
658
659    if verbose: print 'Checking input to separate_points_by_polygon'
660    #Input checks
661    try:
662        points = ensure_numeric(points, Float)
663    except:
664        msg = 'Points could not be converted to Numeric array'
665        raise msg
666
667    #if verbose: print 'Checking input to separate_points_by_polygon 2'
668    try:
669        polygon = ensure_numeric(polygon, Float)
670    except:
671        msg = 'Polygon could not be converted to Numeric array'
672        raise msg
673
674    if verbose: print 'check'
675
676    assert len(polygon.shape) == 2,\
677       'Polygon array must be a 2d array of vertices'
678
679    assert polygon.shape == 2,\
680       'Polygon array must have two columns'
681
682    assert len(points.shape) == 2,\
683       'Points array must be a 2d array'
684
685    assert points.shape == 2,\
686       'Points array must have two columns'
687
688    N = polygon.shape #Number of vertices in polygon
689    M = points.shape  #Number of points
690
691    from util_ext import separate_points_by_polygon
692
693    if verbose: print 'Allocating array for indices'
694
695    indices = zeros( M, Int )
696
697    if verbose: print 'Calling C-version of inside poly'
698    count = separate_points_by_polygon(points, polygon, indices,
699                                       int(closed), int(verbose))
700
701    return indices, count
702
703
704
705class Polygon_function:
706    """Create callable object f: x,y -> z, where a,y,z are vectors and
707    where f will return different values depending on whether x,y belongs
708    to specified polygons.
709
710    To instantiate:
711
712       Polygon_function(polygons)
713
714    where polygons is a list of tuples of the form
715
716      [ (P0, v0), (P1, v1), ...]
717
718      with Pi being lists of vertices defining polygons and vi either
719      constants or functions of x,y to be applied to points with the polygon.
720
721    The function takes an optional argument, default which is the value
722    (or function) to used for points not belonging to any polygon.
723    For example:
724
725       Polygon_function(polygons, default = 0.03)
726
727    If omitted the default value will be 0.0
728
729    Note: If two polygons overlap, the one last in the list takes precedence
730
731    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???
732
733    """
734
735    def __init__(self, regions, default = 0.0):
736
737        try:
738            len(regions)
739        except:
740            msg = 'Polygon_function takes a list of pairs (polygon, value). Got %s' %polygons
741            raise msg
742
743
744        T = regions
745        try:
746            a = len(T)
747        except:
748            msg = 'Polygon_function takes a list of pairs (polygon, value). Got %s' %polygons
749            raise msg
750
751        assert a == 2, 'Must have two component each: %s' %T
752
753        self.regions = regions
754        self.default = default
755
756
757    def __call__(self, x, y):
758        from util import inside_polygon
759        from Numeric import ones, Float, concatenate, array, reshape, choose
760
761        x = array(x).astype(Float)
762        y = array(y).astype(Float)
763
764        N = len(x)
765        assert len(y) == N
766
767        points = concatenate( (reshape(x, (N, 1)),
768                               reshape(y, (N, 1))), axis=1 )
769
770        if callable(self.default):
771            z = self.default(x,y)
772        else:
773            z = ones(N, Float) * self.default
774
775        for polygon, value in self.regions:
776            indices = inside_polygon(points, polygon)
777
778            #FIXME: This needs to be vectorised
779            if callable(value):
780                for i in indices:
781                    xx = array([x[i]])
782                    yy = array([y[i]])
783                    z[i] = value(xx, yy)
784            else:
785                for i in indices:
786                    z[i] = value
787
788        return z
789
791    """Read points assumed to form a polygon
792       There must be exactly two numbers in each line
793    """
794
795    #Get polygon
796    fid = open(filename)
798    fid.close()
799    polygon = []
800    for line in lines:
801        fields = line.split(',')
802        polygon.append( [float(fields), float(fields)] )
803
804    return polygon
805
806def populate_polygon(polygon, number_of_points, seed = None, exclude = None):
807    """Populate given polygon with uniformly distributed points.
808
809    Input:
810       polygon - list of vertices of polygon
811       number_of_points - (optional) number of points
812       seed - seed for random number generator (default=None)
813       exclude - list of polygons (inside main polygon) from where points should be excluded
814
815    Output:
816       points - list of points inside polygon
817
818    Examples:
819       populate_polygon( [[0,0], [1,0], [1,1], [0,1]], 5 )
820       will return five randomly selected points inside the unit square
821    """
822
823    from random import uniform, seed
824
825    seed(seed)
826
827    points = []
828
829    #Find outer extent of polygon
830    max_x = min_x = polygon
831    max_y = min_y = polygon
832    for point in polygon[1:]:
833        x = point
834        if x > max_x: max_x = x
835        if x < min_x: min_x = x
836        y = point
837        if y > max_y: max_y = y
838        if y < min_y: min_y = y
839
840
841    while len(points) < number_of_points:
842        x = uniform(min_x, max_x)
843        y = uniform(min_y, max_y)
844
845        append = False
846        if inside_polygon( [x,y], polygon ):
847
848            append = True
849
850            #Check exclusions
851            if exclude is not None:
852                for ex_poly in exclude:
853                    if inside_polygon( [x,y], ex_poly ):
854                        append = False
855
856
857        if append is True:
858            points.append([x,y])
859
860    return points
861
862
863
864def tsh2sww(filename, verbose=True):
865    """
866    to check if a tsh/msh file 'looks' good.
867    """
868
869    #FIXME Move to data_manager
870    from shallow_water import Domain
871    from pmesh2domain import pmesh_to_domain_instance
872    import time, os
873    from data_manager import get_dataobject
874    from os import sep, path
875    #from util import mean
876
877    if verbose == True:print 'Creating domain from', filename
878    domain = pmesh_to_domain_instance(filename, Domain)
879    if verbose == True:print "Number of triangles = ", len(domain)
880
881    domain.smooth = True
882    domain.format = 'sww'   #Native netcdf visualisation format
883    file_path, filename = path.split(filename)
884    filename, ext = path.splitext(filename)
885    domain.filename = filename
886    domain.reduction = mean
887    if verbose == True:print "file_path",file_path
888    if file_path == "":file_path = "."
890
891    if verbose == True:
892        print "Output written to " + domain.get_datadir() + sep + \
893              domain.filename + "." + domain.format
894    sww = get_dataobject(domain)
895    sww.store_connectivity()
896    sww.store_timestep('stage')
897
898####################################################################
899#Python versions of function that are also implemented in util_gateway.c
900#
901
902def gradient_python(x0, y0, x1, y1, x2, y2, q0, q1, q2):
903    """
904    """
905
906    det = (y2-y0)*(x1-x0) - (y1-y0)*(x2-x0)
907    a = (y2-y0)*(q1-q0) - (y1-y0)*(q2-q0)
908    a /= det
909
910    b = (x1-x0)*(q2-q0) - (x2-x0)*(q1-q0)
911    b /= det
912
913    return a, b
914
915
916def gradient2_python(x0, y0, x1, y1, q0, q1):
918    in the direction orthogonal to (x1-x0), (y1-y0)
919    """
920
921    #Old code
922    #det = x0*y1 - x1*y0
923    #if det != 0.0:
924    #    a = (y1*q0 - y0*q1)/det
925    #    b = (x0*q1 - x1*q0)/det
926
927    #Correct code (ON)
928    det = (x1-x0)**2 + (y1-y0)**2
929    if det != 0.0:
930        a = (x1-x0)*(q1-q0)/det
931        b = (y1-y0)*(q1-q0)/det
932
933    return a, b
934
935
936##############################################
937#Initialise module
938
939import compile
940if compile.can_use_C_extension('util_ext.c'):