SBezier ret = {};

ret.deg = 1;

ret.weight[0] = p0.w;

ret.ctrl [0] = p0.PerspectiveProject();

ret.weight[1] = p1.w;

ret.ctrl [1] = p1.PerspectiveProject();

return ret;

SBezier ret = {};

ret.deg = 2;

ret.weight[0] = p0.w;

ret.ctrl [0] = p0.PerspectiveProject();

ret.weight[1] = p1.w;

ret.ctrl [1] = p1.PerspectiveProject();

ret.weight[2] = p2.w;

ret.ctrl [2] = p2.PerspectiveProject();

return ret;

SBezier ret = {};

ret.deg = 3;

ret.weight[0] = p0.w;

ret.ctrl [0] = p0.PerspectiveProject();

ret.weight[1] = p1.w;

ret.ctrl [1] = p1.PerspectiveProject();

ret.weight[2] = p2.w;

ret.ctrl [2] = p2.PerspectiveProject();

ret.weight[3] = p3.w;

ret.ctrl [3] = p3.PerspectiveProject();

return ret;

return SBezier::From(p0.Project4d(),

p1.Project4d(),

p2.Project4d());

return SBezier::From(p0.Project4d(),

p1.Project4d(),

p2.Project4d(),

p3.Project4d());

int i;

for(i = 0; i < (deg+1)/2; i++) {

swap(ctrl[i], ctrl[deg-i]);

swap(weight[i], weight[deg-i]);

}

int i;

for(i = 0; i <= deg; i++) {

ctrl[i] = ctrl[i].ScaledBy(s);

}

double *umin, double *umax) const

{

int i;

for(i = 0; i <= deg; i++) {

double ut = ((ctrl[i]).Minus(orig)).Dot(u);

if(ut < *umin) *umin = ut;

if(ut > *umax) *umax = ut;

}

SBezier ret = *this;

int i;

for(i = 0; i <= deg; i++) {

ret.ctrl[i] = (ret.ctrl[i]).ScaledBy(scale);

ret.ctrl[i] = (q.Rotate(ret.ctrl[i])).Plus(t);

}

return ret;

int i;

for(i = 0; i <= deg; i++) {

if(fabs((ctrl[i]).Dot(n) - d) > LENGTH_EPS) {

return false;

}

}

return true;

if(deg != 2) return false;

if(ctrl[1].DistanceToLine(ctrl[0], ctrl[2].Minus(ctrl[0])) < LENGTH_EPS) {

// This is almost a line segment. So it's a circle with very large

// radius, which is likely to make code that tries to handle circles

// blow up. So return false.

return false;

}

Vector t0 = (ctrl[0]).Minus(ctrl[1]),

t2 = (ctrl[2]).Minus(ctrl[1]),

r0 = axis.Cross(t0),

r2 = axis.Cross(t2);

*center = Vector::AtIntersectionOfLines(ctrl[0], (ctrl[0]).Plus(r0),

ctrl[2], (ctrl[2]).Plus(r2),

NULL, NULL, NULL);

double rd0 = center->Minus(ctrl[0]).Magnitude(),

rd2 = center->Minus(ctrl[2]).Magnitude();

if(fabs(rd0 - rd2) > LENGTH_EPS) {

return false;

}

*r = rd0;

Vector u = r0.WithMagnitude(1),

v = (axis.Cross(u)).WithMagnitude(1);

Point2d c2 = center->Project2d(u, v),

pa2 = (ctrl[0]).Project2d(u, v).Minus(c2),

pb2 = (ctrl[2]).Project2d(u, v).Minus(c2);

double thetaa = atan2(pa2.y, pa2.x), // in fact always zero due to csys

thetab = atan2(pb2.y, pb2.x),

dtheta = WRAP_NOT_0(thetab - thetaa, 2*PI);

if(dtheta > PI) {

// Not possible with a second order Bezier arc; so we must have

// the points backwards.

dtheta = 2*PI - dtheta;

}

if(fabs(weight[1] - cos(dtheta/2)) > LENGTH_EPS) {

return false;

}

return true;

int i;

for(i = 0; i <= deg; i++) {

if(fabs(weight[i] - 1) > LENGTH_EPS) return true;

}

return false;

Vector origin, double cameraTan) const

