Re: point in polygon test

[Jakob Andreas Baerentzen]

On Mon, 23 Nov 1998, Novak Elliott wrote:

> does anyone have the source code for a point in polygon test ?

Funny - I am just concocting (in an incredible panic) a small ray tracing
exercise for graphics student ... The code below is point in triangle ...
but may be extended to point in polygon without (I think) too much
trouble. The basic strategy is to project the triangle on to an axis
parallel plane - the one that is most perpendicular to the triangle normal
and then see how many times a line from the intersection point to
y=infinite crosses the polygon. One crossing then we are inside. For the
more general case of a polygon I think the number just needs to be odd.

Hope it helps - I know the description is a bit jumbled. The code below
requires some vector stuff (C++) - and of course the plane intersection -  
but you probably have your own?

Anyway. Hope it helps

Andreas
-----------------


void Triangle::intersect(Ray& r) const {
  Vector pos; 
  double t; 
  if(plane.intersection(r,pos,t)) 
    {
      
      int axis1, axis2; 
      
      double n0 = fabs(normal[0]); 
      double n1 = fabs(normal[1]); 
      double n2 = fabs(normal[2]); 
        
      // Find the two axes most perpendicular to
      // the normal
      if (n0 > n1) 
        {
          axis1 = 1; 
          axis2 = (n0>n2)?2:0; 
        }
      else
        {
          axis1 = 0; 
          axis2 = (n1 > n2)?2:1; 
        }

      // Project the plane intersection point
      // and the three vertices of the
      // triangle onto the plane defined by these
      // two axes. 
      
      double Point2Dx = pos[axis1]; 
      double Point2Dy = pos[axis2]; 

      double v02Dx = v0[axis1]; 
      double v02Dy = v0[axis2]; 

      double v12Dx = v1[axis1]; 
      double v12Dy = v1[axis2]; 

      double v22Dx = v2[axis1]; 
      double v22Dy = v2[axis2]; 

      // Find the slope of the edges connecting the
      // three vertices. 
      double slope01 = (v12Dy - v02Dy) / (v12Dx - v02Dx); 
      double slope12 = (v22Dy - v12Dy) / (v22Dx - v12Dx); 
      double slope20 = (v02Dy - v22Dy) / (v02Dx - v22Dx); 

      int cuts = 0; 

      // For each of the three edges we test whether
      // the x coordinate of the intersection point lies
      // between the x coordinates of the two connected
      // vertices. 
      // If this is the case, we test whether the halfline
      // from the point and the direction of positive y intersects
      // the edge. 
      // If that is also the case, `cuts' is incremented. 
      // cuts == 1 means that we are inside the triangle
      // cuts == 2 or 0 means that we are outside. 
      // cuts == something else shouldn't occur. 

      if ( Point2Dx > d_min(v02Dx, v12Dx) &&
           Point2Dx < d_max(v02Dx, v12Dx) &&
           (v02Dy + (Point2Dx - v02Dx) * slope01) > Point2Dy) 
        cuts ++; 

      if ( Point2Dx > d_min(v12Dx, v22Dx) &&
           Point2Dx < d_max(v12Dx, v22Dx) &&
           (v12Dy + (Point2Dx - v12Dx) * slope12) > Point2Dy) 
        cuts ++; 

      if (Point2Dx > d_min(v02Dx, v22Dx) &&
          Point2Dx < d_max(v02Dx, v22Dx) &&
          (v22Dy + (Point2Dx - v22Dx) * slope20) > Point2Dy) 
        cuts ++; 

      if (cuts == 1) 
        r.cond_set_parameter(t,this); 

    } }