Talk:Ray-casting algorithm
Fortran: The point of pts?
What is pts used for in the Fortran example? --Paddy3118 01:17, 26 May 2009 (UTC)
- Ah, It seems there are one set of points for all polygons. --Paddy3118 01:41, 26 May 2009 (UTC)
- Yes, I preferred to use a single set, and to index it (it's not an odd technics if a subset of points is shared among several polygons we want to define). --ShinTakezou 16:53, 26 May 2009 (UTC)
AutoHotkey version
The AutoHotkey solution refers to ray_intersects_segment without defining it, but it seems to me that defining ray_intersects_segment is a key part of this task. --Mr2001 07:31, 31 October 2010 (UTC)
- Help:ENA Templates will be helpful for flagging this kind of thing. --Michael Mol 13:23, 31 October 2010 (UTC)
"Coherent results"
The C example claims to "reveal the meaning of coherent results", whereby it considers a point lying on the left edge of a square outside, but a point on right side inside. This must be a meaning of "coherent" that I wasn't aware of. Also, the intersection code is terrible: it uses an arbitrarily defined value as a floating point precision limit while not rigorously enforcing it, which will cause all kinds of unexpected behaviors when used in real work. --Ledrug 22:10, 25 July 2011 (UTC)
D code generate false positive
We port the D code in our C++ project and we found false positive, I almost sure all other implementations could have exactly the same issue.
Here is my modifications of the D code showing the issue (show comments) : <lang d> import std.stdio, std.math, std.algorithm;
immutable struct Point { float x, y; } // Using float instead of double (reducing precision type) immutable struct Edge { Point a, b; } immutable struct Figure {
string name; Edge[] edges;
}
bool contains(in Figure poly, in Point p) pure nothrow @safe {
static bool raySegI(in Point p, in Edge edge)
pure nothrow @safe { enum float epsilon = 0.001; // Increase error tolerance with (edge) { if (a.y > b.y) //swap(a, b); // if edge is mutable return raySegI(p, Edge(b, a)); if (p.y == a.y || p.y == b.y) //p.y += epsilon; // if p is mutable return raySegI(Point(p.x, p.y + epsilon), edge); if (p.y > b.y || p.y < a.y || p.x > max(a.x, b.x)) return false; if (p.x < min(a.x, b.x)) return true; immutable blue = (abs(a.x - p.x) > float.min_normal) ? ((p.y - a.y) / (p.x - a.x)) : float.max; immutable red = (abs(a.x - b.x) > float.min_normal) ? ((b.y - a.y) / (b.x - a.x)) : float.max; return blue >= red; } }
return poly.edges.count!(e => raySegI(p, e)) % 2;
}
void main() {
immutable Figure[] polys = [
{"Square", [ {{ 0.0, 0.0}, {10.0, 0.0}}, {{10.0, 0.0}, {10.0, 10.0}}, {{10.0, 10.0}, { 0.0, 10.0}}, {{ 0.0, 10.0}, { 0.0, 0.0}}]}, {"Square hole", [ {{ 0.0, 0.0}, {10.0, 0.0}}, {{10.0, 0.0}, {10.0, 10.0}}, {{10.0, 10.0}, { 0.0, 10.0}}, {{ 0.0, 10.0}, { 0.0, 0.0}}, {{ 2.5, 2.5}, { 7.5, 2.5}}, {{ 7.5, 2.5}, { 7.5, 7.5}}, {{ 7.5, 7.5}, { 2.5, 7.5}}, {{ 2.5, 7.5}, { 2.5, 2.5}}]}, {"Strange", [ {{ 0.0, 0.0}, { 2.5, 2.5}}, {{ 2.5, 2.5}, { 0.0, 10.0}}, {{ 0.0, 10.0}, { 2.5, 7.5}}, {{ 2.5, 7.5}, { 7.5, 7.5}}, {{ 7.5, 7.5}, {10.0, 10.0}}, {{10.0, 10.0}, {10.0, 0.0}}, {{10.0, 0}, { 2.5, 2.5}}]}, {"Exagon", [ {{ 3.0, 0.0}, { 7.0, 0.0}}, {{ 7.0, 0.0}, {10.0, 5.0}}, {{10.0, 5.0}, { 7.0, 10.0}}, {{ 7.0, 10.0}, { 3.0, 10.0}}, {{ 3.0, 10.0}, { 0.0, 5.0}}, {{ 0.0, 5.0}, { 3.0, 0.0}}]}, {"Strange2", [ // Our testing polygon {{-0.047600,3.619000},{-0.147595,3.618007}}, {{-0.147595,3.618007},{-0.111895,0.022807}}, {{-0.111895,0.022807},{-0.011900,0.023800}}, {{-0.011900,0.023800},{0.088095,0.024793}}, {{0.088095,0.024793},{0.052395,3.619993}}, {{0.052395,3.619993},{-0.047600,3.619000}}]} ];
immutable Point[] testPoints = [{ 5, 5}, {5, 8}, {-10, 5}, {0, 5}, {10, 5}, {8, 5}, { 10, 10}, {-5.0,0.0238}]; // Last point need to be out of Strange2 polygon, but the result is true
foreach (immutable poly; polys)
{
writefln(`Is point inside figure "%s"?`, poly.name);
foreach (immutable p; testPoints)
writefln(" (%3s, %2s): %s", p.x, p.y, contains(poly, p));
writeln;
}
readln(); } </lang>
To see our polygon and test point you can use geogebra web :
* go to http://www.geogebra.org/webstart/geogebra.html * enter in the input field : (-5.000000f,0.023800f) * enter in the input field : Polygon[(-0.047600,3.619000),(-0.147595,3.618007),(-0.111895,0.022807),(-0.011900,0.023800),(0.088095,0.024793),(0.052395,3.619993)]
We fix this issue by shifting the y coordinate of input point (adding epsilon) while it match with a vertex.
Notice we are shifting the point before testing it against any segment, this increase the stability in the computation of the number of intersections.
Here is our working C++ implementation : <lang c++>
template<class T, class Alloc>
bool Polygon2<T, Alloc>::contains(const Vector2<T>& point, T epsilon) const
{
H3D_ASSERT(this->size() > 1);
Vector2<T> currPt;
// Insure point is not equal to one of the vertex
Vector2<T> shiftedPoint(point);
for (auto it = this->begin(); it != this->end();)
{
currPt = *it;
if (math::epsilonEquals(currPt.y, shiftedPoint.y, epsilon))
{
shiftedPoint.y += epsilon;
it = this->begin(); // Shift the point and recheck all the vertex (we could have shifted the point on a previous vertex)
}
else
{
++it;
}
}
int count = 0;
Vector2<T> lastPt = this->back();
for (auto it = this->begin(); it != this->end(); ++it)
{
currPt = *it;
if (raySegI(shiftedPoint, lastPt, currPt, epsilon))
++count;
lastPt = currPt;
}
return (bool)(count % 2);
}
template<class T, class Alloc>
bool Polygon2<T, Alloc>::raySegI(const Vector2<T>& p, const Vector2<T>& a, const Vector2<T>& b, T epsilon) const
{
if (a.y > b.y)
return raySegI(p, b, a, epsilon);
H3D_ASSERT(math::epsilonEquals(p.y, a.y, epsilon) == false && math::epsilonEquals(p.y, b.y, epsilon) == false);
if (p.y > b.y || p.y < a.y || p.x > std::max(a.x, b.x))
return false;
if (p.x < std::min(a.x, b.x))
return true;
T blue = (abs(a.x - p.x) > std::numeric_limits<T>::min()) ?
((p.y - a.y) / (p.x - a.x)) :
std::numeric_limits<T>::max();
T red = (abs(a.x - b.x) > std::numeric_limits<T>::min()) ?
((b.y - a.y) / (b.x - a.x)) :
std::numeric_limits<T>::max();
return blue >= red;
}
</lang>