4
$\begingroup$

I have a problem that I have been attempting to solve for a few days now. I was wondering if I would be able to get some assistance from the community.

In order to detect if a point is in a polygon, I am using the winding number method. The algorithm is based on what I found on this site here: http://geomalgorithms.com/a03-_inclusion.html. In short, the algorithm will return a result that will indicate if the point is to the left or right of a line segment. For CCW orientation of the line, positive results mean that the point is to the left of the line and negative values indicate that the point is to the right of the line. The signs are switched for CW orientation. If the number of "lefts" is greater then the number of "rights", then the point is inside of the polygon. Lines whose endpoints are between the point are only considered in order to speed up the algorithm.

Now my polygons (I use this term loosely) are composed of both arcs and line segments. These polygons are simple and do not loop back on themselves. The polygons can either be convex or concave. As a side note, arcs in my program are represented as circular arcs.

The winding number algorithm works perfectly for line segments. But when we get to arcs, things can be a little bit more complicated. My current implementation is to consider any point that is underneath the curve of the arc would be inside of the arc. (Note: Source code for the implementation of my algorithm is given at the end of the post)

This implementation mostly works but, from the picture below, the red dot would be given as a false positive. The winding number algorithm would indicate that the red dot is inside the polygon when it clearly is not.

False Positive From the algorithm, line 1 will indicate that the point is to the left and line 4 will indicate that the point is to the right. But arc segment 2 will indicate that the point is also to the left and thus the point is inside of the polygon.

Now, the implementation with the arcs has another fault where the algorithm will return a false negative. Below is a picture of geometry where the point in the center would be outside of the polygon and the point off to the side is inside.

False Negative example

This polygon is composed of only arcs.

For a few days now, I have been considering different implementations but so far, nothing has panned out. Another thing about the arcs is that they are all in the CCW orientation. Whether the arc is concave/convex depends on which node is selected as the first node.

I am considering moving back to the ray crossing algorithm but I would like to see if I can get this to work first. I am fine with if I have 2 different algorithms for lines and arcs (in fact, I think that I will). I ideally, I would like to utilize the winding number algorithm but I am open to additional solutions.

Code:

I am posting a link to the source file for the code in case you want to see a more complete solution. I am also posting snippets that are relevant to convince. If you have any questions regarding the code, please do not hesitate to ask. Thank you

Point in polygon test: https://github.com/philm001/Omni-FEM/blob/master/Include/common/ClosedPath.h

