Anonymous user
Talk:Ray-casting algorithm: Difference between revisions
no edit summary
(→"Coherent results": new section) |
No edit summary |
||
Line 12:
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. --[[User:Ledrug|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>
| |||