点多边形算法

Question

我看到下面的算法用于检查点是否在此链接的给定多边形中:

int pnpoly(int nvert, float *vertx, float *verty, float testx, float testy)
{
  int i, j, c = 0;
  for (i = 0, j = nvert-1; i < nvert; j = i++) {
    if ( ((verty[i]>testy) != (verty[j]>testy)) &&
     (testx < (vertx[j]-vertx[i]) * (testy-verty[i]) / (verty[j]-verty[i]) + vertx[i]) )
       c = !c;
  }
  return c;
}

我试过这个算法,它实际上工作得很完美.但遗憾的是,在花了一些时间试图了解它之后,我无法理解它.

因此,如果有人能够理解这个算法,请向我解释一下.

谢谢.

Answer 1

该算法向右射线投射.循环的每次迭代,都会根据多边形的一个边缘检查测试点.如果点的y-coord在边缘范围内,则if测试的第一行成功.第二行检查测试点是否在线的左侧(我想 - 我没有任何废纸可供检查).如果确实如此,则从测试点向右绘制的线穿过该边缘.

通过反复反转该值c,算法计算向右线穿过多边形的次数.如果它经过了奇数次,则该点在内部; 如果是偶数,则该点在外面.

我会关注a)浮点运算的准确性,以及b)具有水平边缘的效果,或具有相同y坐标的测试点作为顶点的效果.

Answer 2

Chowlett在各方面,形状和形式上都是正确的.该算法假设如果您的点位于多边形的线上,那么这是在外面 - 在某些情况下,这是错误的.将两个'>'运算符更改为'> ='并将'<'更改为'<='将解决此问题.

bool PointInPolygon(Point point, Polygon polygon) {
  vector<Point> points = polygon.getPoints();
  int i, j, nvert = points.size();
  bool c = false;

  for(i = 0, j = nvert - 1; i < nvert; j = i++) {
    if( ( (points[i].y >= point.y ) != (points[j].y >= point.y) ) &&
        (point.x <= (points[j].x - points[i].x) * (point.y - points[i].y) / (points[j].y - points[i].y) + points[i].x)
      )
      c = !c;
  }

  return c;
}

Answer 3

这可能与解释实际代码中的光线跟踪算法一样详细.它可能没有被优化,但必须始终在完全掌握系统之后.

    //method to check if a Coordinate is located in a polygon
public boolean checkIsInPolygon(ArrayList<Coordinate> poly){
    //this method uses the ray tracing algorithm to determine if the point is in the polygon
    int nPoints=poly.size();
    int j=-999;
    int i=-999;
    boolean locatedInPolygon=false;
    for(i=0;i<(nPoints);i++){
        //repeat loop for all sets of points
        if(i==(nPoints-1)){
            //if i is the last vertex, let j be the first vertex
            j= 0;
        }else{
            //for all-else, let j=(i+1)th vertex
            j=i+1;
        }

        float vertY_i= (float)poly.get(i).getY();
        float vertX_i= (float)poly.get(i).getX();
        float vertY_j= (float)poly.get(j).getY();
        float vertX_j= (float)poly.get(j).getX();
        float testX  = (float)this.getX();
        float testY  = (float)this.getY();

        // following statement checks if testPoint.Y is below Y-coord of i-th vertex
        boolean belowLowY=vertY_i>testY;
        // following statement checks if testPoint.Y is below Y-coord of i+1-th vertex
        boolean belowHighY=vertY_j>testY;

        /* following statement is true if testPoint.Y satisfies either (only one is possible) 
        -->(i).Y < testPoint.Y < (i+1).Y        OR  
        -->(i).Y > testPoint.Y > (i+1).Y

        (Note)
        Both of the conditions indicate that a point is located within the edges of the Y-th coordinate
        of the (i)-th and the (i+1)- th vertices of the polygon. If neither of the above
        conditions is satisfied, then it is assured that a semi-infinite horizontal line draw 
        to the right from the testpoint will NOT cross the line that connects vertices i and i+1 
        of the polygon
        */
        boolean withinYsEdges= belowLowY != belowHighY;

        if( withinYsEdges){
            // this is the slope of the line that connects vertices i and i+1 of the polygon
            float slopeOfLine   = ( vertX_j-vertX_i )/ (vertY_j-vertY_i) ;

            // this looks up the x-coord of a point lying on the above line, given its y-coord
            float pointOnLine   = ( slopeOfLine* (testY - vertY_i) )+vertX_i;

            //checks to see if x-coord of testPoint is smaller than the point on the line with the same y-coord
            boolean isLeftToLine= testX < pointOnLine;

            if(isLeftToLine){
                //this statement changes true to false (and vice-versa)
                locatedInPolygon= !locatedInPolygon;
            }//end if (isLeftToLine)
        }//end if (withinYsEdges
    }

    return locatedInPolygon;
}

关于优化的一句话:最短(和/或最简洁)的代码实现得最快是不正确的.从数组中读取和存储元素并在代码块的执行中(可能)多次使用它比在每次需要时访问数组要快得多.如果阵列非常大,这尤其重要.在我看来,通过将一个数组的每个术语存储在一个命名良好的变量中,它也更容易评估其目的,从而形成一个更易读的代码.只是我的两分钱......

Answer 4

我更改了原始代码以使其更具可读性（这也使用了 Eigen）。算法是一样的。

// This uses the ray-casting algorithm to decide whether the point is inside
// the given polygon. See https://en.wikipedia.org/wiki/Point_in_polygon#Ray_casting_algorithm
bool pnpoly(const Eigen::MatrixX2d &poly, float x, float y)
{
    // If we never cross any lines we're inside.
    bool inside = false;

    // Loop through all the edges.
    for (int i = 0; i < poly.rows(); ++i)
    {
        // i is the index of the first vertex, j is the next one.
        // The original code uses a too-clever trick for this.
        int j = (i + 1) % poly.rows();

        // The vertices of the edge we are checking.
        double xp0 = poly(i, 0);
        double yp0 = poly(i, 1);
        double xp1 = poly(j, 0);
        double yp1 = poly(j, 1);

        // Check whether the edge intersects a line from (-inf,y) to (x,y).

        // First check if the line crosses the horizontal line at y in either direction.
        if ((yp0 <= y) && (yp1 > y) || (yp1 <= y) && (yp0 > y))
        {
            // If so, get the point where it crosses that line. This is a simple solution
            // to a linear equation. Note that we can't get a division by zero here -
            // if yp1 == yp0 then the above if be false.
            double cross = (xp1 - xp0) * (y - yp0) / (yp1 - yp0) + xp0;

            // Finally check if it crosses to the left of our test point. You could equally
            // do right and it should give the same result.
            if (cross < x)
                inside = !inside;
        }
    }
    return inside;
}