Rosetta Code

Ray-casting algorithm: Difference between revisions

Page Discussion

Ray-casting algorithm (view source)

Revision as of 04:00, 5 August 2011

583 bytes added ,  12 years ago
→‎{{header|C}}: replaced with much more rigorous code
(Updated D code)
(→‎{{header|C}}: replaced with much more rigorous code)
Line 328:
 
=={{header|C}}==
Required includes and definitions:
<lang c>#include <stdio.h>
#include <stdlib.h>
#include <stdboolstdarg.h>
#include <math.h>
 
typedef struct { double x, y; } vec;
typedef struct { int n; vec* v; } polygon_t, *polygon;
double x, y;
} point_t;
 
#define BIN_V(op, xx, yy) vec v##op(vec a,vec b){vec c;c.x=xx;c.y=yy;return c;}
typedef struct {
#define BIN_S(op, r) double v##op(vec a, vec b){ return r; }
point_t *vertices;
BIN_V(sub, a.x - b.x, a.y - b.y);
int edges[];
BIN_V(add, a.x + b.x, a.y + b.y);
} polygon_t;</lang>
BIN_S(dot, a.x * b.x + a.y * b.y);
BIN_S(cross, a.x * b.y - a.y * b.x);
 
/* return a + s * b */
Polygons for testing:
vec vmadd(vec a, double s, vec b)
{
vec c;
c.x = a.x + s * b.x;
c.y = a.y + s * b.y;
return c;
}
 
/* check if x0->x1 edge crosses y0->y1 edge. dx = x1 - x0, dy = y1 - y0, then
<lang c>point_t square_v[] =
solve x0 + a * dx == y0 + b * dy with a, b in real
{
cross both sides with dx, then: (remember, cross product is a scalar)
{0,0}, {10,0}, {10,10}, {0,10},
x0 X dx = y0 X dx + b * (dy X dx)
{2.5,2.5}, {7.5,0.1}, {7.5,7.5}, {2.5,7.5}
similarly,
};
x0 X dy + a * (dx X dy) == y0 X dy
there is an intersection iff 0 <= a <= 1 and 0 <= b <= 1
 
returns: 1 for intersect, -1 for not, 0 for hard to say (if the intersect
point_t esa_v[] =
point is too close to y0 or y1)
*/
int intersect(vec x0, vec x1, vec y0, vec y1, double tol, vec *sect)
{
vec dx = vsub(x1, x0), dy = vsub(y1, y0);
{3,0}, {7,0}, {10,5}, {7,10}, {3,10}, {0,5}
double d = vcross(dy, dx), a;
};
if (!d) return 0; /* edges are parallel */
 
a = (vcross(x0, dx) - vcross(y0, dx)) / d;
polygon_t esa =
if (sect)
{
*sect = vmadd(y0, a, dy);
esa_v,
{ 0,1, 1,2, 2,3, 3,4, 4,5, 5,0, -1 }
};
 
if (a < -tol || a > 1 + tol) return -1;
polygon_t square =
if (a < tol || a > 1 - tol) return 0;
{
square_v,
{ 0,1, 1,2, 2,3, 3,0, -1 }
};
 
a = (vcross(x0, dy) - vcross(y0, dy)) / d;
polygon_t squarehole =
if (a < 0 || a > 1) return -1;
{
square_v,
{ 0,1, 1,2, 2,3, 3,0,
4,5, 5,6, 6,7, 7,4, -1 }
};
 
return 1;
polygon_t strange =
}
 
/* distance between x and nearest point on y0->y1 segment. if the point
lies outside the segment, returns infinity */
double dist(vec x, vec y0, vec y1, double tol)
{
vec dy = vsub(y1, y0);
square_v,
vec x1, s;
{ 0,4, 4,3, 3,7, 7,6, 6,2, 2,1, 1,5, 5,0, -1 }
int r;
};</lang>
 
x1.x = x.x + dy.y; x1.y = x.y - dy.x;
Check for intersection code:
r = intersect(x, x1, y0, y1, tol, &s);
if (r == -1) return HUGE_VAL;
s = vsub(s, x);
return sqrt(vdot(s, s));
}
 