{

Quaternion q = Quaternion::From(u, v);

q = q.Inverse();

// we want Q*(p - o) = Q*p - Q*o

SBezier ret = this->TransformedBy(q.Rotate(origin).ScaledBy(-1), q, 1.0);

int i;

for(i = 0; i <= deg; i++) {

Vector4 ct = Vector4::From(ret.weight[i], ret.ctrl[i]);

// so the desired curve, before perspective, is

// (x/w, y/w, z/w)

// and after perspective is

// ((x/w)/(1 - (z/w)*cameraTan, ...

// = (x/(w - z*cameraTan), ...

// so we want to let w' = w - z*cameraTan

ct.w = ct.w - ct.z*cameraTan;

ret.ctrl[i] = ct.PerspectiveProject();

ret.weight[i] = ct.w;

}

return ret;

// We just test of identical degree and control points, even though two

// curves could still be coincident (even sharing endpoints).

if(deg != b->deg) return false;

int i;

for(i = 0; i <= deg; i++) {

if(!(ctrl[i]).Equals(b->ctrl[i])) return false;

if(fabs(weight[i] - b->weight[i]) > LENGTH_EPS) return false;

}

return true;

SBezier *sb;

for(sb = l.First(); sb; sb = l.NextAfter(sb)) {

sb->ScaleSelfBy(s);

}

int i, j;

l.ClearTags();

for(i = 0; i < l.n; i++) {

SBezier *bi = &(l.elem[i]), bir;

bir = *bi;

bir.Reverse();

for(j = i + 1; j < l.n; j++) {

SBezier *bj = &(l.elem[j]);

if(bj->Equals(bi) ||

bj->Equals(&bir))

{

bi->tag = 1;

bj->tag = 1;

}

}

}

l.RemoveTagged();

for(const SBezier *sba = l.First(); sba; sba = l.NextAfter(sba)) {

for(const SBezier *sbb = sblb->l.First(); sbb; sbb = sblb->l.NextAfter(sbb)) {

sbb->AllIntersectionsWith(sba, spl);

}

}

SPointList splRaw = {};

SEdgeList sea, seb;

sea = {};

seb = {};

this->MakePwlInto(&sea);

sbb ->MakePwlInto(&seb);

SEdge *se;

for(se = sea.l.First(); se; se = sea.l.NextAfter(se)) {

// This isn't quite correct, since AnyEdgeCrossings doesn't count

// the case where two pairs of line segments intersect at their

// vertices. So this isn't robust, although that case isn't very

// likely.

seb.AnyEdgeCrossings(se->a, se->b, NULL, &splRaw);

}

SPoint *sp;

for(sp = splRaw.l.First(); sp; sp = splRaw.l.NextAfter(sp)) {

Vector p = sp->p;

if(PointOnThisAndCurve(sbb, &p)) {

if(!spl->ContainsPoint(p)) spl->Add(p);

}

}

sea.Clear();

seb.Clear();

splRaw.Clear();

Vector *notCoplanarAt) const

{

Vector pt, ptFar, ptOffLine, dp, n;

double farMax, offLineMax;

int i;

// Get any point on any Bezier; or an arbitrary point if list is empty.

if(l.n > 0) {

pt = l.elem[0].Start();

} else {

pt = Vector::From(0, 0, 0);

}

ptFar = ptOffLine = pt;

// Get the point farthest from our arbitrary point.

farMax = VERY_NEGATIVE;

for(const SBezier *sb = l.First(); sb; sb = l.NextAfter(sb)) {

for(i = 0; i <= sb->deg; i++) {

double m = (pt.Minus(sb->ctrl[i])).Magnitude();

if(m > farMax) {

ptFar = sb->ctrl[i];

farMax = m;

}

}

}

if(ptFar.Equals(pt)) {

// The points are all coincident. So neither basis vector matters.

*p = pt;

*u = Vector::From(1, 0, 0);

*v = Vector::From(0, 1, 0);

return true;

}

// Get the point farthest from the line between pt and ptFar

dp = ptFar.Minus(pt);

offLineMax = VERY_NEGATIVE;

for(const SBezier *sb = l.First(); sb; sb = l.NextAfter(sb)) {

for(i = 0; i <= sb->deg; i++) {

double m = (sb->ctrl[i]).DistanceToLine(pt, dp);

if(m > offLineMax) {

ptOffLine = sb->ctrl[i];

offLineMax = m;

}

}

}

*p = pt;

if(offLineMax < LENGTH_EPS) {

// The points are all colinear; so choose the second basis vector

// arbitrarily.

*u = (ptFar.Minus(pt)).WithMagnitude(1);

*v = (u->Normal(0)).WithMagnitude(1);

} else {

// The points actually define a plane.

n = (ptFar.Minus(pt)).Cross(ptOffLine.Minus(pt));

*u = (n.Normal(0)).WithMagnitude(1);

*v = (n.Normal(1)).WithMagnitude(1);

}

// So we have a plane; but check whether all of the points lie in that

// plane.

n = u->Cross(*v);

n = n.WithMagnitude(1);

double d = p->Dot(n);

for(const SBezier *sb = l.First(); sb; sb = l.NextAfter(sb)) {

for(i = 0; i <= sb->deg; i++) {

if(fabs(n.Dot(sb->ctrl[i]) - d) > LENGTH_EPS) {

if(notCoplanarAt) *notCoplanarAt = sb->ctrl[i];

return false;

}

}

}

return true;

bool *allClosed, SEdge *errorAt)

