Changeset 5403

Timestamp:

Jun 13, 2008, 3:48:16 PM (15 years ago)

Author:

ole

Message:

Added tolerances for point_on_line and changed interpolate_polyline to use a general relative tolerance.

Location:

anuga_core/source/anuga

Files:

• fit_interpolate/interpolate.py (modified) (3 diffs)
• shallow_water/data_manager.py (modified) (1 diff)
• utilities/polygon.py (modified) (2 diffs)
• utilities/polygon_ext.c (modified) (7 diffs)
• utilities/test_polygon.py (modified) (1 diff)

anuga_core/source/anuga/fit_interpolate/interpolate.py

@@ -96 +96 @@
                              point_coordinates=None,
                              verbose=False):
- 
+        """Interpolate linearly between values f on nodes (vertex coordinates) of a polyline to midpoints of triangles
+        of boundary.
+
+        f is the data on the polyline nodes.
+
+        The mesh values representing a smooth surface are
+        assumed to be specified in f.
+
+        Inputs:
+          f: Vector or array of data at the polyline nodes.
+              If f is an array, interpolation will be done for each column as
+              per underlying matrix-matrix multiplication
+          point_coordinates: Interpolate polyline data to these positions.
+              List of coordinate pairs [x, y] of
+              data points or an nx2 Numeric array or a Geospatial_data object
+              
+        Output:
+          Interpolated values at inputted points (z).
+        """
+        
+
+        # FIXME: There is an option of passing a tolerance into this
+        
         if isinstance(point_coordinates, Geospatial_data):
             point_coordinates = point_coordinates.get_data_points( \
                 absolute = True)
  
-        from utilities.polygon import point_on_line,point_on_line_py
+        from utilities.polygon import point_on_line
         from Numeric import ones
         z=ones(len(point_coordinates),Float)
@@ -108 +130 @@
         assert point_coordinates is not None, msg
         msg='function value must be specified at every interpolation node'
-        assert f.shape[0]==vertex_coordinates.shape[0],msg
+        assert f.shape[0]==vertex_coordinates.shape[0], msg
         msg='Must define function value at one or more nodes'
-        assert f.shape[0]>0,msg
+        assert f.shape[0]>0, msg
 
         n=f.shape[0]
         if n==1:
             z=f*z
+            msg = 'Polyline contained only one point. I need more. ', str(f)
+            raise Exception, msg
+            
+        # FIXME (John): add unit test for only One vertex point. Exception should be thrown.
+
+        
         elif n>1:
             for i in range(len(point_coordinates)):
@@ -120 +148 @@
                 for j in range(n):
                     if gauge_neighbour_id[j]>=0:
-                        if point_on_line_py(point_coordinates[i],[vertex_coordinates[j],vertex_coordinates[gauge_neighbour_id[j]]]):
+                        if point_on_line(point_coordinates[i],
+                                         [vertex_coordinates[j], vertex_coordinates[gauge_neighbour_id[j]]],
+                                         rtol=1.0e-6):
                             found=True
                             x0=vertex_coordinates[j][0]

anuga_core/source/anuga/shallow_water/data_manager.py

