Changeset 6189

Timestamp:

Jan 17, 2009, 7:36:16 PM (16 years ago)

Author:

ole

Message:

Implemented interpolate_polyline in C. Hopefully, the STS file boundary will now run much faster.

Location:

anuga_core/source/anuga

Files:

• abstract_2d_finite_volumes/util.py (modified) (2 diffs)
• fit_interpolate/interpolate.py (modified) (3 diffs)
• fit_interpolate/test_interpolate.py (modified) (1 diff)
• utilities/polygon.py (modified) (3 diffs)
• utilities/polygon_ext.c (modified) (5 diffs)
• utilities/test_polygon.py (modified) (1 diff)

anuga_core/source/anuga/abstract_2d_finite_volumes/util.py

@@ -393 +393 @@
         if boundary_polygon is not None:
             #removes sts points that do not lie on boundary
+            # FIXME(Ole): Why don't we just remove such points from the list of points and associated data?
+            # I am actually convinced we can get rid of neighbour_gauge_id altogether as the sts file is produced using the ordering file.
+            # All sts points are therefore always present in the boundary. In fact, they *define* parts of the boundary.
             boundary_polygon=ensure_numeric(boundary_polygon)
             boundary_polygon[:,0] -= xllcorner
@@ -434 +437 @@
             gauge_neighbour_id=None
 
+        #print gauge_neighbour_id
+        
         if interpolation_points is not None:
             # Adjust for georef

anuga_core/source/anuga/fit_interpolate/interpolate.py

@@ -38 +38 @@
 from anuga.abstract_2d_finite_volumes.util import file_function
 from anuga.config import netcdf_mode_r, netcdf_mode_w, netcdf_mode_a
