Changeset 6189


Ignore:
Timestamp:
Jan 17, 2009, 7:36:16 PM (10 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:
6 edited

Legend:

Unmodified
Added
Removed
  • anuga_core/source/anuga/abstract_2d_finite_volumes/util.py

    r6178 r6189  
    393393        if boundary_polygon is not None:
    394394            #removes sts points that do not lie on boundary
     395            # FIXME(Ole): Why don't we just remove such points from the list of points and associated data?
     396            # I am actually convinced we can get rid of neighbour_gauge_id altogether as the sts file is produced using the ordering file.
     397            # All sts points are therefore always present in the boundary. In fact, they *define* parts of the boundary.
    395398            boundary_polygon=ensure_numeric(boundary_polygon)
    396399            boundary_polygon[:,0] -= xllcorner
     
    434437            gauge_neighbour_id=None
    435438
     439        #print gauge_neighbour_id
     440       
    436441        if interpolation_points is not None:
    437442            # Adjust for georef
  • anuga_core/source/anuga/fit_interpolate/interpolate.py

    r6187 r6189  
    3838from anuga.abstract_2d_finite_volumes.util import file_function
    3939from anuga.config import netcdf_mode_r, netcdf_mode_w, netcdf_mode_a
    40 from utilities.polygon import point_on_line
     40from utilities.polygon import interpolate_polyline
    4141
    4242
     
    507507       
    508508       
    509 
    510 ##
    511 # @brief Interpolate linearly from polyline nodes to midpoints of triangles.
    512 # @param f The data on the polyline nodes.
    513 # @param vertex_coordinates ??
    514 # @param gauge_neighbour_id ??
    515 # @param point_coordinates ??
    516 # @param verbose True if this function is to be verbose.
    517 def interpolate_polyline(f,
    518                          polyline_nodes,
    519                          gauge_neighbour_id,
    520                          point_coordinates=None,
    521                          verbose=False):
    522     """Interpolate linearly between values f on polyline nodes
    523     of a polyline to midpoints of triangles of boundary.
    524 
    525     f is the data on the polyline nodes.
    526 
    527     The mesh values representing a smooth surface are
    528     assumed to be specified in f.
    529 
    530     Inputs:
    531       f: Vector or array of data at the polyline nodes.
    532           If f is an array, interpolation will be done for each column as
    533           per underlying matrix-matrix multiplication
    534       point_coordinates: Interpolate polyline data to these positions.
    535           List of coordinate pairs [x, y] of
    536           data points or an nx2 Numeric array or a Geospatial_data object
    537          
    538     Output:
    539       Interpolated values at inputted points (z).
    540     """
    541     # FIXME: There is an option of passing a tolerance into this
    542    
    543     if isinstance(point_coordinates, Geospatial_data):
    544         point_coordinates = point_coordinates.get_data_points(absolute=True)
    545 
    546     z = num.zeros(len(point_coordinates), num.Float)
    547 
    548     f = ensure_numeric(f, num.Float)
    549     polyline_nodes = ensure_numeric(polyline_nodes, num.Float)
    550     point_coordinates = ensure_numeric(point_coordinates, num.Float)
    551     gauge_neighbour_id = ensure_numeric(gauge_neighbour_id, num.Int)
    552 
    553     n = polyline_nodes.shape[0] # Number of nodes in polyline       
    554     # Input sanity check
    555     msg = 'point coordinates are not given (interpolate.py)'
    556     assert point_coordinates is not None, msg
    557     msg = 'function value must be specified at every interpolation node'
    558     assert f.shape[0]==polyline_nodes.shape[0], msg
    559     msg = 'Must define function value at one or more nodes'
    560     assert f.shape[0]>0, msg
    561 
    562 
    563     if n == 1:
    564         msg = 'Polyline contained only one point. I need more. ' + str(f)
    565         raise Exception, msg
    566     elif n > 1:
    567         _interpolate_polyline_aux(point_coordinates, len(point_coordinates),
    568                                   polyline_nodes, len(polyline_nodes),
    569                                   f,
    570                                   gauge_neighbour_id,
    571                                   z)
    572        
    573     return z
    574 
    575        
    576 def _interpolate_polyline_aux(point_coordinates, number_of_points,
    577                               polyline_nodes, number_of_nodes,
    578                               f, gauge_neighbour_id, z):
    579     """Auxiliary function used by interpolate_polyline
    580     """
    581 
    582     for j in range(number_of_nodes):               
    583         neighbour_id = gauge_neighbour_id[j]
    584        
    585        
    586         if neighbour_id >= 0:
    587             x0, y0 = polyline_nodes[j,:]
    588             x1, y1 = polyline_nodes[neighbour_id,:]
    589            
    590             segment_len = sqrt((x1-x0)**2 + (y1-y0)**2)
    591             segment_delta = f[neighbour_id] - f[j]           
    592             slope = segment_delta/segment_len
    593            
    594                
    595             for i in range(number_of_points):               
    596                
    597                 x2, y2 = point_coordinates[i,:]
    598                 if point_on_line([x2, y2],
    599                                  [[x0, y0], [x1, y1]],
    600                                  rtol=1.0e-6):
    601                                  
    602 
    603                     dist = sqrt((x2-x0)**2 + (y2-y0)**2)
    604                     z[i] = slope*dist + f[j]
    605      
    606 
    607        
    608509       
    609510##
     
    1052953                                                      vertex_coordinates,
    1053954                                                      gauge_neighbour_id,
    1054                                                       point_coordinates=\
     955                                                      interpolation_points=\
    1055956                                                          self.interpolation_points)
    1056957                       
  • anuga_core/source/anuga/fit_interpolate/test_interpolate.py

    r6186 r6189  
    18461846        assert num.allclose(z, answer)
    18471847       
    1848        
    1849     def test_interpolate_polyline(self):
    1850         """test_interpolate_polyline(self):
    1851        
    1852         This test is added under the assumption that the function interpolate_polyline implemented by
    1853         John Jakeman works. It has been exercised somewhat by tests of sts boundary, but never before separately.
    1854         """
    1855        
    1856         f = num.array([58.06150614, 58.06150614, 58.06150614])
    1857         vertex_coordinates = num.array([[0., 0., ],
    1858                                         [4.04092634, 1106.11074699],
    1859                                         [8.08836552, 2212.16910609]])
    1860         gauge_neighbour_id = [1, 2, -1]
    1861         point_coordinates = num.array([[2.21870766e+03, 1.09802864e+03],
    1862                                        [1.62739645e+03, 2.20626983e+03],
    1863                                        [5.20084967e+02, 2.21030386e+03],
    1864                                        [6.06464546e+00, 1.65913993e+03],
    1865                                        [1.61934862e+03, -5.88143836e+00],
    1866                                        [5.11996623e+02, -1.85956061e+00],
    1867                                        [2.02046270e+00, 5.53055373e+02]])
    1868                              
    1869         z_ref = [0., 0., 0., 58.06150614, 0., 0., 58.06150614]
    1870        
    1871         z = interpolate_polyline(f, vertex_coordinates, gauge_neighbour_id, point_coordinates)
    1872         assert num.allclose(z, z_ref)
    1873        
    1874         # Another f
    1875         f = num.array([58.06150614, 158.06150614, 258.06150614])
    1876         z_ref = [0., 0., 0., 208.06150645, 0., 0., 108.0615061]       
    1877         z = interpolate_polyline(f, vertex_coordinates, gauge_neighbour_id, point_coordinates)
    1878         assert num.allclose(z, z_ref)       
    1879                        
    1880 
    1881         # Other and simpler numbers
    1882         f = num.array([1, 5, 13])       
    1883         vertex_coordinates = num.array([[0., 0., ],
    1884                                         [4., 4.],
    1885                                         [8., 8.]])       
    1886         point_coordinates = num.array([[0.1, 0.1],
    1887                                        [3.5, 3.5],
    1888                                        [4.0, 4.0],
    1889                                        [5.2, 5.2],
    1890                                        [7.0, 7.0],
    1891                                        [8.3, 8.3]])
    1892         gauge_neighbour_id = [1, 2, -1]
    1893                                                
    1894         z = interpolate_polyline(f, vertex_coordinates, gauge_neighbour_id, point_coordinates)
    1895         z_ref = [1.1, 4.5, 5., 7.4, 11., 0.]
    1896         #print z
    1897         assert num.allclose(z, z_ref)               
    1898        
    1899         # Test exception thrown for one point
    1900         f = num.array([5])               
    1901         vertex_coordinates = num.array([[4., 4.]])             
    1902         try: 
    1903             z = interpolate_polyline(f, vertex_coordinates, gauge_neighbour_id, point_coordinates)
    1904         except Exception:
    1905             pass
    1906         else:
    1907             raise Exception, 'One point should have raised exception'
    19081848                       
    19091849#-------------------------------------------------------------
  • anuga_core/source/anuga/utilities/polygon.py

    r6166 r6189  
    77from math import sqrt
    88from anuga.utilities.numerical_tools import ensure_numeric
    9 from anuga.geospatial_data.geospatial_data import ensure_absolute
     9from anuga.geospatial_data.geospatial_data import ensure_absolute, Geospatial_data
    1010
    1111
     
    957957
    958958    return reduced_polygon
     959   
     960   
     961       
     962
     963##
     964# @brief Interpolate linearly from polyline nodes to midpoints of triangles.
     965# @param data The data on the polyline nodes.
     966# @param polyline_nodes ??
     967# @param gauge_neighbour_id ??  FIXME(Ole): I want to get rid of this
     968# @param point_coordinates ??
     969# @param verbose True if this function is to be verbose.
     970def interpolate_polyline(data,
     971                         polyline_nodes,
     972                         gauge_neighbour_id,
     973                         interpolation_points=None,
     974                         rtol=1.0e-6,
     975                         atol=1.0e-8,
     976                         verbose=False):
     977    """Interpolate linearly between values data on polyline nodes
     978    of a polyline to list of interpolation points.
     979
     980    data is the data on the polyline nodes.
     981
     982
     983    Inputs:
     984      data: Vector or array of data at the polyline nodes.
     985      polyline_nodes: Location of nodes where data is available.     
     986      gauge_neighbour_id: ?
     987      interpolation_points: Interpolate polyline data to these positions.
     988          List of coordinate pairs [x, y] of
     989          data points or an nx2 Numeric array or a Geospatial_data object
     990      rtol, atol: Used to determine whether a point is on the polyline or not. See point_on_line.
     991
     992    Output:
     993      Interpolated values at interpolation points
     994    """
     995   
     996    if isinstance(interpolation_points, Geospatial_data):
     997        interpolation_points = interpolation_points.get_data_points(absolute=True)
     998
     999    interpolated_values = num.zeros(len(interpolation_points), num.Float)
     1000
     1001    data = ensure_numeric(data, num.Float)
     1002    polyline_nodes = ensure_numeric(polyline_nodes, num.Float)
     1003    interpolation_points = ensure_numeric(interpolation_points, num.Float)
     1004    gauge_neighbour_id = ensure_numeric(gauge_neighbour_id, num.Int)
     1005
     1006    n = polyline_nodes.shape[0] # Number of nodes in polyline       
     1007    # Input sanity check
     1008    msg = 'interpolation_points are not given (interpolate.py)'
     1009    assert interpolation_points is not None, msg
     1010    msg = 'function value must be specified at every interpolation node'
     1011    assert data.shape[0]==polyline_nodes.shape[0], msg
     1012    msg = 'Must define function value at one or more nodes'
     1013    assert data.shape[0]>0, msg
     1014
     1015    if n == 1:
     1016        msg = 'Polyline contained only one point. I need more. ' + str(data)
     1017        raise Exception, msg
     1018    elif n > 1:
     1019        _interpolate_polyline(data,
     1020                              polyline_nodes,
     1021                              gauge_neighbour_id,
     1022                              interpolation_points,                               
     1023                              interpolated_values,
     1024                              rtol,
     1025                              atol)
     1026       
     1027    return interpolated_values
     1028
     1029       
     1030def _interpolate_polyline(data,
     1031                          polyline_nodes,
     1032                          gauge_neighbour_id,
     1033                          interpolation_points,
     1034                          interpolated_values,
     1035                          rtol=1.0e-6,
     1036                          atol=1.0e-8):
     1037    """Auxiliary function used by interpolate_polyline
     1038   
     1039    NOTE: OBSOLETED BY C-EXTENSION
     1040    """
     1041   
     1042    number_of_nodes = len(polyline_nodes)               
     1043    number_of_points = len(interpolation_points)
     1044   
     1045    for j in range(number_of_nodes):               
     1046        neighbour_id = gauge_neighbour_id[j]
     1047       
     1048        # FIXME(Ole): I am convinced that gauge_neighbour_id can be discarded, but need to check with John J.
     1049        # Keep it for now (17 Jan 2009)
     1050        # When gone, we can simply interpolate between neighbouring nodes, i.e. neighbour_id = j+1.
     1051        # and the test below becomes something like: if j < number_of_nodes... 
     1052       
     1053        if neighbour_id >= 0:
     1054            x0, y0 = polyline_nodes[j,:]
     1055            x1, y1 = polyline_nodes[neighbour_id,:]
     1056           
     1057            segment_len = sqrt((x1-x0)**2 + (y1-y0)**2)
     1058            segment_delta = data[neighbour_id] - data[j]           
     1059            slope = segment_delta/segment_len
     1060           
     1061               
     1062            for i in range(number_of_points):               
     1063               
     1064                x, y = interpolation_points[i,:]
     1065                if point_on_line([x, y],
     1066                                 [[x0, y0], [x1, y1]],
     1067                                 rtol=rtol,
     1068                                 atol=atol):
     1069                                 
     1070
     1071                    dist = sqrt((x-x0)**2 + (y-y0)**2)
     1072                    interpolated_values[i] = slope*dist + data[j]
     1073     
     1074
     1075   
    9591076
    9601077##############################################
     
    9661083    from polygon_ext import _point_on_line
    9671084    from polygon_ext import _separate_points_by_polygon
     1085    from polygon_ext import _interpolate_polyline   
    9681086    #from polygon_ext import _intersection
    9691087
  • anuga_core/source/anuga/utilities/polygon_ext.c

    r5897 r6189  
    1818#include "math.h"
    1919
     20double dist(double x,
     21            double y) {
     22 
     23  return sqrt(x*x + y*y);
     24}
     25
    2026
    2127int __point_on_line(double x, double y,
     
    6672    // subject to specified absolute tolerance
    6773
    68     len_a = sqrt(a0*a0 + a1*a1);
    69     len_b = sqrt(b0*b0 + b1*b1);
     74    len_a = dist(a0, a1); //sqrt(a0*a0 + a1*a1);
     75    len_b = dist(b0, b1); //sqrt(b0*b0 + b1*b1);
    7076
    7177    if (a0*b0 + a1*b1 >= 0 && len_a <= len_b) {
     
    184190}
    185191*/
     192
     193
     194
     195int __interpolate_polyline(int number_of_nodes,
     196                           int number_of_points,
     197                           double* data,
     198                           double* polyline_nodes,
     199                           long* gauge_neighbour_id,
     200                           double* interpolation_points,                               
     201                           double* interpolated_values,
     202                           double rtol,
     203                           double atol) {
     204                           
     205  int j, i, neighbour_id;
     206  double x0, y0, x1, y1, x, y;
     207  double segment_len, segment_delta, slope, alpha;
     208
     209  for (j=0; j<number_of_nodes; j++) { 
     210
     211    neighbour_id = gauge_neighbour_id[j];
     212       
     213    // FIXME(Ole): I am convinced that gauge_neighbour_id can be discarded, but need to check with John J.
     214    // Keep it for now (17 Jan 2009)
     215    // When gone, we can simply interpolate between neighbouring nodes, i.e. neighbour_id = j+1.
     216    // and the test below becomes something like: if j < number_of_nodes... 
     217       
     218    if (neighbour_id >= 0) {
     219      x0 = polyline_nodes[2*j];
     220      y0 = polyline_nodes[2*j+1];
     221     
     222      x1 = polyline_nodes[2*neighbour_id];
     223      y1 = polyline_nodes[2*neighbour_id+1];     
     224     
     225           
     226      segment_len = dist(x1-x0, y1-y0);
     227      segment_delta = data[neighbour_id] - data[j];           
     228      slope = segment_delta/segment_len;
     229           
     230      for (i=0; i<number_of_points; i++) {               
     231        x = interpolation_points[2*i];
     232        y = interpolation_points[2*i+1];       
     233       
     234        if (__point_on_line(x, y, x0, y0, x1, y1, rtol, atol)) {
     235          alpha = dist(x-x0, y-y0);
     236          interpolated_values[i] = slope*alpha + data[j];
     237        }
     238      }
     239    }
     240  }
     241                           
     242  return 0;                         
     243}                                                     
    186244
    187245
     
    304362
    305363
     364
     365// Gateways to Python
     366PyObject *_interpolate_polyline(PyObject *self, PyObject *args) {
     367  //
     368  // _interpolate_polyline(data, polyline_nodes, gauge_neighbour_id, interpolation_points
     369  //                       interpolated_values):
     370  //
     371
     372 
     373  PyArrayObject
     374    *data,
     375    *polyline_nodes,
     376    *gauge_neighbour_id,
     377    *interpolation_points,
     378    *interpolated_values;
     379
     380  double rtol, atol; 
     381  int number_of_nodes, number_of_points, res;
     382 
     383  // Convert Python arguments to C
     384  if (!PyArg_ParseTuple(args, "OOOOOdd",
     385                        &data,
     386                        &polyline_nodes,
     387                        &gauge_neighbour_id,
     388                        &interpolation_points,
     389                        &interpolated_values,
     390                        &rtol,
     391                        &atol)) {
     392   
     393    PyErr_SetString(PyExc_RuntimeError,
     394                    "_interpolate_polyline could not parse input");
     395    return NULL;
     396  }
     397
     398  number_of_nodes = polyline_nodes -> dimensions[0];  // Number of nodes in polyline
     399  number_of_points = interpolation_points -> dimensions[0];   //Number of points
     400 
     401
     402  // Call underlying routine
     403  res = __interpolate_polyline(number_of_nodes,
     404                               number_of_points,
     405                               (double*) data -> data,
     406                               (double*) polyline_nodes -> data,
     407                               (long*) gauge_neighbour_id -> data,
     408                               (double*) interpolation_points -> data,                         
     409                               (double*) interpolated_values -> data,
     410                               rtol,
     411                               atol);                                                 
     412
     413  // Return
     414  return Py_BuildValue(""); 
     415}
     416
     417
     418
    306419/*
    307420PyObject *_intersection(PyObject *self, PyObject *args) {
     
    424537  {"_separate_points_by_polygon", _separate_points_by_polygon,
    425538                                 METH_VARARGS, "Print out"},
     539  {"_interpolate_polyline", _interpolate_polyline,
     540                                 METH_VARARGS, "Print out"},                             
    426541  {NULL, NULL, 0, NULL}   /* sentinel */
    427542};
  • anuga_core/source/anuga/utilities/test_polygon.py

    r6158 r6189  
    16901690
    16911691       
     1692    def test_interpolate_polyline(self):
     1693        """test_interpolate_polyline(self):
     1694       
     1695        This test is added under the assumption that the function interpolate_polyline implemented by
     1696        John Jakeman works. It has been exercised somewhat by tests of sts boundary, but never before separately.
     1697        """
     1698       
     1699        f = num.array([58.06150614, 58.06150614, 58.06150614])
     1700        vertex_coordinates = num.array([[0., 0., ],
     1701                                        [4.04092634, 1106.11074699],
     1702                                        [8.08836552, 2212.16910609]])
     1703        gauge_neighbour_id = [1, 2, -1]
     1704        point_coordinates = num.array([[2.21870766e+03, 1.09802864e+03],
     1705                                       [1.62739645e+03, 2.20626983e+03],
     1706                                       [5.20084967e+02, 2.21030386e+03],
     1707                                       [6.06464546e+00, 1.65913993e+03],
     1708                                       [1.61934862e+03, -5.88143836e+00],
     1709                                       [5.11996623e+02, -1.85956061e+00],
     1710                                       [2.02046270e+00, 5.53055373e+02]])
     1711                             
     1712        z_ref = [0., 0., 0., 58.06150614, 0., 0., 58.06150614]
     1713       
     1714        z = interpolate_polyline(f, vertex_coordinates, gauge_neighbour_id, point_coordinates)
     1715        assert num.allclose(z, z_ref)
     1716       
     1717        # Another f
     1718        f = num.array([58.06150614, 158.06150614, 258.06150614])
     1719        z_ref = [0., 0., 0., 208.06150645, 0., 0., 108.0615061]       
     1720        z = interpolate_polyline(f, vertex_coordinates, gauge_neighbour_id, point_coordinates)
     1721        assert num.allclose(z, z_ref)       
     1722                       
     1723
     1724        # Other and simpler numbers
     1725        f = num.array([1, 5, 13])       
     1726        vertex_coordinates = num.array([[0., 0., ],
     1727                                        [4., 4.],
     1728                                        [8., 8.]])       
     1729        point_coordinates = num.array([[0.1, 0.1],
     1730                                       [3.5, 3.5],
     1731                                       [4.0, 4.0],
     1732                                       [5.2, 5.2],
     1733                                       [7.0, 7.0],
     1734                                       [8.3, 8.3]])
     1735        gauge_neighbour_id = [1, 2, -1]
     1736                                               
     1737        z = interpolate_polyline(f, vertex_coordinates, gauge_neighbour_id, point_coordinates)
     1738        z_ref = [1.1, 4.5, 5., 7.4, 11., 0.]
     1739        #print z
     1740        assert num.allclose(z, z_ref)               
     1741       
     1742        # Test exception thrown for one point
     1743        f = num.array([5])               
     1744        vertex_coordinates = num.array([[4., 4.]])             
     1745        try: 
     1746            z = interpolate_polyline(f, vertex_coordinates, gauge_neighbour_id, point_coordinates)
     1747        except Exception:
     1748            pass
     1749        else:
     1750            raise Exception, 'One point should have raised exception'
     1751       
     1752
     1753        # More polyline nodes
     1754        data = num.array([1, 5, 13, 12, 6, 29])       
     1755        polyline_nodes = num.array([[0., 0.],
     1756                                    [4., 4.],
     1757                                    [8., 8.],
     1758                                    [10., 10.],
     1759                                    [10., 5.],
     1760                                    [10., 0.]])
     1761                                               
     1762        point_coordinates = num.array([[0.1, 0.1],
     1763                                       [3.5, 3.5],
     1764                                       [4.0, 4.0],
     1765                                       [5.2, 5.2],
     1766                                       [7.0, 7.0],
     1767                                       [8.3, 8.3],
     1768                                       [10., 10.],
     1769                                       [10., 9.],
     1770                                       [10., 7.1],
     1771                                       [10., 4.3],
     1772                                       [10., 1.0]])
     1773                                       
     1774        gauge_neighbour_id = [1, 2, 3, 4, 5, -1]
     1775                                               
     1776        z = interpolate_polyline(data, polyline_nodes, gauge_neighbour_id, point_coordinates)
     1777        z_ref = [1.1, 4.5, 5., 7.4, 11., 12.85, 12., 10.8, 8.52, 9.22, 24.4]
     1778        assert num.allclose(z, z_ref)               
     1779       
     1780       
     1781               
    16921782#-------------------------------------------------------------
    16931783if __name__ == "__main__":
Note: See TracChangeset for help on using the changeset viewer.