buggy intersectwith() implementations?
buggy intersectwith() implementations?
hello to everyone. i recently tried to find the intersection between a polyline and a boundless line perpendiculal to the x-axis in a known point, in order to compute its distance from the first vertex along the polyline and its y-coordinate. no way to get it.
i first tried declaring an odgeline3d and searching for intersections with each odgelineseg3d belonging to the polyline (almost rude, i guess). the y-coordinate of the retrieved point was completely out of my expectations and i didn't even try to get the distance from the start.
then i tried to replace odgeline3d with odgelineseg3d; to calculate the intersection of the line with the proper segment of the polyline and, vice versa, of the segment with the line. the results (in one case the method made the application crash ) induced me to change the approach to the problem, considering the closest vertex of the polyline to the actual intersection.
i know odgefoofooint#d classes aren't implemented yet - i guess they would be in the 2.0 release. may i ask a temporary workaround?
thanks in advance for your suggestions.
mirko pontrelli.
hopefully the following sgi code can help, i had to write my own classes to do this kind of thing as our application is not just a dwg viewer but also a pdf viewer.
i found the following code was a fantastic start to place and understand 3d line intersection, however, the rest is up to moi.
cheers
jason
#include <stdio.h>
#include <math.h>
// demonstrate c function line_line_closest_points3d
// on an sgi, compile as cc -o linelinecp -xcpluscomm linelinecp.c -lm
// then run as ./linelinecp
#define eps (1.0e-7f)
#define minlen (eps)
#define maxreciplen (1.0f / minlen)
// non-homogeneous 3d
typedef struct pointstruct {
float px;
float py;
float pz;
} point;
// plane ax + by + cz + d == 0
typedef struct planestruct {
float a;
float b;
float c;
float d;
} plane;
//
typedef struct vecstruct {
float dx;
float dy;
float dz;
} vector;
j.a.d
float
edist (
point *a,
point *b)
{
float rval;
rval = fsqrt ((b->px - a->px) * (b->px - a->px)
+ (b->py - a->py) * (b->py - a->py)
+ (b->pz - a->pz) * (b->pz - a->pz));
return (rval);
}
j.a.d
// a = p + v
point *
pplusv (
point *a,
point *p,
vector *v)
{
a->px = p->px + v->dx;
a->py = p->py + v->dy;
a->pz = p->pz + v->dz;
return (a);
}
j.a.d
// c = a x b
// careful: c cannot == a or b
vector *
vcross (
vector *c,
vector *a,
vector *b)
{
c->dx = a->dy * b->dz - a->dz * b->dy;
c->dy = a->dz * b->dx - a->dx * b->dz;
c->dz = a->dx * b->dy - a->dy * b->dx;
return (c);
}
j.a.d
float
vdot(
vector *a,
vector *b)
{
return (a->dx * b->dx + a->dy * b->dy + a->dz * b->dz);
}
float
vmag(
vector *a)
{
return (fsqrt(vdot(a, a)));
}
float
recip_vmag(
vector *a)
{
return (1.0 / fsqrt(vdot(a, a)));
}
// a /= ||a||
vector *
vnorm (
vector *a)
{
float d;
if ((d = recip_vmag(a)) > maxreciplen) {
fprintf (stderr, "\nvector at 0x%x: %g %g %g?",
(int)a, a->dx, a->dy, a->dz);
a->dx = a->dy = a->dz = 0.0;
return (null);
}
else {
a->dx *= d;
a->dy *= d;
a->dz *= d;
return (a);
}
}
float
point_plane_dist (
point *a,
plane *p)
{
float rval;
rval = (a->px * p->a + a->py * p->b + a->pz * p->c + p->d);
return (rval);
}
// p = (1-t)a + tb
point *
plerp (
point *p,
point *a,
point *b,
float t)
{
p->px = t * b->px + (1-t) * a->px;
p->py = t * b->py + (1-t) * a->py;
p->pz = t * b->pz + (1-t) * a->pz;
return (p);
}
point *
intersect_line_plane (
point *p, // returned as intersection point
point *a, point *b, // the line to intersect
plane *m) // the plane to intersect
{
float mdota, mdotb, denom, t;
int ip,iq;
mdota = point_plane_dist (a, m);
mdotb = point_plane_dist (b, m);
denom = mdota - mdotb;
if (fabsf(denom / (fabsf(mdota) + fabsf(mdotb))) < eps) {
fprintf (stderr, "int_line_plane(): no intersection?\n");
p->px = p->py = p->pz = 0.0;
return (null);
}
else {
t = mdota / denom;
plerp (p, a, b, t);
return (p);
}
}
// <point,direction> line form
point *
intersect_dline_plane (
point *p, // returned as intersection point
point *a, vector *adir, // the line to intersect
plane *m) // the plane in question
{
static point b, *b = &b;
pplusv (b, a, adir);
return (intersect_line_plane (p, a, b, m));
}
plane *
plane_from_two_vectors_and_point (
plane *m,
vector *u, vector *v, point *p)
{
vnorm (vcross ((vector *)m, u, v));
// plane contains p
m->d = -(m->a * p->px + m->b * p->py + m->c * p->pz);
return (m);
}
// computes pb on line b closest to line a
// computes pa on line a closest to line b
// return: 0 if parallel; 1 if coincident; 2 if generic
int
line_line_closest_points3d (
point *pa, point *pb,
point *a, vector *adir,
point *b, vector *bdir)
{
static vector cdir, *cdir = &cdir;
static plane ac, *ac = &ac, bc, *bc = &bc;
// connecting line is perpendicular to both
vcross (cdir, adir, bdir);
// lines are near-parallel -- all points are closest
if (!vnorm(cdir)) {
*pa = *a;
*pb = *b;
return (0);
}
// form plane containing line a, parallel to cdir
plane_from_two_vectors_and_point (ac, cdir, adir, a);
// form plane containing line b, parallel to cdir
plane_from_two_vectors_and_point (bc, cdir, bdir, b);
// closest point on a is line a ^ bc
intersect_dline_plane (pa, a, adir, bc);
// closest point on b is line b ^ ac
intersect_dline_plane (pb, b, bdir, ac);
if (edist (pa, pb) < eps)
return (1); // coincident
else
return (2); // distinct
}
main ()
{
static point a = { 0.0, 0.0, 0.0 };
static vector da = { 1.0, 1.0, 1.0 };
static point b = { 1.0, 0.0, 0.0 };
static vector db = { 0.0, 0.0, 1.0 };
static point pa, pb;
int rval;
rval = line_line_closest_points3d ( &pa, &pb, &a, &da, &b, &db );
printf ( "l_l_c_p returns %d; pa = %g %g %g ; pb = %g %g %g\n",
rval, pa.px, pa.py, pa.pz, pb.px, pb.py, pb.pz );
}
jason, thanks a lot for your piece of code.
dwgdirect defines the basic data structures in there, so i guess converting that code won't be a hard task i will study it carefully.
greetings.
mirko pontrelli