bool pointInContour(wxRealPoint point)
    {
        int windingNumber = 0;
        bool isFirstRun = true;
        bool reverseWindingResult = false;

        double additionTerms = 0;
        double subtractionTerms = 0;

        // First, the algorithm will need to ensure that all of the lines are oriented in either CCW or CW.
        // Each contour, the order of the lines is CCW or CW but the order of the nodes is different.
        // For example, the start node of the 1st line could be "connected" to the end node of the next line.
        // In which case, we need to swap the nodes of the 1st line.

        for(auto lineIterator = p_closedPath.begin(); lineIterator != p_closedPath.end(); lineIterator++)
        {
            auto nextLineIterator = lineIterator + 1;

            if(nextLineIterator == p_closedPath.end())
                nextLineIterator = p_closedPath.begin();

            if(isFirstRun)
            {
                if(*(*lineIterator)->getSecondNode() != *(*nextLineIterator)->getFirstNode())
                {
                    if(*(*lineIterator)->getFirstNode() == *(*nextLineIterator)->getFirstNode())
                        (*lineIterator)->swap();
                    else if(*(*lineIterator)->getFirstNode() == *(*nextLineIterator)->getSecondNode())
                    {
                        (*lineIterator)->swap();
                        (*nextLineIterator)->swap();
                    }
                    else if(*(*lineIterator)->getSecondNode() == *(*nextLineIterator)->getSecondNode())
                        (*nextLineIterator)->swap();
                }

                isFirstRun = false;
            }
            else
            {
                if(*(*lineIterator)->getSecondNode() == *(*nextLineIterator)->getSecondNode())
                    (*nextLineIterator)->swap();
            }

            /* 
             * This is the actual shoelace algorithm that is implemented. In short, we have two terms, addition terms
             * and subtraction terms. THe addition terms are the summation of the Xpoint of the first node multiplied by
             * the Ypoint of the second node for all edges (arcs are a special case)
             * 
             * The subtraction term is the summation of the Xpoint of the second node multiplied by the
             * Ypoint of the first node (arcs are a special case)
             */ 
            if((*lineIterator)->isArc())
            {
                /* This algorithm is not being used to accurately determine the arc within a closed contour
                 * with that said, we only need to approximately determine the area. For arcs, we shall draw
                 * a line from the first node to the mid point and then from the mid point to the second node.
                 * We will use these two lines as the calculation for the shoelace algorithm.
                 */ 
                additionTerms += ((*lineIterator)->getFirstNode()->getCenter().x) * ((*lineIterator)->getMidPoint().y);
                subtractionTerms += ((*lineIterator)->getMidPoint().x) * ((*lineIterator)->getFirstNode()->getCenter().y);

                additionTerms += ((*lineIterator)->getMidPoint().x) * ((*lineIterator)->getSecondNode()->getCenter().y);
                subtractionTerms += ((*lineIterator)->getSecondNode()->getCenter().x) * ((*lineIterator)->getMidPoint().y);
            }
            else
            {
                additionTerms += ((*lineIterator)->getFirstNode()->getCenter().x) * ((*lineIterator)->getSecondNode()->getCenter().y);
                subtractionTerms += ((*lineIterator)->getSecondNode()->getCenter().x) * ((*lineIterator)->getFirstNode()->getCenter().y);
            }
        }

        // Next, we need to run the shoe-lace algorithm to determine if the ordering is CCW or CW. If CW, then we will need to swap
        // the start node and end node of all of the edges in order to ensure the polygon edge is in CCW
        double shoelaceResult = additionTerms - subtractionTerms;

        if(shoelaceResult < 0)
            reverseWindingResult = true;

        // Now we can perform the winding number algorithm
        for(auto lineIterator = p_closedPath.begin(); lineIterator != p_closedPath.end(); lineIterator++)
        {

            if((*lineIterator)->isArc() && (*lineIterator)->getSwappedState())
                (*lineIterator)->swap();

            double isLeftResult = (*lineIterator)->isLeft(point);

            if(!(*lineIterator)->isArc())
            {
                if(reverseWindingResult)
                    isLeftResult *= -1;

                if((*lineIterator)->getFirstNode()->getCenterYCoordinate() <= point.y)
                {
                    if((*lineIterator)->getSecondNode()->getCenterYCoordinate() > point.y)
                    {
                        if(isLeftResult > 0)
                            windingNumber++;
                        else if(isLeftResult < 0)
                            windingNumber--;
                    }
                }
                else if((*lineIterator)->getSecondNode()->getCenterYCoordinate() <= point.y)
                {
                    if(isLeftResult < 0)
                        windingNumber--;
                    else if(isLeftResult > 0)
                        windingNumber++;
                }

            }
            else
            {
                if((*lineIterator)->getSwappedState())
                {
                    isLeftResult *= -1;
                    // For arcs, we need to preserve the orientation of the first and second
                    // node for drawing purposes. Lines do not matter as much for drawing but arcs,
                    // it matters
                    (*lineIterator)->swap();
                }

                if(isLeftResult > 0)
                    windingNumber++;
                else
                    windingNumber--;
            }
        }

        if(windingNumber > 0)
            return true;
        else
            return false;
}

IsLeft function implantation (code begins on line 896): https://github.com/philm001/Omni-FEM/blob/master/Include/UI/geometryShapes.h

