/*------------------------------------------------------------------------- Create a contour slice through a 3 vertex facet pa,pb,pc The contour "level" is a horizontal plane perpendicular to the z axis, ie: The equation of the contour plane Ax + By + Cz + D = 0 has A = 0, B = 0, C = 1, D = -level Return 0 if the contour plane doesn't cut the facet 2 if it does cut the facet -1 for an unexpected occurrence If a vertex touches the contour plane nothing need to be drawn!? */ int ContourFace(XYZ pa,XYZ pb,XYZ pc,double level,XYZ *p1,XYZ *p2) { double sidea,sideb,sidec; /* Evaluate the equation of the plane for each vertex sidea = A * pa.x + B * pa.y + C * pa.z + D; sideb = A * pb.x + B * pb.y + C * pb.z + D; sidec = A * pc.x + B * pc.y + C * pc.z + D; */ sidea = pa.z - level; sideb = pb.z - level; sidec = pc.z - level; /* Are all the vertices on one side */ if (sidea >= 0 && sideb >= 0 && sidec >= 0) return(0); if (sidea <= 0 && sideb <= 0 && sidec <= 0) return(0); /* Is p0 the only point on a side by itself */ if ((SIGN(sidea) != SIGN(sideb)) && (SIGN(sidea) != SIGN(sidec))) { p1->x = pa.x - sidea * (pc.x - pa.x) / (sidec - sidea); p1->y = pa.y - sidea * (pc.y - pa.y) / (sidec - sidea); p1->z = pa.z - sidea * (pc.z - pa.z) / (sidec - sidea); p2->x = pa.x - sidea * (pb.x - pa.x) / (sideb - sidea); p2->y = pa.y - sidea * (pb.y - pa.y) / (sideb - sidea); p2->z = pa.z - sidea * (pb.z - pa.z) / (sideb - sidea); return(2); } /* Is p1 the only point on a side by itself */ if ((SIGN(sideb) != SIGN(sidea)) && (SIGN(sideb) != SIGN(sidec))) { p1->x = pb.x - sideb * (pc.x - pb.x) / (sidec - sideb); p1->y = pb.y - sideb * (pc.y - pb.y) / (sidec - sideb); p1->z = pb.z - sideb * (pc.z - pb.z) / (sidec - sideb); p2->x = pb.x - sideb * (pa.x - pb.x) / (sidea - sideb); p2->y = pb.y - sideb * (pa.y - pb.y) / (sidea - sideb); p2->z = pb.z - sideb * (pa.z - pb.z) / (sidea - sideb); return(2); } /* Is p2 the only point on a side by itself */ if ((SIGN(sidec) != SIGN(sidea)) && (SIGN(sidec) != SIGN(sideb))) { p1->x = pc.x - sidec * (pa.x - pc.x) / (sidea - sidec); p1->y = pc.y - sidec * (pa.y - pc.y) / (sidea - sidec); p1->z = pc.z - sidec * (pa.z - pc.z) / (sidea - sidec); p2->x = pc.x - sidec * (pb.x - pc.x) / (sideb - sidec); p2->y = pc.y - sidec * (pb.y - pc.y) / (sideb - sidec); p2->z = pc.z - sidec * (pb.z - pc.z) / (sideb - sidec); return(2); } /* Shouldn't get here */ return(-1); }