@@ -5413 +5413 @@
     #Title
     fid.write('location_index'+d+'lat'+d+ 'long' +d+ 'Easting' +d+ \
-                  'Northing' +d+ 'depth m' + "\n")
+              'Northing' +d+ 'depth m' + "\n")
     i = 0
     for depth, point_utm, lat, long in map(None, depths,

anuga_core/source/anuga/utilities/polygon.py

@@ -17 +17 @@
 from anuga.geospatial_data.geospatial_data import ensure_absolute
 
-def point_on_line_py(point, line):
-    from Numeric import fabs
-    point = ensure_numeric(point)
-    line = ensure_numeric(line)
-
-
-    x=point[0];y=point[1]
-    x0=line[0,0];y0=line[0,1]
-    x1=line[1,0];y1=line[1,1]
-    #from pylab import plot,show
-    #plot(line[:,0],line[:,1])
-    #plot([x],[y],'o')
-    #show()
-    a0 = x - x0;
-    a1 = y - y0;
-
-    a_normal0 = a1;
-    a_normal1 = -a0;
-
-    b0 = x1 - x0;
-    b1 = y1 - y0;
-
-    #if ( a_normal0*b0 + a_normal1*b1 == 0.0 ):
-    eps=200
-    #print 'normal',a_normal0*b0 + a_normal1*b1
-    if ( fabs(a_normal0*b0 + a_normal1*b1) < eps ):
-        #for some reason (perhaps catastrophic cancellation) urs_boundary.py
-        #example only works for eps=2
-        #point is somewhere on the infinite extension of the line
-        # FIXME (Ole): Perhaps add a tolerance here instead of 0.0 
-
-        len_a = sqrt(a0*a0 + a1*a1);
-        len_b = sqrt(b0*b0 + b1*b1);
-
-        if (a0*b0 + a1*b1 >= 0 and len_a <= len_b):
-            return 1
-        else:
-            return 0
-    else:
-        return 0
-
-
-def point_on_line(point, line):
+
+def point_on_line(point, line, rtol=0.0, atol=0.0):
     """Determine whether a point is on a line segment
 
@@ -68 +27 @@
 
     Output:
-        
-    """
+
+    Note: Line can be degenerate and function still works to discern coinciding points from non-coinciding.
+    """
+
+    # FIXME(Ole): Perhaps make defaults as in allclose: rtol=1.0e-5, atol=1.0e-8
 
     point = ensure_numeric(point)
     line = ensure_numeric(line)
 
-
     res = _point_on_line(point[0], point[1],
                          line[0,0], line[0,1],
-                         line[1,0], line[1,1])
+                         line[1,0], line[1,1],
+                         rtol, atol)
     
     return bool(res)

anuga_core/source/anuga/utilities/polygon_ext.c

@@ -21 +21 @@
 int __point_on_line(double x, double y,
                     double x0, double y0,
-                    double x1, double y1) {
+                    double x1, double y1,
+                    double rtol,
+                    double atol) {
   /*Determine whether a point is on a line segment
 
@@ -31 +33 @@
 
   double a0, a1, a_normal0, a_normal1, b0, b1, len_a, len_b;
+  double nominator, denominator;
+  int is_parallel;
 
   a0 = x - x0;
@@ -41 +45 @@
   b1 = y1 - y0;
 
-  if ( a_normal0*b0 + a_normal1*b1 == 0.0 ) {
-    //Point is somewhere on the infinite extension of the line
-    // FIXME (Ole): Perhaps add a tolerance here instead of 0.0 
+  nominator = fabs(a_normal0*b0 + a_normal1*b1);
+  denominator = b0*b0 + b1*b1;
+  
+  // Determine if line is parallel to point vector up to a tolerance
+  is_parallel = 0;
+  if (denominator == 0.0) {
+    // Use absolute tolerance
+    if (nominator <= atol) {
+      is_parallel = 1;
+    }
+  } else {
+    // Denominator is positive - use relative tolerance
+    if (nominator/denominator <= rtol) {
+      is_parallel = 1;
+    }    
+  }
+    
+  if (is_parallel) {
+    // Point is somewhere on the infinite extension of the line
+    // subject to specified absolute tolerance
 
     len_a = sqrt(a0*a0 + a1*a1);
@@ -173 +194 @@
                                 int verbose) {
 
-  double minpx, maxpx, minpy, maxpy, x, y, px_i, py_i, px_j, py_j;
+  double minpx, maxpx, minpy, maxpy, x, y, px_i, py_i, px_j, py_j, rtol=0.0, atol=0.0;
   int i, j, k, outside_index, inside_index, inside;
 
@@ -223 +244 @@
 
         //Check for case where point is contained in line segment
-        if (__point_on_line(x, y, px_i, py_i, px_j, py_j)) {
+        if (__point_on_line(x, y, px_i, py_i, px_j, py_j, rtol, atol)) {
           if (closed == 1) {
             inside = 1;
@@ -260 +281 @@
   //
 
-  double x, y, x0, y0, x1, y1;
+  double x, y, x0, y0, x1, y1, rtol, atol;
   int res;
   PyObject *result;
 
   // Convert Python arguments to C
-  if (!PyArg_ParseTuple(args, "dddddd", &x, &y, &x0, &y0, &x1, &y1)) {
+  if (!PyArg_ParseTuple(args, "dddddddd", &x, &y, &x0, &y0, &x1, &y1, &rtol, &atol)) {
     PyErr_SetString(PyExc_RuntimeError, 
                     "point_on_line could not parse input");    
@@ -273 +294 @@
 
   // Call underlying routine
-  res = __point_on_line(x, y, x0, y0, x1, y1);
+  res = __point_on_line(x, y, x0, y0, x1, y1, rtol, atol);
 
   // Return values a and b

anuga_core/source/anuga/utilities/test_polygon.py

@@ -165 +165 @@
         assert not point_on_line( [40,19], [[40,20], [40,40]] )
 
+        # Degenerate line
+        assert point_on_line( [40,19], [[40,19], [40,19]] )
+        assert not point_on_line( [40,19], [[40,40], [40,40]] )