#ifndef lint static char *sccsid = "@(#)linesegment.c 1.2 (UKC) 9/12/87"; #endif lint /* * Plot the line segment (originally an edge of a polygon), * as far as it is not hidden by the triangles in trilist * Trilist may be half way down the list as this routine is * recursive. */ #include /* for stderr */ #include "types.h" /* * Fudge factors to allow for machine inaccuracy */ #define eps (1.0e-10) /* should be a small no like 1e-5 */ #define meps (-eps) #define oneminus (1-eps) #define oneplus (1+eps) /* * statistics for debugging / optimising * Number of triangles rejected at various tests. */ int ntest1 = 0; int ntest2 = 0; int ntest3 = 0; int ntest4 = 0; int ntest5 = 0; int ntest6 = 0; int ntest6a = 0; int ntest6b = 0; int nplain = 0; /* untouched by all tests */ linesegment(xP, yP, zP, xQ, yQ, zQ, tl) double xP, yP, zP, xQ, yQ, zQ; trilist_t *tl; { trilist_t *tp; int worktodo = 1; /* has this line survived all triangles? */ double a, b, c, h; /* local copies from tp->triangle */ double xA, yA, xB, yB, xC, yC; /* ditto */ double hP, hQ; double dA, dB, dC; /* for test 2 */ double xN, yN; /* normal to PQ */ int Poutside, Qoutside; /* bools for test 3 */ double xPQ, yPQ, zPQ; /* for test 3 */ double xAB, yAB; double temp; /* for constant calcs outside loops and */ /* swapping and shuffling of variables. */ double MIN, MAX; /* lambda values for I and J */ double eps1; double mu, lab; int minset, maxset; int i; double xmin, ymin, zmin; double xmax, ymax, zmax; for (tp = tl; tp != NULLTRI; tp = tp->link) { a = tp->a; b = tp->b; c = tp->c; h = tp->h; eps1 = eps + eps*h; /* test 1: * If PQ is in front of or in ABC, ABC does not obscure it. * includes PQ is an edge of ABC. */ hP = a*xP + b*yP + c*zP; hQ = a*xQ + b*yQ + c*zQ; if (hP <= h+eps1 && hQ <= h+eps1) { ntest1++; continue; } /* need these now */ xA = tp->A.x; yA = tp->A.y; xB = tp->B.x; yB = tp->B.y; xC = tp->C.x; yC = tp->C.y; /* test 2: * If A, B and C lie on the same side of the * infinite line through PQ, PQ is visually outside the * triangle, hence not obscured by it. * * Lying on it counts as not being on the same side, * otherwise a line exactly behind a crack shows through. * * For a flat projection, this is a 2D test for [ABC][xy] * being on the same side of line [PQ][xy] */ /* compute a normal to PQ */ xN = yQ - yP; yN = xP - xQ; /* sign of dot product tells us which side of PQ it is on. */ dA = xN * (xA - xP) + yN * (yA - yP); dB = xN * (xB - xP) + yN * (yB - yP); dC = xN * (xC - xP) + yN * (yC - yP); /* if all points are on the same side of the line, * then the triangle does not hide the line. */ if ((dA > eps && dB > eps && dC > eps) || (dA < meps && dB < meps && dC < meps )) { ntest2++; continue; } /* test 3: * infinite line through PQ definitely goes though ABC. * * If P and Q are both outside the triangle on the same side, * the triangle does not obscure the line. * * Again, 2D calculations because of flat projection instead * of perspective. */ Poutside = Qoutside = 0; /* process the three sides of the triangle */ for (i = 0; i < 3; i++) { int Pbeyond, Qbeyond; /* are P,Q outside this line? */ /* compute a normal to AB, pointing outwards. */ xN = yB - yA; yN = xA - xB; /* If both P and Q are outside AB, * this triangle does not obscure the line. * If AB and PQ are coincident, we still need to clip * the line otherwise lines show through cracks. */ Pbeyond = (xN * (xP - xA) + yN * (yP - yA) > eps); Qbeyond = (xN * (xQ - xA) + yN * (yQ - yA) > eps); if (Pbeyond && Qbeyond) { ntest3++; goto nexttriangle; } if (Pbeyond) Poutside = 1; if (Qbeyond) Qoutside = 1; /* shuffle ABC round (they end up in the same place) */ temp = xA; xA = xB; xB = xC; xC = temp; temp = yA; yA = yB; yB = yC; yC = temp; } /* test 4 * If both P & Q are inside triangle, the line is hidden. */ if (! (Poutside || Qoutside)) { worktodo = 0; ntest4++; break; } /* test 5: * Find I and J and if they are before the plane, no prob. */ xPQ = xQ - xP; yPQ = yQ - yP; MIN = 2.0; minset = 0; MAX = -1.0; maxset = 0; for (i=0; i<3; i++) { /* compute intersection of PQ and AB */ xAB = xB - xA; yAB = yB - yA; temp = xAB * yPQ - yAB * xPQ ; if (temp == 0) { /* PQ and AB are parallel */ } else { mu = ((xP-xA) * yPQ - (yP-yA) * xPQ) / temp; lab= ((xA-xP) * yAB - (yA-yP) * xAB) / -temp; /* stash I and J for later tests */ if ( mu >= meps && mu <= oneplus && lab >= meps && lab <= oneplus ) { if (lab < MIN) { MIN = lab; minset = 1; } if (lab > MAX) { MAX = lab; maxset = 1; } } } /* shuffle ABC round */ temp = xA; xA = xB; xB = xC; xC = temp; temp = yA; yA = yB; yB = yC; yC = temp; } /* just check... */ if (!maxset) { if (!minset) error("MAX nor MIN set"); else error("MAX not set"); } else if (!minset) error("MIN not set"); else goto noboob; /* verily, what wondrous flow control constructs! Ugh! */ goto impossible; noboob: zPQ = zQ - zP; xmin = xP + MIN * xPQ; ymin = yP + MIN * yPQ; zmin = zP + MIN * zPQ; if (a*xmin + b*ymin + c*zmin < h-eps1) { ntest5++; continue; } xmax = xP + MAX * xPQ; ymax = yP + MAX * yPQ; zmax = zP + MAX * zPQ; if (a*xmax + b*ymax + c*zmax < h-eps1) { error("Hum! Have you been giving me intersecting polygons?"); impossible: fprintf(stderr, "while trying to clip (%g %g %g) (%g %g %g)\nagainst (%g %g %g) (%g %g %g) (%g %g %g)\n", xP, yP, zP, xQ, yQ, zQ, xA, yA, tp->A.z, xB, yB, tp->B.z, xC, yC, tp->C.z); continue; } /* test 6 * don't really understand this. */ if (Poutside || hP < h-eps1) { linesegment(xP, yP, zP, xmin, ymin, zmin, tp->link); ntest6a++; } if (Qoutside || hQ < h-eps1) { linesegment(xmax, ymax, zmax, xQ, yQ, zQ, tp->link); ntest6b++; } ntest6++; worktodo = 0; break; nexttriangle: ; } if (worktodo) { nplain++; outputline(xP, yP, xQ, yQ); } } printstats() { fprintf(stderr, "1:%d 2:%d 3:%d 4:%d 5:%d 6:%d 6a:%d 6b:%d p:%d\n", ntest1, ntest2, ntest3, ntest4, ntest5, ntest6, ntest6a, ntest6b, nplain); }