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.
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.
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;
}
}
}