<lang c>#define MAXfor_v(Ai,B z, p) for(((A)i = 0, z = p->(B))?(A):(B)v; i < p->n; i++, z++)
/* returns 1 for inside, -1 for outside, 0 for on edge */
#define MIN(A,B) (((A)>(B))?(B):(A))
int inside(vec v, polygon p, double tol)
#define minP(A,B,C) ( (((A)->C) > ((B)->C)) ? (B) : (A) )
#define coeff_ang(PA,PB) ( ((PB)->y - (PA)->y) / ((PB)->x - (PA)->x) )
#define EPS 0.00001
inline bool hseg_intersect_seg(point_t *s, point_t *a, point_t *b)
{
/* should assert p->n > 1 */
double eps = 0.0;
int i, k, crosses;
if ( s->y == MAX(a->y, b->y) ||
vec *pv;
s->y == MIN(a->y, b->y) ) eps = EPS;
double min_x, max_x, min_y, max_y;
if ( (s->y + eps) > MAX(a->y, b->y) ||
(s->y + eps) < MIN(a->y, b->y) ||
s->x > MAX(a->x, b->x) ) return false;
if ( s->x <= MIN(a->x, b->x) ) return true;
double ca = (a->x != b->x) ? coeff_ang(a,b) :
HUGE_VAL;
point_t *my = minP(a,b,y);
double cp = (s->x - my->x) ? coeff_ang(my, s) :
HUGE_VAL;
return cp >= ca;
}</lang>
 
for (i = 0; i < p->n; i++) {
The ray-casting algorithm:
k = i + 1 % p->n;
min_x = dist(v, p->v[i], p->v[k], tol);
if (min_x < tol) return 0;
}
 
min_x = max_x = p->v[0].x;
<lang c>bool point_is_inside(polygon_t *poly, point_t *pt)
min_y = max_y = p->v[1].y;
{
int cross = 0, i;
for(i=0; poly->edges[i] != -1 ; i+=2) {
if ( hseg_intersect_seg(pt,
&poly->vertices[ poly->edges[i] ],
&poly->vertices[ poly->edges[i+1] ]) )
cross++;
}
return cross%2 != 0;
}</lang>
 
/* calculate extent of polygon */
Testing:
for_v(i, pv, p) {
if (pv->x > max_x) max_x = pv->x;
if (pv->x < min_x) min_x = pv->x;
if (pv->y > max_y) max_y = pv->y;
if (pv->y < min_y) min_y = pv->y;
}
if (v.x < min_x || v.x > max_x || v.y < min_y || v.y > max_y)
return -1;
 
max_x -= min_x; max_x *= 2;
<lang c>#define MAKE_TEST(S) do { \
max_y -= min_y; max_y *= 2;
printf("point (%.5f,%.5f) is ", test_points[i].x, test_points[i].y); \
max_x += max_y;
if ( point_is_inside(&S, &test_points[i]) ) \
 
printf("INSIDE " #S "\n"); \
vec e;
else \
while (1) {
printf("OUTSIDE " #S "\n"); \
crosses = 0;
} while(0);
/* pick a rand point far enough to be outside polygon */
e.x = v.x + (1 + rand() / (RAND_MAX + 1.)) * max_x;
e.y = v.y + (1 + rand() / (RAND_MAX + 1.)) * max_x;
 
for (i = 0; i < p->n; i++) {
k = (i + 1) % p->n;
k = intersect(v, e, p->v[i], p->v[k], tol, 0);
 
/* picked a bad point, ray got too close to vertex.
re-pick */
if (!k) break;
 
if (k == 1) crosses++;
}
if (i == p->n) break;
}
return (crosses & 1) ? 1 : -1;
}
 
int main()
{
int i;
point_t test_points[] = { {5,5}, {5,8}, {2,2},
vec vsq[] = { {0,0}, {10,100}, {2.510,2.510}, {0,10},
{02.015,2.5}, {27.25,70.41}, {07.5,7.5}, {10,2.5}, {-4,107.5}};
 
int i;
polygon_t sq = { 4, vsq }, /* outer square */
sq_hole = { 8, vsq }; /* outer and inner square, ie hole */
for(i=0; i < sizeof(test_points)/sizeof(point_t); i++) {
 
MAKE_TEST(square);
vec c = { 10, 5 }; /* on edge */
MAKE_TEST(squarehole);
vec d = { 5, 5 };
MAKE_TEST(strange);
 
MAKE_TEST(esa);
printf("%d\n", inside(c, &sq, 1e-10));
printf("%d\n", inside(c, &sq_hole, 1e-10));
}
 
printf("%d\n", inside(d, &sq, 1e-10)); /* in */
printf("%d\n", inside(d, &sq_hole, 1e-10)); /* out (in the hole) */
 
return EXIT_SUCCESS0;
}</lang>
 
The test's output reveals the meaning of <cite>coherent results</cite>: a point on the leftmost vertical side of the square (coordinate 0,5) is considered outside; while a point on the rightmost vertical side of the square (coordinate 10,5) is considered inside, but on the top-right vertex (coordinate 10,10), it is considered outside again.
=={{header|D}}==
<lang d>import std.stdio, std.math, std.algorithm, std.conv;
Anonymous user
Retrieved from "https://rosettacode.org/wiki/Special:MobileDiff/111307"