{

SBezierLoop loop = {};

if(sbl->l.n < 1) return loop;

sbl->l.ClearTags();

SBezier *first = &(sbl->l.elem[0]);

first->tag = 1;

loop.l.Add(first);

Vector start = first->Start();

Vector hanging = first->Finish();

int auxA = first->auxA;

sbl->l.RemoveTagged();

while(sbl->l.n > 0 && !hanging.Equals(start)) {

int i;

bool foundNext = false;

for(i = 0; i < sbl->l.n; i++) {

SBezier *test = &(sbl->l.elem[i]);

if((test->Finish()).Equals(hanging) && test->auxA == auxA) {

test->Reverse();

// and let the next test catch it

}

if((test->Start()).Equals(hanging) && test->auxA == auxA) {

test->tag = 1;

loop.l.Add(test);

hanging = test->Finish();

sbl->l.RemoveTagged();

foundNext = true;

break;

}

}

if(!foundNext) {

// The loop completed without finding the hanging edge, so

// it's an open loop

errorAt->a = hanging;

errorAt->b = start;

*allClosed = false;

return loop;

}

}

if(hanging.Equals(start)) {

*allClosed = true;

} else {

// We ran out of edges without forming a closed loop.

errorAt->a = hanging;

errorAt->b = start;

*allClosed = false;

}

return loop;

l.Reverse();

SBezier *sb;

for(sb = l.First(); sb; sb = l.NextAfter(sb)) {

// If we didn't reverse each curve, then the next curve in list would

// share your start, not your finish.

sb->Reverse();

}

double *umin, double *umax) const

{

for(const SBezier *sb = l.First(); sb; sb = l.NextAfter(sb)) {

sb->GetBoundingProjd(u, orig, umin, umax);

}

for(const SBezier *sb = l.First(); sb; sb = l.NextAfter(sb)) {

sb->MakePwlInto(sc, chordTol);

// Avoid double points at join between Beziers; except that

// first and last points should be identical.

if(l.NextAfter(sb) != NULL) {

sc->l.RemoveLast(1);

}

}

// Ensure that it's exactly closed, not just within a numerical tolerance.

if((sc->l.elem[sc->l.n - 1].p).Equals(sc->l.elem[0].p)) {

sc->l.elem[sc->l.n - 1] = sc->l.elem[0];

}

if(l.n < 1) return false;

Vector s = l.elem[0].Start(),

f = l.elem[l.n-1].Finish();

return s.Equals(f);

double chordTol,

bool *allClosed, SEdge *errorAt,

SBezierList *openContours)

{

SBezierLoopSet ret = {};

*allClosed = true;

while(sbl->l.n > 0) {

bool thisClosed;

SBezierLoop loop;

loop = SBezierLoop::FromCurves(sbl, &thisClosed, errorAt);

if(!thisClosed) {

// Record open loops in a separate list, if requested.

*allClosed = false;

if(openContours) {

SBezier *sb;

for(sb = loop.l.First(); sb; sb = loop.l.NextAfter(sb)) {

openContours->l.Add(sb);

}

}

loop.Clear();

} else {

ret.l.Add(&loop);

poly->AddEmptyContour();

loop.MakePwlInto(&(poly->l.elem[poly->l.n-1]), chordTol);

}

}

poly->normal = poly->ComputeNormal();

ret.normal = poly->normal;

if(poly->l.n > 0) {

ret.point = poly->AnyPoint();

} else {

ret.point = Vector::From(0, 0, 0);

}

return ret;

double *umin, double *umax) const

