Talk:Ray-casting algorithm: Difference between revisions
Content added Content deleted
(→"Coherent results": new section) |
No edit summary |
||
Line 12: | 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) |
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> |