Algoritmi: Punto interno o esterno ad un poligono

#define ZERO     1e-8
#define PERTURB  1e-6

#define P_ABS(x) (((x) > 0) ? (x) : (-(x)))

/*
* This function determines, whether the 'node' lies inside or
* outside of the polygon, which is given in 'list'.
* There 'n' nodes in 'list'.
* The coordinates are given in 'Coords'.
*
* The function assumes, that there are no two nodes which are as
* close to each other as PERTURB !!!
*
* It uses the ray-tracing method. From the given node it draws a
* horizontal line to the right. It counts how many intersection
* occurs. If a node would lie exactly on the ray, then the node
* is moved up a little bit (perturbed). If the 'TO_BE_SURE_CHECK'
* macro is defined the function will check that the length of the
* boundary edges are bigger than the perturbation value 'PERTURB'.
* This check is for the insane. :-)
*
* Returns:  1 - node is inside
*           0 - node is outside
*          -1 - node is part of the list, better to chooose another node
*/

int Check2DNodeInside(double **Coords, int n, int *list,
                     int node, double x, double y)
{
 int    j;
 double a, b; /* parameters of the line which is formed by an edge */
 double ix;   /* x coord. of the intersection point of a segment and the ray */
 double px, py, cx, cy; /* segment nodes, previous, current */
 int    ci; /* counts intersections */
 
 /* start with the last node */
 px = Coords[list[n-1]-1][0];
 py = Coords[list[n-1]-1][1];

 /* nodes are identical */
 if(P_ABS(y-py) < ZERO && P_ABS(x-px) < ZERO)
   return 1;

 /* check whether the node lies on the trace exactly */
 if(P_ABS(y-py) < ZERO)
 {
   /* then we have to perturb it */
   py += PERTURB;
 }

 ci = 0;
 for(j = 0; j < n; j++)
 {
   /* if the nodes are identical, quit */
   if(node == list[j])
     break;

   cx = Coords[list[j]-1][0];
   cy = Coords[list[j]-1][1];

#ifdef TO_BE_SURE_CHECK
   {
     double ll;
     ll = sqrt((cx - px)*(cx - px) + (cy - py)*(cy - py));
     if(ll <= PERTURB)
       e_GeneralError("Check2DNodeInside: Boundary edge is smaller than PERTURB !");
   }
#endif

   /* nodes are identical */
   if(P_ABS(y-py) < ZERO && P_ABS(x-px) < ZERO)
     return 1;

   /* check whether current node lies on the trace exactly */
   if(P_ABS(y-cy)  < ZERO)
   {
     /* then we have to perturb it */
     cy += PERTURB;
   }
     
   /* vertical line ? */
   if(P_ABS(cx-px) > ZERO) /* no */
   {
     /* horizontal line ? */
     if(P_ABS(cy-py) > ZERO) /* no, general case */
     {
       a = (cy-py) / (cx-px);
       b = cy - a * cx;
       ix = (y - b) / a;
     }
     else /* yes */
     {
       /* side is horizontal there cannot be intersection then */
       ix = x-1;
     }
   }
   else /* yes */
   {
     ix = cx;
   }
     
   /* only nodes on the left hand side are counted */
   if(ix > x)
   {
     /* this is the general case */
     if(   ((py < y) && (y < cy))
        || ((cy < y) && (y < py)) )
     {
       ci++;
     }
   }
     
   px = cx;
   py = cy;
 }

 if(j == n)
 {
   return(ci % 2);
 }
 else
   return -1;
}