{

for(const SBezierLoop *sbl = l.First(); sbl; sbl = l.NextAfter(sbl)) {

sbl->GetBoundingProjd(u, orig, umin, umax);

}

for(const SBezierLoop *sbl = l.First(); sbl; sbl = l.NextAfter(sbl)) {

sp->AddEmptyContour();

sbl->MakePwlInto(&(sp->l.elem[sp->l.n - 1]));

}

int i;

for(i = 0; i < l.n; i++) {

(l.elem[i]).Clear();

}

l.Clear();

SSurface *srfuv,

double chordTol,

bool *allClosed, SEdge *notClosedAt,

bool *allCoplanar, Vector *notCoplanarAt,

SBezierList *openContours)

{

SSurface srfPlane;

if(!srfuv) {

Vector p, u, v;

*allCoplanar =

sbl->GetPlaneContainingBeziers(&p, &u, &v, notCoplanarAt);

if(!*allCoplanar) {

// Don't even try to assemble them into loops if they're not

// all coplanar.

if(openContours) {

SBezier *sb;

for(sb = sbl->l.First(); sb; sb = sbl->l.NextAfter(sb)) {

openContours->l.Add(sb);

}

}

return;

}

// All the curves lie in a plane through p with basis vectors u and v.

srfPlane = SSurface::FromPlane(p, u, v);

srfuv = &srfPlane;

}

int i, j;

// Assemble the Bezier trim curves into closed loops; we also get the

// piecewise linearization of the curves (in the SPolygon spxyz), as a

// calculation aid for the loop direction.

SBezierLoopSet sbls = SBezierLoopSet::From(sbl, spxyz, chordTol,

allClosed, notClosedAt,

openContours);

if(sbls.l.n != spxyz->l.n) return;

// Convert the xyz piecewise linear to uv piecewise linear.

SPolygon spuv = {};

SContour *sc;

for(sc = spxyz->l.First(); sc; sc = spxyz->l.NextAfter(sc)) {

spuv.AddEmptyContour();

SPoint *pt;

for(pt = sc->l.First(); pt; pt = sc->l.NextAfter(pt)) {

double u, v;

srfuv->ClosestPointTo(pt->p, &u, &v);

spuv.l.elem[spuv.l.n - 1].AddPoint(Vector::From(u, v, 0));

}

}

spuv.normal = Vector::From(0, 0, 1); // must be, since it's in xy plane now

static const int OUTER_LOOP = 10;

static const int INNER_LOOP = 20;

static const int USED_LOOP = 30;

// Fix the contour directions; we do this properly, in uv space, so it

// works for curved surfaces too (important for STEP export).

spuv.FixContourDirections();

for(i = 0; i < spuv.l.n; i++) {

SContour *contour = &(spuv.l.elem[i]);

SBezierLoop *bl = &(sbls.l.elem[i]);

if(contour->tag) {

// This contour got reversed in the polygon to make the directions

// consistent, so the same must be necessary for the Bezier loop.

bl->Reverse();

}

if(contour->IsClockwiseProjdToNormal(spuv.normal)) {

bl->tag = INNER_LOOP;

} else {

bl->tag = OUTER_LOOP;

}

}

bool loopsRemaining = true;

while(loopsRemaining) {

loopsRemaining = false;

for(i = 0; i < sbls.l.n; i++) {

SBezierLoop *loop = &(sbls.l.elem[i]);

if(loop->tag != OUTER_LOOP) continue;

// Check if this contour contains any outer loops; if it does, then

// we should do those "inner outer loops" first; otherwise we

// will steal their holes, since their holes also lie inside this

// contour.

for(j = 0; j < sbls.l.n; j++) {

SBezierLoop *outer = &(sbls.l.elem[j]);

if(i == j) continue;

if(outer->tag != OUTER_LOOP) continue;

Vector p = spuv.l.elem[j].AnyEdgeMidpoint();

if(spuv.l.elem[i].ContainsPointProjdToNormal(spuv.normal, p)) {

break;

}

}

if(j < sbls.l.n) {

// It does, can't do this one yet.

continue;

}

SBezierLoopSet outerAndInners = {};

loopsRemaining = true;

loop->tag = USED_LOOP;

outerAndInners.l.Add(loop);

int auxA = 0;

if(loop->l.n > 0) auxA = loop->l.elem[0].auxA;

for(j = 0; j < sbls.l.n; j++) {

SBezierLoop *inner = &(sbls.l.elem[j]);

if(inner->tag != INNER_LOOP) continue;

if(inner->l.n < 1) continue;

if(inner->l.elem[0].auxA != auxA) continue;

Vector p = spuv.l.elem[j].AnyEdgeMidpoint();

if(spuv.l.elem[i].ContainsPointProjdToNormal(spuv.normal, p)) {

outerAndInners.l.Add(inner);

inner->tag = USED_LOOP;

}

}

outerAndInners.point = srfuv->PointAt(0, 0);

outerAndInners.normal = srfuv->NormalAt(0, 0);

l.Add(&outerAndInners);

}

}

// If we have poorly-formed loops--for example, overlapping zero-area

// stuff--then we can end up with leftovers. We use this function to

// group stuff into closed paths for export when possible, so it's bad

// to screw up on that stuff. So just add them onto the open curve list.

// Very ugly, but better than losing curves.

for(i = 0; i < sbls.l.n; i++) {

SBezierLoop *loop = &(sbls.l.elem[i]);

if(loop->tag == USED_LOOP) continue;

if(openContours) {

SBezier *sb;

for(sb = loop->l.First(); sb; sb = loop->l.NextAfter(sb)) {

openContours->l.Add(sb);

}

}

loop->Clear();

// but don't free the used loops, since we shallow-copied them to

// ourself

}

sbls.l.Clear(); // not sbls.Clear(), since that would deep-clear

spuv.Clear();

SBezierLoop sbl = {};

sbl.l.Add(sb);

SBezierLoopSet sbls = {};

sbls.l.Add(&sbl);

l.Add(&sbls);

SBezierLoopSet *sbls;

for(sbls = l.First(); sbls; sbls = l.NextAfter(sbls)) {

sbls->Clear();

}

l.Clear();

Quaternion q, double scale)