-from utilities.polygon import point_on_line
+from utilities.polygon import interpolate_polyline
 
 
@@ -507 +507 @@
         
         
-
-##
-# @brief Interpolate linearly from polyline nodes to midpoints of triangles.
-# @param f The data on the polyline nodes.
-# @param vertex_coordinates ??
-# @param gauge_neighbour_id ??
-# @param point_coordinates ??
-# @param verbose True if this function is to be verbose.
-def interpolate_polyline(f,
-                         polyline_nodes,
-                         gauge_neighbour_id,
-                         point_coordinates=None,
-                         verbose=False):
-    """Interpolate linearly between values f on polyline nodes
-    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)
-
-    z = num.zeros(len(point_coordinates), num.Float)
-
-    f = ensure_numeric(f, num.Float)
-    polyline_nodes = ensure_numeric(polyline_nodes, num.Float)
-    point_coordinates = ensure_numeric(point_coordinates, num.Float)
-    gauge_neighbour_id = ensure_numeric(gauge_neighbour_id, num.Int)
-
-    n = polyline_nodes.shape[0] # Number of nodes in polyline        
-    # Input sanity check
-    msg = 'point coordinates are not given (interpolate.py)'
-    assert point_coordinates is not None, msg
-    msg = 'function value must be specified at every interpolation node'
-    assert f.shape[0]==polyline_nodes.shape[0], msg
-    msg = 'Must define function value at one or more nodes'
-    assert f.shape[0]>0, msg
-
-
-    if n == 1:
-        msg = 'Polyline contained only one point. I need more. ' + str(f)
-        raise Exception, msg
-    elif n > 1:
-        _interpolate_polyline_aux(point_coordinates, len(point_coordinates), 
-                                  polyline_nodes, len(polyline_nodes),
-                                  f,
-                                  gauge_neighbour_id,
-                                  z)
-        
-    return z
-
-        
-def _interpolate_polyline_aux(point_coordinates, number_of_points, 
-                              polyline_nodes, number_of_nodes, 
-                              f, gauge_neighbour_id, z):
-    """Auxiliary function used by interpolate_polyline
-    """
-
-    for j in range(number_of_nodes):                
-        neighbour_id = gauge_neighbour_id[j]
-        
-        
-        if neighbour_id >= 0:
-            x0, y0 = polyline_nodes[j,:]
-            x1, y1 = polyline_nodes[neighbour_id,:]
-            
-            segment_len = sqrt((x1-x0)**2 + (y1-y0)**2)
-            segment_delta = f[neighbour_id] - f[j]            
-            slope = segment_delta/segment_len
-            
-                
-            for i in range(number_of_points):                
-                
-                x2, y2 = point_coordinates[i,:]
-                if point_on_line([x2, y2], 
-                                 [[x0, y0], [x1, y1]], 
-                                 rtol=1.0e-6):
-                                 
-
-                    dist = sqrt((x2-x0)**2 + (y2-y0)**2)
-                    z[i] = slope*dist + f[j]
-      
-
-        
         
 ##
@@ -1052 +953 @@
                                                       vertex_coordinates,
                                                       gauge_neighbour_id,
-                                                      point_coordinates=\
+                                                      interpolation_points=\
                                                           self.interpolation_points)
                         

anuga_core/source/anuga/fit_interpolate/test_interpolate.py

@@ -1846 +1846 @@
         assert num.allclose(z, answer)
         
-        
-    def test_interpolate_polyline(self):
-        """test_interpolate_polyline(self):
-        
-        This test is added under the assumption that the function interpolate_polyline implemented by
-        John Jakeman works. It has been exercised somewhat by tests of sts boundary, but never before separately.
-        """
-        
-        f = num.array([58.06150614, 58.06150614, 58.06150614])
-        vertex_coordinates = num.array([[0., 0., ],
-                                        [4.04092634, 1106.11074699],
-                                        [8.08836552, 2212.16910609]])
-        gauge_neighbour_id = [1, 2, -1]
-        point_coordinates = num.array([[2.21870766e+03, 1.09802864e+03],
-                                       [1.62739645e+03, 2.20626983e+03],
-                                       [5.20084967e+02, 2.21030386e+03],
-                                       [6.06464546e+00, 1.65913993e+03],
-                                       [1.61934862e+03, -5.88143836e+00],
-                                       [5.11996623e+02, -1.85956061e+00],
-                                       [2.02046270e+00, 5.53055373e+02]])
-                             
-        z_ref = [0., 0., 0., 58.06150614, 0., 0., 58.06150614]
-        
-        z = interpolate_polyline(f, vertex_coordinates, gauge_neighbour_id, point_coordinates)
-        assert num.allclose(z, z_ref)
-        
-        # Another f
-        f = num.array([58.06150614, 158.06150614, 258.06150614])
-        z_ref = [0., 0., 0., 208.06150645, 0., 0., 108.0615061]        
-        z = interpolate_polyline(f, vertex_coordinates, gauge_neighbour_id, point_coordinates)
-        assert num.allclose(z, z_ref)        
-                       
-
-        # Other and simpler numbers
-        f = num.array([1, 5, 13])        
-        vertex_coordinates = num.array([[0., 0., ],
-                                        [4., 4.],
-                                        [8., 8.]])        
-        point_coordinates = num.array([[0.1, 0.1],
-                                       [3.5, 3.5],
-                                       [4.0, 4.0],
-                                       [5.2, 5.2],
-                                       [7.0, 7.0],
-                                       [8.3, 8.3]])
-        gauge_neighbour_id = [1, 2, -1]
-                                               
-        z = interpolate_polyline(f, vertex_coordinates, gauge_neighbour_id, point_coordinates)
-        z_ref = [1.1, 4.5, 5., 7.4, 11., 0.]
-        #print z
-        assert num.allclose(z, z_ref)                
-        
-        # Test exception thrown for one point
-        f = num.array([5])                
-        vertex_coordinates = num.array([[4., 4.]])              
-        try:  
-            z = interpolate_polyline(f, vertex_coordinates, gauge_neighbour_id, point_coordinates)
-        except Exception:
-            pass
-        else:
-            raise Exception, 'One point should have raised exception'
                        
 #-------------------------------------------------------------

anuga_core/source/anuga/utilities/polygon.py

@@ -7 +7 @@
 from math import sqrt
 from anuga.utilities.numerical_tools import ensure_numeric
-from anuga.geospatial_data.geospatial_data import ensure_absolute
+from anuga.geospatial_data.geospatial_data import ensure_absolute, Geospatial_data
 
 
@@ -957 +957 @@
 
     return reduced_polygon
+    
+    
+        
+
+##
+# @brief Interpolate linearly from polyline nodes to midpoints of triangles.
+# @param data The data on the polyline nodes.
+# @param polyline_nodes ??
+# @param gauge_neighbour_id ??  FIXME(Ole): I want to get rid of this
+# @param point_coordinates ??
+# @param verbose True if this function is to be verbose.
+def interpolate_polyline(data,
+                         polyline_nodes,
+                         gauge_neighbour_id,
+                         interpolation_points=None,
+                         rtol=1.0e-6,
+                         atol=1.0e-8,
+                         verbose=False):
+    """Interpolate linearly between values data on polyline nodes
+    of a polyline to list of interpolation points. 
+
+    data is the data on the polyline nodes.
+
+
+    Inputs:
+      data: Vector or array of data at the polyline nodes.
+      polyline_nodes: Location of nodes where data is available.      
+      gauge_neighbour_id: ?
+      interpolation_points: 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
+      rtol, atol: Used to determine whether a point is on the polyline or not. See point_on_line.
+
+    Output:
+      Interpolated values at interpolation points
+    """
+    
+    if isinstance(interpolation_points, Geospatial_data):
+        interpolation_points = interpolation_points.get_data_points(absolute=True)
+
+    interpolated_values = num.zeros(len(interpolation_points), num.Float)
+
+    data = ensure_numeric(data, num.Float)
+    polyline_nodes = ensure_numeric(polyline_nodes, num.Float)
+    interpolation_points = ensure_numeric(interpolation_points, num.Float)
+    gauge_neighbour_id = ensure_numeric(gauge_neighbour_id, num.Int)
+
+    n = polyline_nodes.shape[0] # Number of nodes in polyline        
+    # Input sanity check
+    msg = 'interpolation_points are not given (interpolate.py)'
+    assert interpolation_points is not None, msg
+    msg = 'function value must be specified at every interpolation node'
+    assert data.shape[0]==polyline_nodes.shape[0], msg
+    msg = 'Must define function value at one or more nodes'
+    assert data.shape[0]>0, msg
+
+    if n == 1:
+        msg = 'Polyline contained only one point. I need more. ' + str(data)
+        raise Exception, msg
+    elif n > 1:
+        _interpolate_polyline(data,
+                              polyline_nodes,
+                              gauge_neighbour_id,
+                              interpolation_points,                               
+                              interpolated_values,
+                              rtol,
+                              atol)
+        
+    return interpolated_values
+
+        
+def _interpolate_polyline(data,
+                          polyline_nodes, 
+                          gauge_neighbour_id, 
+                          interpolation_points, 
+                          interpolated_values,
+                          rtol=1.0e-6,
+                          atol=1.0e-8):
+    """Auxiliary function used by interpolate_polyline
+    
+    NOTE: OBSOLETED BY C-EXTENSION
+    """
+    
+    number_of_nodes = len(polyline_nodes)                
+    number_of_points = len(interpolation_points)
+    
+    for j in range(number_of_nodes):                
+        neighbour_id = gauge_neighbour_id[j]
+        
+        # FIXME(Ole): I am convinced that gauge_neighbour_id can be discarded, but need to check with John J.
+        # Keep it for now (17 Jan 2009)
+        # When gone, we can simply interpolate between neighbouring nodes, i.e. neighbour_id = j+1.
+        # and the test below becomes something like: if j < number_of_nodes...  
+        
+        if neighbour_id >= 0:
+            x0, y0 = polyline_nodes[j,:]
+            x1, y1 = polyline_nodes[neighbour_id,:]
+            
+            segment_len = sqrt((x1-x0)**2 + (y1-y0)**2)
+            segment_delta = data[neighbour_id] - data[j]            
+            slope = segment_delta/segment_len
+            
+                
+            for i in range(number_of_points):                
+                
+                x, y = interpolation_points[i,:]
+                if point_on_line([x, y], 
+                                 [[x0, y0], [x1, y1]], 
+                                 rtol=rtol,
+                                 atol=atol):
+                                 
+
+                    dist = sqrt((x-x0)**2 + (y-y0)**2)
+                    interpolated_values[i] = slope*dist + data[j]
+      
+
+    
 
 ##############################################
@@ -966 +1083 @@
     from polygon_ext import _point_on_line
     from polygon_ext import _separate_points_by_polygon
+    from polygon_ext import _interpolate_polyline    
     #from polygon_ext import _intersection
 

anuga_core/source/anuga/utilities/polygon_ext.c

@@ -18 +18 @@
 #include "math.h"
 
+double dist(double x,
+            double y) {
+  
+  return sqrt(x*x + y*y);
+}
+
 
 int __point_on_line(double x, double y,
@@ -66 +72 @@
     // subject to specified absolute tolerance
 
-    len_a = sqrt(a0*a0 + a1*a1);
-    len_b = sqrt(b0*b0 + b1*b1);
+    len_a = dist(a0, a1); //sqrt(a0*a0 + a1*a1);
+    len_b = dist(b0, b1); //sqrt(b0*b0 + b1*b1);
 
     if (a0*b0 + a1*b1 >= 0 && len_a <= len_b) {
@@ -184 +190 @@
 } 
 */
+
+
+
+int __interpolate_polyline(int number_of_nodes,
+                           int number_of_points,
+                           double* data,
+                           double* polyline_nodes,
+                           long* gauge_neighbour_id,
+                           double* interpolation_points,                               
+                           double* interpolated_values,
+                           double rtol,
+                           double atol) {
+                           
+  int j, i, neighbour_id;
+  double x0, y0, x1, y1, x, y;
+  double segment_len, segment_delta, slope, alpha;
+
+  for (j=0; j<number_of_nodes; j++) {  
+
+    neighbour_id = gauge_neighbour_id[j];
+        
+    // FIXME(Ole): I am convinced that gauge_neighbour_id can be discarded, but need to check with John J.
+    // Keep it for now (17 Jan 2009)
+    // When gone, we can simply interpolate between neighbouring nodes, i.e. neighbour_id = j+1.
+    // and the test below becomes something like: if j < number_of_nodes...  
+        
+    if (neighbour_id >= 0) {
+      x0 = polyline_nodes[2*j];
+      y0 = polyline_nodes[2*j+1];
+      
+      x1 = polyline_nodes[2*neighbour_id];
+      y1 = polyline_nodes[2*neighbour_id+1];      
+      
+            
+      segment_len = dist(x1-x0, y1-y0);
+      segment_delta = data[neighbour_id] - data[j];            
+      slope = segment_delta/segment_len;
+            
+      for (i=0; i<number_of_points; i++) {                
+        x = interpolation_points[2*i];
+        y = interpolation_points[2*i+1];        
+        
+        if (__point_on_line(x, y, x0, y0, x1, y1, rtol, atol)) {
+          alpha = dist(x-x0, y-y0);
+          interpolated_values[i] = slope*alpha + data[j];
+        }
+      }
+    }
+  }
+                           
+  return 0;                          
+}                                                      
 
 
@@ -304 +362 @@
 
 
+
+// Gateways to Python
+PyObject *_interpolate_polyline(PyObject *self, PyObject *args) {
+  //
+  // _interpolate_polyline(data, polyline_nodes, gauge_neighbour_id, interpolation_points
+  //                       interpolated_values):
+  //
+
+  
+  PyArrayObject
+    *data,
+    *polyline_nodes,
+    *gauge_neighbour_id,
+    *interpolation_points,
+    *interpolated_values;
+
+  double rtol, atol;  
+  int number_of_nodes, number_of_points, res;
+  
+  // Convert Python arguments to C
+  if (!PyArg_ParseTuple(args, "OOOOOdd",
+                        &data,
+                        &polyline_nodes,
+                        &gauge_neighbour_id,
+                        &interpolation_points,
+                        &interpolated_values,
+                        &rtol,
+                        &atol)) {
+    
+    PyErr_SetString(PyExc_RuntimeError, 
+                    "_interpolate_polyline could not parse input");
+    return NULL;
+  }
+
+  number_of_nodes = polyline_nodes -> dimensions[0];  // Number of nodes in polyline
+  number_of_points = interpolation_points -> dimensions[0];   //Number of points
+  
+
+  // Call underlying routine
+  res = __interpolate_polyline(number_of_nodes,
+                               number_of_points,
+                               (double*) data -> data,
+                               (double*) polyline_nodes -> data,
+                               (long*) gauge_neighbour_id -> data,
+                               (double*) interpolation_points -> data,                         
+                               (double*) interpolated_values -> data,
+                               rtol,
+                               atol);                                                  
+
+  // Return
+  return Py_BuildValue("");  
+}
+
+
+
 /*
 PyObject *_intersection(PyObject *self, PyObject *args) {
@@ -424 +537 @@
   {"_separate_points_by_polygon", _separate_points_by_polygon, 
                                  METH_VARARGS, "Print out"},
+  {"_interpolate_polyline", _interpolate_polyline, 
+                                 METH_VARARGS, "Print out"},                             
   {NULL, NULL, 0, NULL}   /* sentinel */
 };

anuga_core/source/anuga/utilities/test_polygon.py

@@ -1690 +1690 @@
 
         
+    def test_interpolate_polyline(self):
+        """test_interpolate_polyline(self):
+        
+        This test is added under the assumption that the function interpolate_polyline implemented by
+        John Jakeman works. It has been exercised somewhat by tests of sts boundary, but never before separately.
+        """
+        
+        f = num.array([58.06150614, 58.06150614, 58.06150614])
+        vertex_coordinates = num.array([[0., 0., ],
+                                        [4.04092634, 1106.11074699],
+                                        [8.08836552, 2212.16910609]])
+        gauge_neighbour_id = [1, 2, -1]
+        point_coordinates = num.array([[2.21870766e+03, 1.09802864e+03],
+                                       [1.62739645e+03, 2.20626983e+03],
+                                       [5.20084967e+02, 2.21030386e+03],
+                                       [6.06464546e+00, 1.65913993e+03],
+                                       [1.61934862e+03, -5.88143836e+00],
+                                       [5.11996623e+02, -1.85956061e+00],
+                                       [2.02046270e+00, 5.53055373e+02]])
+                             
+        z_ref = [0., 0., 0., 58.06150614, 0., 0., 58.06150614]
+        
+        z = interpolate_polyline(f, vertex_coordinates, gauge_neighbour_id, point_coordinates)
+        assert num.allclose(z, z_ref)
+        
+        # Another f
+        f = num.array([58.06150614, 158.06150614, 258.06150614])
+        z_ref = [0., 0., 0., 208.06150645, 0., 0., 108.0615061]        
+        z = interpolate_polyline(f, vertex_coordinates, gauge_neighbour_id, point_coordinates)
+        assert num.allclose(z, z_ref)        
+                       
+
+        # Other and simpler numbers
+        f = num.array([1, 5, 13])        
+        vertex_coordinates = num.array([[0., 0., ],
+                                        [4., 4.],
+                                        [8., 8.]])        
+        point_coordinates = num.array([[0.1, 0.1],
+                                       [3.5, 3.5],
+                                       [4.0, 4.0],
+                                       [5.2, 5.2],
+                                       [7.0, 7.0],
+                                       [8.3, 8.3]])
+        gauge_neighbour_id = [1, 2, -1]
+                                               
+        z = interpolate_polyline(f, vertex_coordinates, gauge_neighbour_id, point_coordinates)
+        z_ref = [1.1, 4.5, 5., 7.4, 11., 0.]
+        #print z
+        assert num.allclose(z, z_ref)                
+        
+        # Test exception thrown for one point
+        f = num.array([5])                
+        vertex_coordinates = num.array([[4., 4.]])              
+        try:  
+            z = interpolate_polyline(f, vertex_coordinates, gauge_neighbour_id, point_coordinates)
+        except Exception:
+            pass
+        else:
+            raise Exception, 'One point should have raised exception'
+        
+
+        # More polyline nodes
+        data = num.array([1, 5, 13, 12, 6, 29])        
+        polyline_nodes = num.array([[0., 0.],
+                                    [4., 4.],
+                                    [8., 8.],
+                                    [10., 10.],
+                                    [10., 5.],
+                                    [10., 0.]])
+                                                
+        point_coordinates = num.array([[0.1, 0.1],
+                                       [3.5, 3.5],
+                                       [4.0, 4.0],
+                                       [5.2, 5.2],
+                                       [7.0, 7.0],
+                                       [8.3, 8.3],
+                                       [10., 10.],
+                                       [10., 9.],
+                                       [10., 7.1],
+                                       [10., 4.3],
+                                       [10., 1.0]])
+                                       
+        gauge_neighbour_id = [1, 2, 3, 4, 5, -1]
+                                               
+        z = interpolate_polyline(data, polyline_nodes, gauge_neighbour_id, point_coordinates)
+        z_ref = [1.1, 4.5, 5., 7.4, 11., 12.85, 12., 10.8, 8.52, 9.22, 24.4]
+        assert num.allclose(z, z_ref)                
+        
+        
+                
 #-------------------------------------------------------------
 if __name__ == "__main__":