double isLeft(wxRealPoint point)
    {
        if(!p_isArc)
        {
            return ((_secondNode->getCenterXCoordinate() - _firstNode->getCenterXCoordinate()) * (point.y - _firstNode->getCenterYCoordinate()) 
                        - (point.x - _firstNode->getCenterXCoordinate()) * (_secondNode->getCenterYCoordinate() - _firstNode->getCenterYCoordinate())); 
        }
        else
        {
            double result = 0;

            if(isSameSign(crossProduct(point - _firstNode->getCenter(), _secondNode->getCenter() - _firstNode->getCenter()), 
                            crossProduct(this->getCenter() - _firstNode->getCenter(), _secondNode->getCenter() - _firstNode->getCenter())))
            {
                result = dotProduct(point - _firstNode->getCenter(), _secondNode->getCenter() - _firstNode->getCenter()) / 
                            dotProduct(_secondNode->getCenter() - _firstNode->getCenter(), _secondNode->getCenter() - _firstNode->getCenter());

                if(result >= 0 && result <= 1)
                    return 1.0;
                else
                    return -1.0;
            }
            else
            {
                result = dotProduct(point - this->getCenter(), point - this->getCenter());

                if(result <= pow(_radius, 2))
                    return 1.0;
                else
                    return -1.0;
            }


        }
    }
$\endgroup$
1
$\begingroup$

Enclose each one of your arcs by a box that connects the two end points on one diagonal and choose the other two points to make things into a box. Then replace the original polygon-with-arcs $P$ by the polygon $P'$ that contains all of the straight line segments plus, say, the inside two sides of all of the boxes. (Or the outsides, or choose randomly.)

Then, given a point to test, you can first figure out whether the point is inside one of the boxes. If it is not, then if it is inside (outside) $P'$ (easily testable) then it must also be inside (outside) $P$. On the other hand, if the point is in one of the boxes that enclose the arcs, then you'll have to test whether it is on one side of the arc or the other -- but that's easy because an arc has such a simple structure.

$\endgroup$
  • $\begingroup$ Hello Dr. Bangerth, Thank you for your response. To clarify some details, I do have a few questions regarding your solution 1) for the box that encloses the arc, how large should it be? I am imagining a box whose side is the line segment of the 2 points. The midpoint of the opposite side would intersect the midpoint of the arc itself. This way, the box is just big enough to fit the entire arc. Or, should I make the box slightly bigger? If this is the case, then I am confused on where the inside 2 sides of the box would be at? $\endgroup$ – philm Apr 2 '18 at 3:25
  • $\begingroup$ Additional, I assuming that the time that I would test if a point is enclosed by an arc is if the point lies in P' (sorry, I am not sure how i can format the text to have it look nicer) and in one of the arc-boxes? I understand that you are very busy but would it be possible to draw a very simple picture of the solution with all of the boxes (in fact, would you be able to show this on my 1st picture?) I have a feeling that I mostly understand what you are describing but I feel that a picture would clear up my thinking (I am more of a visualize person with these types of problems). $\endgroup$ – philm Apr 2 '18 at 3:30
  • $\begingroup$ Also, the polygon that encompasses the arc, does it need to specifically be a box or would a triangle be sufficient? $\endgroup$ – philm Apr 2 '18 at 15:21
  • $\begingroup$ It doesn't have to be a box -- what I'm suggesting is to enclose the arc by a polygon so that (i) the end points of the arc are vertices of the polygon, (ii) it is easy to test whether a point that's in this polygon is on one side of the arc or the other. $\endgroup$ – Wolfgang Bangerth Apr 2 '18 at 18:57
  • $\begingroup$ It's true that I don't have the time to draw a picture. But for each arc with that covers an angle of less than 90 degrees, you can take the triangle that connects the two end points and the other two sides are chosen as the tangent lines to the arc. For between 90 and 180 degrees, follow the tangent lines for a while and then draw a line between the two tangent lines in such a way that you don't intersect the arc. $\endgroup$ – Wolfgang Bangerth Apr 2 '18 at 19:07

Your Answer

By clicking "Post Your Answer", you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.