{

bool needRotate = !EXACT(q.vx == 0.0 && q.vy == 0.0 && q.vz == 0.0 && q.w == 1.0);

bool needTranslate = !EXACT(t.x == 0.0 && t.y == 0.0 && t.z == 0.0);

bool needScale = !EXACT(scale == 1.0);

SCurve ret = {};

ret.h = a->h;

ret.isExact = a->isExact;

ret.exact = (a->exact).TransformedBy(t, q, scale);

ret.surfA = a->surfA;

ret.surfB = a->surfB;

SCurvePt *p;

for(p = a->pts.First(); p; p = a->pts.NextAfter(p)) {

SCurvePt pp = *p;

if(needScale) {

pp.p = (pp.p).ScaledBy(scale);

}

if(needRotate) {

pp.p = q.Rotate(pp.p);

}

if(needTranslate) {

pp.p = pp.p.Plus(t);

}

ret.pts.Add(&pp);

}

return ret;

// Three, not two; curves are pwl'd to at least two edges (three points)

// even if not necessary, to avoid square holes.

if(pts.n <= 3) return;

pts.ClearTags();

Vector prev = pts.elem[0].p;

int i, a;

for(i = 1; i < pts.n - 1; i++) {

SCurvePt *sct = &(pts.elem[i]),

*scn = &(pts.elem[i+1]);

if(sct->vertex) {

prev = sct->p;

continue;

}

bool mustKeep = false;

// We must check against both surfaces; the piecewise linear edge

// may have a different chord tolerance in the two surfaces. (For

// example, a circle in the surface of a cylinder is just a straight

// line, so it always has perfect chord tol, but that circle in

// a plane is a circle so it doesn't).

for(a = 0; a < 2; a++) {

SSurface *srf = (a == 0) ? srfA : srfB;

Vector puv, nuv;

srf->ClosestPointTo(prev, &(puv.x), &(puv.y));

srf->ClosestPointTo(scn->p, &(nuv.x), &(nuv.y));

if(srf->ChordToleranceForEdge(nuv, puv) > SS.ChordTolMm()) {

mustKeep = true;

}

}

if(mustKeep) {

prev = sct->p;

} else {

sct->tag = 1;

// and prev is unchanged, since there's no longer any point

// in between

}

}

pts.RemoveTagged();

STrimBy stb = {};

stb.curve = hsc;

SCurve *sc = shell->curve.FindById(hsc);

if(backwards) {

stb.finish = sc->pts.elem[0].p;

stb.start = sc->pts.elem[sc->pts.n - 1].p;

stb.backwards = true;

} else {

stb.start = sc->pts.elem[0].p;

stb.finish = sc->pts.elem[sc->pts.n - 1].p;

stb.backwards = false;

}

return stb;