Changeset 2510


Ignore:
Timestamp:
Mar 9, 2006, 2:34:32 PM (18 years ago)
Author:
ole
Message:

Removed obsolete code

Location:
inundation/utilities
Files:
2 edited

Legend:

Unmodified
Added
Removed

  • inundation/utilities/polygon_ext.c

    r1910 r2510  
    269269
    270270
    271 
    272 //OBSOLETE
    273 
    274 int _inside_polygon(int M,           // Number of points
    275                     int N,           // Number of polygon vertices
    276                     double* points,
    277                     double* polygon,
    278                     long* indices,    // M-Array for storage indices
    279                     int closed,
    280                     int verbose) {
    281 
    282   double minpx, maxpx, minpy, maxpy, x, y, px_i, py_i, px_j, py_j;
    283   int i, j, k, count, inside;
    284 
    285   //Find min and max of poly used for optimisation when points
    286   //are far away from polygon
    287 
    288   minpx = polygon[0]; maxpx = minpx;
    289   minpy = polygon[1]; maxpy = minpy;
    290 
    291   for (i=0; i<N; i++) {
    292     px_i = polygon[2*i];
    293     py_i = polygon[2*i + 1];
    294 
    295     if (px_i < minpx) minpx = px_i;
    296     if (px_i > maxpx) maxpx = px_i;
    297     if (py_i < minpy) minpy = py_i;
    298     if (py_i > maxpy) maxpy = py_i;
    299   }
    300 
    301   //Begin main loop (for each point)
    302   count = 0;
    303   for (k=0; k<M; k++) {
    304     if (verbose){
    305       if (k %((M+10)/10)==0) printf("Doing %d of %d\n", k, M);
    306     }
    307    
    308     x = points[2*k];
    309     y = points[2*k + 1];
    310 
    311     inside = 0;
    312 
    313     //Optimisation
    314     if ((x > maxpx) || (x < minpx)) continue;
    315     if ((y > maxpy) || (y < minpy)) continue;
    316 
    317     //Check polygon
    318     for (i=0; i<N; i++) {
    319       //printf("k,i=%d,%d\n", k, i);
    320       j = (i+1)%N;
    321 
    322       px_i = polygon[2*i];
    323       py_i = polygon[2*i+1];
    324       px_j = polygon[2*j];
    325       py_j = polygon[2*j+1];
    326 
    327       //Check for case where point is contained in line segment
    328       if (_point_on_line(x, y, px_i, py_i, px_j, py_j)) {
    329         if (closed == 1) {
    330           inside = 1;
    331         } else {
    332           inside = 0;
    333         }
    334         break;
    335       } else {
    336         //Check if truly inside polygon
    337         if ( ((py_i < y) && (py_j >= y)) ||
    338              ((py_j < y) && (py_i >= y)) ) {
    339           if (px_i + (y-py_i)/(py_j-py_i)*(px_j-px_i) < x)
    340             inside = 1-inside;
    341         }
    342       }
    343     }
    344     if (inside == 1) {
    345       indices[count] = k;
    346       count++;
    347     }
    348   } // End k
    349 
    350   return count;
    351 }
    352 
    353 
    354 
    355 PyObject *inside_polygon(PyObject *self, PyObject *args) {
    356   //def inside_polygon(point, polygon, closed, verbose, one_point):
    357   //  """Determine whether points are inside or outside a polygon
    358   //
    359   //  Input:
    360   //     point - Tuple of (x, y) coordinates, or list of tuples
    361   //     polygon - list of vertices of polygon
    362   //     one_poin - If True Boolean value should be returned
    363   //                If False, indices of points inside returned
    364   //     closed - (optional) determine whether points on boundary should be
    365   //     regarded as belonging to the polygon (closed = True)
    366   //     or not (closed = False)
    367 
    368   //
    369   //  Output:
    370   //     If one point is considered, True or False is returned.
    371   //     If multiple points are passed in, the function returns indices
    372   //     of those points that fall inside the polygon
    373   //
    374   //  Examples:
    375   //     inside_polygon( [0.5, 0.5], [[0,0], [1,0], [1,1], [0,1]] )
    376   //     will evaluate to True as the point 0.5, 0.5 is inside the unit square
    377   //
    378   //     inside_polygon( [[0.5, 0.5], [1, -0.5], [0.3, 0.2]] )
    379   //     will return the indices [0, 2] as only the first and the last point
    380   //     is inside the unit square
    381   //
    382   //  Remarks:
    383   //     The vertices may be listed clockwise or counterclockwise and
    384   //     the first point may optionally be repeated.
    385   //     Polygons do not need to be convex.
    386   //     Polygons can have holes in them and points inside a hole is
    387   //     regarded as being outside the polygon.
    388   //
    389   //
    390   //  Algorithm is based on work by Darel Finley,
    391   //  http://www.alienryderflex.com/polygon/
    392   //
    393   //
    394 
    395   PyArrayObject
    396     *point,
    397     *polygon,
    398     *indices;
    399 
    400   int closed, verbose; //Flags
    401   int count, M, N;
    402 
    403   // Convert Python arguments to C
    404   if (!PyArg_ParseTuple(args, "OOOii",
    405                         &point,
    406                         &polygon,
    407                         &indices,
    408                         &closed,
    409                         &verbose))
    410     return NULL;
    411 
    412   M = point -> dimensions[0];    //Number of points
    413   N = polygon -> dimensions[0];  //Number of vertices in polygon
    414 
    415   if (verbose) printf("Got %d points and %d polygon vertices\n", M, N);
    416  
    417   //Call underlying routine
    418   count = _inside_polygon(M, N,
    419                           (double*) point -> data,
    420                           (double*) polygon -> data,
    421                           (long*) indices -> data,
    422                           closed, verbose);
    423 
    424   //NOTE: return number of points inside..
    425   //printf("COunt=%d\n", count);
    426   return Py_BuildValue("i", count);
    427 }

  • inundation/utilities/util_ext.c

    r2508 r2510  
    2121#include "util_ext.h"
    2222
    23 
    24 
    25 int _point_on_line(double x, double y,
    26                    double x0, double y0,
    27                    double x1, double y1) {
    28   /*Determine whether a point is on a line segment
    29 
    30     Input: x, y, x0, x0, x1, y1: where
    31         point is given by x, y
    32         line is given by (x0, y0) and (x1, y1)
    33 
    34   */
    35 
    36   double a0, a1, a_normal0, a_normal1, b0, b1, len_a, len_b;
    37 
    38   a0 = x - x0;
    39   a1 = y - y0;
    40 
    41   a_normal0 = a1;
    42   a_normal1 = -a0;
    43 
    44   b0 = x1 - x0;
    45   b1 = y1 - y0;
    46 
    47   if ( a_normal0*b0 + a_normal1*b1 == 0 ) {
    48     //Point is somewhere on the infinite extension of the line
    49 
    50     len_a = sqrt(a0*a0 + a1*a1);
    51     len_b = sqrt(b0*b0 + b1*b1);
    52 
    53     if (a0*b0 + a1*b1 >= 0 && len_a <= len_b) {
    54       return 1;
    55     } else {
    56       return 0;
    57     }
    58   } else {
    59     return 0;
    60   }
    61 }
    62 
    63 
    64 
    65 int _separate_points_by_polygon(int M,     // Number of points
    66                                 int N,     // Number of polygon vertices
    67                                 double* points,
    68                                 double* polygon,
    69                                 long* indices,  // M-Array for storage indices
    70                                 int closed,
    71                                 int verbose) {
    72 
    73   double minpx, maxpx, minpy, maxpy, x, y, px_i, py_i, px_j, py_j;
    74   int i, j, k, outside_index, inside_index, inside;
    75 
    76   //Find min and max of poly used for optimisation when points
    77   //are far away from polygon
    78 
    79   minpx = polygon[0]; maxpx = minpx;
    80   minpy = polygon[1]; maxpy = minpy;
    81 
    82   for (i=0; i<N; i++) {
    83     px_i = polygon[2*i];
    84     py_i = polygon[2*i + 1];
    85 
    86     if (px_i < minpx) minpx = px_i;
    87     if (px_i > maxpx) maxpx = px_i;
    88     if (py_i < minpy) minpy = py_i;
    89     if (py_i > maxpy) maxpy = py_i;
    90   }
    91 
    92   //Begin main loop (for each point)
    93   inside_index = 0;    //Keep track of points inside
    94   outside_index = M-1; //Keep track of points outside (starting from end)   
    95   for (k=0; k<M; k++) {
    96     if (verbose){
    97       if (k %((M+10)/10)==0) printf("Doing %d of %d\n", k, M);
    98     }
    99    
    100     x = points[2*k];
    101     y = points[2*k + 1];
    102 
    103     inside = 0;
    104 
    105     //Optimisation
    106     if ((x > maxpx) || (x < minpx) || (y > maxpy) || (y < minpy)) {
    107       //Nothing
    108     } else {   
    109       //Check polygon
    110       for (i=0; i<N; i++) {
    111         //printf("k,i=%d,%d\n", k, i);
    112         j = (i+1)%N;
    113 
    114         px_i = polygon[2*i];
    115         py_i = polygon[2*i+1];
    116         px_j = polygon[2*j];
    117         py_j = polygon[2*j+1];
    118 
    119         //Check for case where point is contained in line segment
    120         if (_point_on_line(x, y, px_i, py_i, px_j, py_j)) {
    121           if (closed == 1) {
    122             inside = 1;
    123           } else {
    124             inside = 0;
    125           }
    126           break;
    127         } else {
    128           //Check if truly inside polygon
    129           if ( ((py_i < y) && (py_j >= y)) ||
    130                ((py_j < y) && (py_i >= y)) ) {
    131             if (px_i + (y-py_i)/(py_j-py_i)*(px_j-px_i) < x)
    132               inside = 1-inside;
    133           }
    134         }
    135       }
    136     }
    137     if (inside == 1) {
    138       indices[inside_index] = k;
    139       inside_index += 1;
    140     } else {
    141       indices[outside_index] = k;
    142       outside_index -= 1;   
    143     }
    144   } // End k
    145 
    146   return inside_index;
    147 }
    148 
    149 
    150 
    151 // Gateways to Python
    152 PyObject *point_on_line(PyObject *self, PyObject *args) {
    153   //
    154   // point_on_line(x, y, x0, y0, x1, y1)
    155   //
    156 
    157   double x, y, x0, y0, x1, y1;
    158   int res;
    159   PyObject *result;
    160 
    161   // Convert Python arguments to C
    162   if (!PyArg_ParseTuple(args, "dddddd", &x, &y, &x0, &y0, &x1, &y1))
    163     return NULL;
    164 
    165 
    166   // Call underlying routine
    167   res = _point_on_line(x, y, x0, y0, x1, y1);
    168 
    169   // Return values a and b
    170   result = Py_BuildValue("i", res);
    171   return result;
    172 }
    173 
    174 
    175 
    176 PyObject *separate_points_by_polygon(PyObject *self, PyObject *args) {
    177   //def separate_points_by_polygon(points, polygon, closed, verbose, one_point):
    178   //  """Determine whether points are inside or outside a polygon
    179   //
    180   //  Input:
    181   //     point - Tuple of (x, y) coordinates, or list of tuples
    182   //     polygon - list of vertices of polygon
    183   //     closed - (optional) determine whether points on boundary should be
    184   //     regarded as belonging to the polygon (closed = True)
    185   //     or not (closed = False)
    186 
    187   //
    188   //  Output:
    189   //     indices: array of same length as points with indices of points falling
    190   //     inside the polygon listed from the beginning and indices of points
    191   //     falling outside listed from the end.
    192   //     
    193   //     count: count of points falling inside the polygon
    194   //     
    195   //     The indices of points inside are obtained as indices[:count]
    196   //     The indices of points outside are obtained as indices[count:]       
    197   //
    198   //  Examples:
    199   //     separate_polygon( [[0.5, 0.5], [1, -0.5], [0.3, 0.2]] )
    200   //     will return the indices [0, 2, 1] as only the first and the last point
    201   //     is inside the unit square
    202   //
    203   //  Remarks:
    204   //     The vertices may be listed clockwise or counterclockwise and
    205   //     the first point may optionally be repeated.
    206   //     Polygons do not need to be convex.
    207   //     Polygons can have holes in them and points inside a hole is
    208   //     regarded as being outside the polygon.
    209   //
    210   //
    211   //  Algorithm is based on work by Darel Finley,
    212   //  http://www.alienryderflex.com/polygon/
    213   //
    214   //
    215 
    216   PyArrayObject
    217     *points,
    218     *polygon,
    219     *indices;
    220 
    221   int closed, verbose; //Flags
    222   int count, M, N;
    223 
    224   // Convert Python arguments to C
    225   if (!PyArg_ParseTuple(args, "OOOii",
    226                         &points,
    227                         &polygon,
    228                         &indices,
    229                         &closed,
    230                         &verbose))
    231     return NULL;
    232 
    233   M = points -> dimensions[0];   //Number of points
    234   N = polygon -> dimensions[0];  //Number of vertices in polygon
    235 
    236   if (verbose) printf("Got %d points and %d polygon vertices\n", M, N);
    237  
    238   //Call underlying routine
    239   count = _separate_points_by_polygon(M, N,
    240                                       (double*) points -> data,
    241                                       (double*) polygon -> data,
    242                                       (long*) indices -> data,
    243                                       closed, verbose);
    244  
    245   //NOTE: return number of points inside..
    246   return Py_BuildValue("i", count);
    247 }
    24823
    24924
    … …  
    30378  {"gradient", gradient, METH_VARARGS, "Print out"},
    30479  {"gradient2", gradient2, METH_VARARGS, "Print out"}, 
    305   {"point_on_line", point_on_line, METH_VARARGS, "Print out"},
    306   //{"inside_polygon", inside_polygon, METH_VARARGS, "Print out"},
    307   {"separate_points_by_polygon", separate_points_by_polygon,
    308                                  METH_VARARGS, "Print out"},
    30980  {NULL, NULL, 0, NULL}   /* sentinel */
    31081};
    … …  
    32293
    32394
    324 //OBSOLETE
    325 
    326 int _inside_polygon(int M,           // Number of points
    327                     int N,           // Number of polygon vertices
    328                     double* points,
    329                     double* polygon,
    330                     long* indices,    // M-Array for storage indices
    331                     int closed,
    332                     int verbose) {
    333 
    334   double minpx, maxpx, minpy, maxpy, x, y, px_i, py_i, px_j, py_j;
    335   int i, j, k, count, inside;
    336 
    337   //Find min and max of poly used for optimisation when points
    338   //are far away from polygon
    339 
    340   minpx = polygon[0]; maxpx = minpx;
    341   minpy = polygon[1]; maxpy = minpy;
    342 
    343   for (i=0; i<N; i++) {
    344     px_i = polygon[2*i];
    345     py_i = polygon[2*i + 1];
    346 
    347     if (px_i < minpx) minpx = px_i;
    348     if (px_i > maxpx) maxpx = px_i;
    349     if (py_i < minpy) minpy = py_i;
    350     if (py_i > maxpy) maxpy = py_i;
    351   }
    352 
    353   //Begin main loop (for each point)
    354   count = 0;
    355   for (k=0; k<M; k++) {
    356     if (verbose){
    357       if (k %((M+10)/10)==0) printf("Doing %d of %d\n", k, M);
    358     }
    359    
    360     x = points[2*k];
    361     y = points[2*k + 1];
    362 
    363     inside = 0;
    364 
    365     //Optimisation
    366     if ((x > maxpx) || (x < minpx)) continue;
    367     if ((y > maxpy) || (y < minpy)) continue;
    368 
    369     //Check polygon
    370     for (i=0; i<N; i++) {
    371       //printf("k,i=%d,%d\n", k, i);
    372       j = (i+1)%N;
    373 
    374       px_i = polygon[2*i];
    375       py_i = polygon[2*i+1];
    376       px_j = polygon[2*j];
    377       py_j = polygon[2*j+1];
    378 
    379       //Check for case where point is contained in line segment
    380       if (_point_on_line(x, y, px_i, py_i, px_j, py_j)) {
    381         if (closed == 1) {
    382           inside = 1;
    383         } else {
    384           inside = 0;
    385         }
    386         break;
    387       } else {
    388         //Check if truly inside polygon
    389         if ( ((py_i < y) && (py_j >= y)) ||
    390              ((py_j < y) && (py_i >= y)) ) {
    391           if (px_i + (y-py_i)/(py_j-py_i)*(px_j-px_i) < x)
    392             inside = 1-inside;
    393         }
    394       }
    395     }
    396     if (inside == 1) {
    397       indices[count] = k;
    398       count++;
    399     }
    400   } // End k
    401 
    402   return count;
    403 }
    404 
    405 
    406 
    407 PyObject *inside_polygon(PyObject *self, PyObject *args) {
    408   //def inside_polygon(point, polygon, closed, verbose, one_point):
    409   //  """Determine whether points are inside or outside a polygon
    410   //
    411   //  Input:
    412   //     point - Tuple of (x, y) coordinates, or list of tuples
    413   //     polygon - list of vertices of polygon
    414   //     one_poin - If True Boolean value should be returned
    415   //                If False, indices of points inside returned
    416   //     closed - (optional) determine whether points on boundary should be
    417   //     regarded as belonging to the polygon (closed = True)
    418   //     or not (closed = False)
    419 
    420   //
    421   //  Output:
    422   //     If one point is considered, True or False is returned.
    423   //     If multiple points are passed in, the function returns indices
    424   //     of those points that fall inside the polygon
    425   //
    426   //  Examples:
    427   //     inside_polygon( [0.5, 0.5], [[0,0], [1,0], [1,1], [0,1]] )
    428   //     will evaluate to True as the point 0.5, 0.5 is inside the unit square
    429   //
    430   //     inside_polygon( [[0.5, 0.5], [1, -0.5], [0.3, 0.2]] )
    431   //     will return the indices [0, 2] as only the first and the last point
    432   //     is inside the unit square
    433   //
    434   //  Remarks:
    435   //     The vertices may be listed clockwise or counterclockwise and
    436   //     the first point may optionally be repeated.
    437   //     Polygons do not need to be convex.
    438   //     Polygons can have holes in them and points inside a hole is
    439   //     regarded as being outside the polygon.
    440   //
    441   //
    442   //  Algorithm is based on work by Darel Finley,
    443   //  http://www.alienryderflex.com/polygon/
    444   //
    445   //
    446 
    447   PyArrayObject
    448     *point,
    449     *polygon,
    450     *indices;
    451 
    452   int closed, verbose; //Flags
    453   int count, M, N;
    454 
    455   // Convert Python arguments to C
    456   if (!PyArg_ParseTuple(args, "OOOii",
    457                         &point,
    458                         &polygon,
    459                         &indices,
    460                         &closed,
    461                         &verbose))
    462     return NULL;
    463 
    464   M = point -> dimensions[0];    //Number of points
    465   N = polygon -> dimensions[0];  //Number of vertices in polygon
    466 
    467   if (verbose) printf("Got %d points and %d polygon vertices\n", M, N);
    468  
    469   //Call underlying routine
    470   count = _inside_polygon(M, N,
    471                           (double*) point -> data,
    472                           (double*) polygon -> data,
    473                           (long*) indices -> data,
    474                           closed, verbose);
    475 
    476   //NOTE: return number of points inside..
    477   //printf("COunt=%d\n", count);
    478   return Py_BuildValue("i", count);
    479 }
Note: See TracChangeset for help on using the changeset viewer.