//-----------------------------------------------------------------------------
// How to intersect two surfaces, to get some type of curve. This is either
// an exact special case (e.g., two planes to make a line), or a numerical
// thing.
//
// Copyright 2008-2013 Jonathan Westhues.
//-----------------------------------------------------------------------------
#include "solvespace.h"
extern int FLAG;
void SSurface::AddExactIntersectionCurve(SBezier *sb, SSurface *srfB,
SShell *agnstA, SShell *agnstB, SShell *into)
{
SCurve sc = {};
// Important to keep the order of (surfA, surfB) consistent; when we later
// rewrite the identifiers, we rewrite surfA from A and surfB from B.
sc.surfA = h;
sc.surfB = srfB->h;
sc.exact = *sb;
sc.isExact = true;
// Now we have to piecewise linearize the curve. If there's already an
// identical curve in the shell, then follow that pwl exactly, otherwise
// calculate from scratch.
SCurve split, *existing = NULL;
SBezier sbrev = *sb;
sbrev.Reverse();
bool backwards = false;
#pragma omp critical(into)
{
for(SCurve &se : into->curve) {
if(se.isExact) {
if(sb->Equals(&(se.exact))) {
existing = &se;
break;
}
if(sbrev.Equals(&(se.exact))) {
existing = &se;
backwards = true;
break;
}
}
}
}// end omp critical
if(existing) {
SCurvePt *v;
for(v = existing->pts.First(); v; v = existing->pts.NextAfter(v)) {
sc.pts.Add(v);
}
if(backwards) sc.pts.Reverse();
split = sc;
sc = {};
} else {
sb->MakePwlInto(&(sc.pts));
// and split the line where it intersects our existing surfaces
split = sc.MakeCopySplitAgainst(agnstA, agnstB, this, srfB);
sc.Clear();
}
// Test if the curve lies entirely outside one of the
SCurvePt *scpt;
bool withinA = false, withinB = false;
for(scpt = split.pts.First(); scpt; scpt = split.pts.NextAfter(scpt)) {
double tol = 0.01;
Point2d puv;
ClosestPointTo(scpt->p, &puv);
if(puv.x > -tol && puv.x < 1 + tol &&
puv.y > -tol && puv.y < 1 + tol)
{
withinA = true;
}
srfB->ClosestPointTo(scpt->p, &puv);
if(puv.x > -tol && puv.x < 1 + tol &&
puv.y > -tol && puv.y < 1 + tol)
{
withinB = true;
}
// Break out early, no sense wasting time if we already have the answer.
if(withinA && withinB) break;
}
if(!(withinA && withinB)) {
// Intersection curve lies entirely outside one of the surfaces, so
// it's fake.
split.Clear();
return;
}
#if 0
if(sb->deg == 2) {
dbp(" ");
SCurvePt *prev = NULL, *v;
dbp("split.pts.n = %d", split.pts.n);
for(v = split.pts.First(); v; v = split.pts.NextAfter(v)) {
if(prev) {
Vector e = (prev->p).Minus(v->p).WithMagnitude(0);
SS.nakedEdges.AddEdge((prev->p).Plus(e), (v->p).Minus(e));
}
prev = v;
}
}
#endif // 0
ssassert(!(sb->Start()).Equals(sb->Finish()),
"Unexpected zero-length edge");
split.source = SCurve::Source::INTERSECTION;
#pragma omp critical(into)
{
into->curve.AddAndAssignId(&split);
}
}
void SSurface::IntersectAgainst(SSurface *b, SShell *agnstA, SShell *agnstB,
SShell *into)
{
Vector amax, amin, bmax, bmin;
GetAxisAlignedBounding(&amax, &amin);
b->GetAxisAlignedBounding(&bmax, &bmin);
if(Vector::BoundingBoxesDisjoint(amax, amin, bmax, bmin)) {
// They cannot possibly intersect, no curves to generate
return;
}
Vector alongt, alongb;
SBezier oft, ofb;
bool isExtdt = this->IsExtrusion(&oft, &alongt),
isExtdb = b->IsExtrusion(&ofb, &alongb);
if(degm == 1 && degn == 1 && b->degm == 1 && b->degn == 1) {
// Line-line intersection; it's a plane or nothing.
Vector na = NormalAt(0, 0).WithMagnitude(1),
nb = b->NormalAt(0, 0).WithMagnitude(1);
double da = na.Dot(PointAt(0, 0)),
db = nb.Dot(b->PointAt(0, 0));
Vector dl = na.Cross(nb);
if(dl.Magnitude() < LENGTH_EPS) return; // parallel planes
dl = dl.WithMagnitude(1);
Vector p = Vector::AtIntersectionOfPlanes(na, da, nb, db);
// Trim it to the region 0 <= {u,v} <= 1 for each plane; not strictly
// necessary, since line will be split and excess edges culled, but
// this improves speed and robustness.
int i;
double tmax = VERY_POSITIVE, tmin = VERY_NEGATIVE;
for(i = 0; i < 2; i++) {
SSurface *s = (i == 0) ? this : b;
Vector tu, tv;
s->TangentsAt(0, 0, &tu, &tv);
double up, vp, ud, vd;
s->ClosestPointTo(p, &up, &vp);
ud = (dl.Dot(tu)) / tu.MagSquared();
vd = (dl.Dot(tv)) / tv.MagSquared();
// so u = up + t*ud
// v = vp + t*vd
if(ud > LENGTH_EPS) {
tmin = max(tmin, -up/ud);
tmax = min(tmax, (1 - up)/ud);
} else if(ud < -LENGTH_EPS) {
tmax = min(tmax, -up/ud);
tmin = max(tmin, (1 - up)/ud);
} else {
if(up < -LENGTH_EPS || up > 1 + LENGTH_EPS) {
// u is constant, and outside [0, 1]
tmax = VERY_NEGATIVE;
}
}
if(vd > LENGTH_EPS) {
tmin = max(tmin, -vp/vd);
tmax = min(tmax, (1 - vp)/vd);
} else if(vd < -LENGTH_EPS) {
tmax = min(tmax, -vp/vd);
tmin = max(tmin, (1 - vp)/vd);
} else {
if(vp < -LENGTH_EPS || vp > 1 + LENGTH_EPS) {
// v is constant, and outside [0, 1]
tmax = VERY_NEGATIVE;
}
}
}
if(tmax > tmin + LENGTH_EPS) {
SBezier bezier = SBezier::From(p.Plus(dl.ScaledBy(tmin)),
p.Plus(dl.ScaledBy(tmax)));
AddExactIntersectionCurve(&bezier, b, agnstA, agnstB, into);
}
} else if((degm == 1 && degn == 1 && isExtdb) ||
(b->degm == 1 && b->degn == 1 && isExtdt))
{
// The intersection between a plane and a surface of extrusion
SSurface *splane, *sext;
if(degm == 1 && degn == 1) {
splane = this;
sext = b;
} else {
splane = b;
sext = this;
}
Vector n = splane->NormalAt(0, 0).WithMagnitude(1), along;
double d = n.Dot(splane->PointAt(0, 0));
SBezier bezier;
(void)sext->IsExtrusion(&bezier, &along);
if(fabs(n.Dot(along)) < LENGTH_EPS) {
// Direction of extrusion is parallel to plane; so intersection
// is zero or more lines. Build a line within the plane, and
// normal to the direction of extrusion, and intersect that line
// against the surface; each intersection point corresponds to
// a line.
Vector pm, alu, p0, dp;
// a point halfway along the extrusion
pm = ((sext->ctrl[0][0]).Plus(sext->ctrl[0][1])).ScaledBy(0.5);
alu = along.WithMagnitude(1);
dp = (n.Cross(along)).WithMagnitude(1);
// n, alu, and dp form an orthogonal csys; set n component to
// place it on the plane, alu component to lie halfway along
// extrusion, and dp component doesn't matter so zero
p0 = n.ScaledBy(d).Plus(alu.ScaledBy(pm.Dot(alu)));
List