C++ SSurface类代码示例

//-----------------------------------------------------------------------------
// In our shell, find all surfaces that are coincident with the prototype
// surface (with same or opposite normal, as specified), and copy all of
// their trim polygons into el. The edges are returned in uv coordinates for
// the prototype surface.
//-----------------------------------------------------------------------------
void SShell::MakeCoincidentEdgesInto(SSurface *proto, bool sameNormal,
                                     SEdgeList *el, SShell *useCurvesFrom)
{
    SSurface *ss;
    for(ss = surface.First(); ss; ss = surface.NextAfter(ss)) {
        if(proto->CoincidentWith(ss, sameNormal)) {
            ss->MakeEdgesInto(this, el, SSurface::MakeAs::XYZ, useCurvesFrom);
        }
    }

    SEdge *se;
    for(se = el->l.First(); se; se = el->l.NextAfter(se)) {
        double ua, va, ub, vb;
        proto->ClosestPointTo(se->a, &ua, &va);
        proto->ClosestPointTo(se->b, &ub, &vb);

        if(sameNormal) {
            se->a = Vector::From(ua, va, 0);
            se->b = Vector::From(ub, vb, 0);
        } else {
            // Flip normal, so flip all edge directions
            se->b = Vector::From(ua, va, 0);
            se->a = Vector::From(ub, vb, 0);
        }
    }
}

void SShell::MakeSectionEdgesInto(Vector n, double d, SEdgeList *sel, SBezierList *sbl)
{
    SSurface *s;
    for(s = surface.First(); s; s = surface.NextAfter(s)) {
        if(s->CoincidentWithPlane(n, d)) {
            s->MakeSectionEdgesInto(this, sel, sbl);
        }
    }
}

void SShell::AllPointsIntersecting(Vector a, Vector b,
                                   List<SInter> *il,
                                   bool seg, bool trimmed, bool inclTangent)
{
    SSurface *ss;
    for(ss = surface.First(); ss; ss = surface.NextAfter(ss)) {
        ss->AllPointsIntersecting(a, b, il, seg, trimmed, inclTangent);
    }
}

void SShell::Clear(void) {
    SSurface *s;
    for(s = surface.First(); s; s = surface.NextAfter(s)) {
        s->Clear();
    }
    surface.Clear();

    SCurve *c;
    for(c = curve.First(); c; c = curve.NextAfter(c)) {
        c->Clear();
    }
    curve.Clear();
}

void SSurface::EdgeNormalsWithinSurface(Point2d auv, Point2d buv,
                                        Vector *pt,
                                        Vector *enin, Vector *enout,
                                        Vector *surfn,
                                        uint32_t auxA,
                                        SShell *shell, SShell *sha, SShell *shb)
{
    // the midpoint of the edge
    Point2d muv  = (auv.Plus(buv)).ScaledBy(0.5);

    *pt    = PointAt(muv);

    // If this edge just approximates a curve, then refine our midpoint so
    // so that it actually lies on that curve too. Otherwise stuff like
    // point-on-face tests will fail, since the point won't actually lie
    // on the other face.
    hSCurve hc = { auxA };
    SCurve *sc = shell->curve.FindById(hc);
    if(sc->isExact && sc->exact.deg != 1) {
        double t;
        sc->exact.ClosestPointTo(*pt, &t, false);
        *pt = sc->exact.PointAt(t);
        ClosestPointTo(*pt, &muv);
    } else if(!sc->isExact) {
        SSurface *trimmedA = sc->GetSurfaceA(sha, shb),
                 *trimmedB = sc->GetSurfaceB(sha, shb);
        *pt = trimmedA->ClosestPointOnThisAndSurface(trimmedB, *pt);
        ClosestPointTo(*pt, &muv);
    }

    *surfn = NormalAt(muv.x, muv.y);

    // Compute the edge's inner normal in xyz space.
    Vector ab    = (PointAt(auv)).Minus(PointAt(buv)),
           enxyz = (ab.Cross(*surfn)).WithMagnitude(SS.ChordTolMm());
    // And based on that, compute the edge's inner normal in uv space. This
    // vector is perpendicular to the edge in xyz, but not necessarily in uv.
    Vector tu, tv;
    TangentsAt(muv.x, muv.y, &tu, &tv);
    Point2d enuv;
    enuv.x = enxyz.Dot(tu) / tu.MagSquared();
    enuv.y = enxyz.Dot(tv) / tv.MagSquared();

    // Compute the inner and outer normals of this edge (within the srf),
    // in xyz space. These are not necessarily antiparallel, if the
    // surface is curved.
    Vector pin   = PointAt(muv.Minus(enuv)),
           pout  = PointAt(muv.Plus(enuv));
    *enin  = pin.Minus(*pt),
    *enout = pout.Minus(*pt);
}

//-----------------------------------------------------------------------------
// When we split line segments wherever they intersect a surface, we introduce
// extra pwl points. This may create very short edges that could be removed
// without violating the chord tolerance. Those are ugly, and also break
// stuff in the Booleans. So remove them.
//-----------------------------------------------------------------------------
void SCurve::RemoveShortSegments(SSurface *srfA, SSurface *srfB) {
    // 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();
}

SSurface SSurface::FromTransformationOf(SSurface *a,
                                        Vector t, Quaternion q, double scale,
                                        bool includingTrims)
{
    SSurface ret = {};

    ret.h = a->h;
    ret.color = a->color;
    ret.face = a->face;

    ret.degm = a->degm;
    ret.degn = a->degn;
    int i, j;
    for(i = 0; i <= 3; i++) {
        for(j = 0; j <= 3; j++) {
            ret.ctrl[i][j] = a->ctrl[i][j];
            ret.ctrl[i][j] = (ret.ctrl[i][j]).ScaledBy(scale);
            ret.ctrl[i][j] = (q.Rotate(ret.ctrl[i][j])).Plus(t);

            ret.weight[i][j] = a->weight[i][j];
        }
    }

    if(includingTrims) {
        STrimBy *stb;
        for(stb = a->trim.First(); stb; stb = a->trim.NextAfter(stb)) {
            STrimBy n = *stb;
            n.start  = n.start.ScaledBy(scale);
            n.finish = n.finish.ScaledBy(scale);
            n.start  = (q.Rotate(n.start)) .Plus(t);
            n.finish = (q.Rotate(n.finish)).Plus(t);
            ret.trim.Add(&n);
        }
    }

    if(scale < 0) {
        // If we mirror every surface of a shell, then it will end up inside
        // out. So fix that here.
        ret.Reverse();
    }

    return ret;
}

void SShell::TriangulateInto(SMesh *sm) {
    SSurface *s;
    for(s = surface.First(); s; s = surface.NextAfter(s)) {
        s->TriangulateInto(this, sm);
    }
}

void SShell::MakeEdgesInto(SEdgeList *sel) {
    SSurface *s;
    for(s = surface.First(); s; s = surface.NextAfter(s)) {
        s->MakeEdgesInto(this, sel, SSurface::AS_XYZ);
    }
}

void SShell::MakeFromRevolutionOf(SBezierLoopSet *sbls, Vector pt, Vector axis, RgbaColor color, Group *group)
{
    SBezierLoop *sbl;

    int i0 = surface.n, i;

    // Normalize the axis direction so that the direction of revolution
    // ends up parallel to the normal of the sketch, on the side of the
    // axis where the sketch is.
    Vector pto;
    double md = VERY_NEGATIVE;
    for(sbl = sbls->l.First(); sbl; sbl = sbls->l.NextAfter(sbl)) {
        SBezier *sb;
        for(sb = sbl->l.First(); sb; sb = sbl->l.NextAfter(sb)) {
            // Choose the point farthest from the axis; we'll get garbage
            // if we choose a point that lies on the axis, for example.
            // (And our surface will be self-intersecting if the sketch
            // spans the axis, so don't worry about that.)
            Vector p = sb->Start();
            double d = p.DistanceToLine(pt, axis);
            if(d > md) {
                md = d;
                pto = p;
            }
        }
    }
    Vector ptc = pto.ClosestPointOnLine(pt, axis),
           up  = (pto.Minus(ptc)).WithMagnitude(1),
           vp  = (sbls->normal).Cross(up);
    if(vp.Dot(axis) < 0) {
        axis = axis.ScaledBy(-1);
    }

    // Now we actually build and trim the surfaces.
    for(sbl = sbls->l.First(); sbl; sbl = sbls->l.NextAfter(sbl)) {
        int i, j;
        SBezier *sb, *prev;
        List<Revolved> hsl = {};

        for(sb = sbl->l.First(); sb; sb = sbl->l.NextAfter(sb)) {
            Revolved revs;
            for(j = 0; j < 4; j++) {
                if(sb->deg == 1 &&
                        (sb->ctrl[0]).DistanceToLine(pt, axis) < LENGTH_EPS &&
                        (sb->ctrl[1]).DistanceToLine(pt, axis) < LENGTH_EPS)
                {
                    // This is a line on the axis of revolution; it does
                    // not contribute a surface.
                    revs.d[j].v = 0;
                } else {
                    SSurface ss = SSurface::FromRevolutionOf(sb, pt, axis,
                                  (PI/2)*j,
                                  (PI/2)*(j+1));
                    ss.color = color;
                    if(sb->entity != 0) {
                        hEntity he;
                        he.v = sb->entity;
                        hEntity hface = group->Remap(he, Group::REMAP_LINE_TO_FACE);
                        if(SK.entity.FindByIdNoOops(hface) != NULL) {
                            ss.face = hface.v;
                        }
                    }
                    revs.d[j] = surface.AddAndAssignId(&ss);
                }
            }
            hsl.Add(&revs);
        }

        for(i = 0; i < sbl->l.n; i++) {
            Revolved revs  = hsl.elem[i],
                     revsp = hsl.elem[WRAP(i-1, sbl->l.n)];

            sb   = &(sbl->l.elem[i]);
            prev = &(sbl->l.elem[WRAP(i-1, sbl->l.n)]);

            for(j = 0; j < 4; j++) {
                SCurve sc;
                Quaternion qs = Quaternion::From(axis, (PI/2)*j);
                // we want Q*(x - p) + p = Q*x + (p - Q*p)
                Vector ts = pt.Minus(qs.Rotate(pt));

                // If this input curve generate a surface, then trim that
                // surface with the rotated version of the input curve.
                if(revs.d[j].v) {
                    sc = {};
                    sc.isExact = true;
                    sc.exact = sb->TransformedBy(ts, qs, 1.0);
                    (sc.exact).MakePwlInto(&(sc.pts));
                    sc.surfA = revs.d[j];
                    sc.surfB = revs.d[WRAP(j-1, 4)];

                    hSCurve hcb = curve.AddAndAssignId(&sc);

                    STrimBy stb;
                    stb = STrimBy::EntireCurve(this, hcb, true);
                    (surface.FindById(sc.surfA))->trim.Add(&stb);
                    stb = STrimBy::EntireCurve(this, hcb, false);
                    (surface.FindById(sc.surfB))->trim.Add(&stb);
                }

//.........这里部分代码省略.........

//-----------------------------------------------------------------------------
// Trim this surface against the specified shell, in the way that's appropriate
// for the specified Boolean operation type (and which operand we are). We
// also need a pointer to the shell that contains our own surface, since that
// contains our original trim curves.
//-----------------------------------------------------------------------------
SSurface SSurface::MakeCopyTrimAgainst(SShell *parent,
                                       SShell *sha, SShell *shb,
                                       SShell *into,
                                       int type)
{
    bool opA = (parent == sha);
    SShell *agnst = opA ? shb : sha;

    SSurface ret;
    // The returned surface is identical, just the trim curves change
    ret = *this;
    ret.trim = {};

    // First, build a list of the existing trim curves; update them to use
    // the split curves.
    STrimBy *stb;
    for(stb = trim.First(); stb; stb = trim.NextAfter(stb)) {
        STrimBy stn = *stb;
        stn.curve = (parent->curve.FindById(stn.curve))->newH;
        ret.trim.Add(&stn);
    }

    if(type == SShell::AS_DIFFERENCE && !opA) {
        // The second operand of a Boolean difference gets turned inside out
        ret.Reverse();
    }

    // Build up our original trim polygon; remember the coordinates could
    // be changed if we just flipped the surface normal, and we are using
    // the split curves (not the original curves).
    SEdgeList orig = {};
    ret.MakeEdgesInto(into, &orig, AS_UV);
    ret.trim.Clear();
    // which means that we can't necessarily use the old BSP...
    SBspUv *origBsp = SBspUv::From(&orig, &ret);

    // And now intersect the other shell against us
    SEdgeList inter = {};

    SSurface *ss;
    for(ss = agnst->surface.First(); ss; ss = agnst->surface.NextAfter(ss)) {
        SCurve *sc;
        for(sc = into->curve.First(); sc; sc = into->curve.NextAfter(sc)) {
            if(sc->source != SCurve::FROM_INTERSECTION) continue;
            if(opA) {
                if(sc->surfA.v != h.v || sc->surfB.v != ss->h.v) continue;
            } else {
                if(sc->surfB.v != h.v || sc->surfA.v != ss->h.v) continue;
            }

            int i;
            for(i = 1; i < sc->pts.n; i++) {
                Vector a = sc->pts.elem[i-1].p,
                       b = sc->pts.elem[i].p;

                Point2d auv, buv;
                ss->ClosestPointTo(a, &(auv.x), &(auv.y));
                ss->ClosestPointTo(b, &(buv.x), &(buv.y));

                int c = (ss->bsp) ? ss->bsp->ClassifyEdge(auv, buv, ss) : SBspUv::OUTSIDE;
                if(c != SBspUv::OUTSIDE) {
                    Vector ta = Vector::From(0, 0, 0);
                    Vector tb = Vector::From(0, 0, 0);
                    ret.ClosestPointTo(a, &(ta.x), &(ta.y));
                    ret.ClosestPointTo(b, &(tb.x), &(tb.y));

                    Vector tn = ret.NormalAt(ta.x, ta.y);
                    Vector sn = ss->NormalAt(auv.x, auv.y);

                    // We are subtracting the portion of our surface that
                    // lies in the shell, so the in-plane edge normal should
                    // point opposite to the surface normal.
                    bool bkwds = true;
                    if((tn.Cross(b.Minus(a))).Dot(sn) < 0) bkwds = !bkwds;
                    if(type == SShell::AS_DIFFERENCE && !opA) bkwds = !bkwds;
                    if(bkwds) {
                        inter.AddEdge(tb, ta, sc->h.v, 1);
                    } else {
                        inter.AddEdge(ta, tb, sc->h.v, 0);
                    }
                }
            }
        }
    }

    // Record all the points where more than two edges join, which I will call
    // the choosing points. If two edges join at a non-choosing point, then
    // they must either both be kept or both be discarded (since that would
    // otherwise create an open contour).
    SPointList choosing = {};
    SEdge *se;
    for(se = orig.l.First(); se; se = orig.l.NextAfter(se)) {
        choosing.IncrementTagFor(se->a);
        choosing.IncrementTagFor(se->b);
//.........这里部分代码省略.........

void Group::GenerateShellAndMesh(void) {
    bool prevBooleanFailed = booleanFailed;
    booleanFailed = false;

    Group *srcg = this;

    thisShell.Clear();
    thisMesh.Clear();
    runningShell.Clear();
    runningMesh.Clear();

    // Don't attempt a lathe or extrusion unless the source section is good:
    // planar and not self-intersecting.
    bool haveSrc = true;
    if(type == EXTRUDE || type == LATHE) {
        Group *src = SK.GetGroup(opA);
        if(src->polyError.how != POLY_GOOD) {
            haveSrc = false;
        }
    }

    if(type == TRANSLATE || type == ROTATE) {
        // A step and repeat gets merged against the group's prevous group,
        // not our own previous group.
        srcg = SK.GetGroup(opA);

        GenerateForStepAndRepeat<SShell>(&(srcg->thisShell), &thisShell);
        GenerateForStepAndRepeat<SMesh> (&(srcg->thisMesh),  &thisMesh);
    } else if(type == EXTRUDE && haveSrc) {
        Group *src = SK.GetGroup(opA);
        Vector translate = Vector::From(h.param(0), h.param(1), h.param(2));

        Vector tbot, ttop;
        if(subtype == ONE_SIDED) {
            tbot = Vector::From(0, 0, 0); ttop = translate.ScaledBy(2);
        } else {
            tbot = translate.ScaledBy(-1); ttop = translate.ScaledBy(1);
        }

        SBezierLoopSetSet *sblss = &(src->bezierLoops);
        SBezierLoopSet *sbls;
        for(sbls = sblss->l.First(); sbls; sbls = sblss->l.NextAfter(sbls)) {
            int is = thisShell.surface.n;
            // Extrude this outer contour (plus its inner contours, if present)
            thisShell.MakeFromExtrusionOf(sbls, tbot, ttop, color);

            // And for any plane faces, annotate the model with the entity for
            // that face, so that the user can select them with the mouse.
            Vector onOrig = sbls->point;
            int i;
            for(i = is; i < thisShell.surface.n; i++) {
                SSurface *ss = &(thisShell.surface.elem[i]);
                hEntity face = Entity::NO_ENTITY;

                Vector p = ss->PointAt(0, 0),
                       n = ss->NormalAt(0, 0).WithMagnitude(1);
                double d = n.Dot(p);

                if(i == is || i == (is + 1)) {
                    // These are the top and bottom of the shell.
                    if(fabs((onOrig.Plus(ttop)).Dot(n) - d) < LENGTH_EPS) {
                        face = Remap(Entity::NO_ENTITY, REMAP_TOP);
                        ss->face = face.v;
                    }
                    if(fabs((onOrig.Plus(tbot)).Dot(n) - d) < LENGTH_EPS) {
                        face = Remap(Entity::NO_ENTITY, REMAP_BOTTOM);
                        ss->face = face.v;
                    }
                    continue;
                }

                // So these are the sides
                if(ss->degm != 1 || ss->degn != 1) continue;

                Entity *e;
                for(e = SK.entity.First(); e; e = SK.entity.NextAfter(e)) {
                    if(e->group.v != opA.v) continue;
                    if(e->type != Entity::LINE_SEGMENT) continue;

                    Vector a = SK.GetEntity(e->point[0])->PointGetNum(),
                           b = SK.GetEntity(e->point[1])->PointGetNum();
                    a = a.Plus(ttop);
                    b = b.Plus(ttop);
                    // Could get taken backwards, so check all cases.
                    if((a.Equals(ss->ctrl[0][0]) && b.Equals(ss->ctrl[1][0])) ||
                       (b.Equals(ss->ctrl[0][0]) && a.Equals(ss->ctrl[1][0])) ||
                       (a.Equals(ss->ctrl[0][1]) && b.Equals(ss->ctrl[1][1])) ||
                       (b.Equals(ss->ctrl[0][1]) && a.Equals(ss->ctrl[1][1])))
                    {
                        face = Remap(e->h, REMAP_LINE_TO_FACE);
                        ss->face = face.v;
                        break;
                    }
                }
            }
        }
    } else if(type == LATHE && haveSrc) {
        Group *src = SK.GetGroup(opA);

        Vector pt   = SK.GetEntity(predef.origin)->PointGetNum(),
//.........这里部分代码省略.........

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
//.........这里部分代码省略.........

//-----------------------------------------------------------------------------
// Does the given point lie on our shell? There are many cases; inside and
// outside are obvious, but then there's all the edge-on-edge and edge-on-face
// possibilities.
//
// To calculate, we intersect a ray through p with our shell, and classify
// using the closest intersection point. If the ray hits a surface on edge,
// then just reattempt in a different random direction.
//-----------------------------------------------------------------------------
bool SShell::ClassifyEdge(int *indir, int *outdir,
                          Vector ea, Vector eb,
                          Vector p,
                          Vector edge_n_in, Vector edge_n_out, Vector surf_n)
{
    List<SInter> l;
    ZERO(&l);

    srand(0);

    // First, check for edge-on-edge
    int edge_inters = 0;
    Vector inter_surf_n[2], inter_edge_n[2];
    SSurface *srf;
    for(srf = surface.First(); srf; srf = surface.NextAfter(srf)) {
        if(srf->LineEntirelyOutsideBbox(ea, eb, true)) continue;

        SEdgeList *sel = &(srf->edges);
        SEdge *se;
        for(se = sel->l.First(); se; se = sel->l.NextAfter(se)) {
            if((ea.Equals(se->a) && eb.Equals(se->b)) ||
               (eb.Equals(se->a) && ea.Equals(se->b)) ||
                p.OnLineSegment(se->a, se->b))
            {
                if(edge_inters < 2) {
                    // Edge-on-edge case
                    Point2d pm;
                    srf->ClosestPointTo(p,  &pm, false);
                    // A vector normal to the surface, at the intersection point
                    inter_surf_n[edge_inters] = srf->NormalAt(pm);
                    // A vector normal to the intersecting edge (but within the
                    // intersecting surface) at the intersection point, pointing
                    // out.
                    inter_edge_n[edge_inters] =
                      (inter_surf_n[edge_inters]).Cross((se->b).Minus((se->a)));
                }

                edge_inters++;
            }
        }
    }

    if(edge_inters == 2) {
        // TODO, make this use the appropriate curved normals
        double dotp[2];
        for(int i = 0; i < 2; i++) {
            dotp[i] = edge_n_out.DirectionCosineWith(inter_surf_n[i]);
        }

        if(fabs(dotp[1]) < DOTP_TOL) {
            SWAP(double, dotp[0],         dotp[1]);
            SWAP(Vector, inter_surf_n[0], inter_surf_n[1]);
            SWAP(Vector, inter_edge_n[0], inter_edge_n[1]);
        }

        int coinc = (surf_n.Dot(inter_surf_n[0])) > 0 ? COINC_SAME : COINC_OPP;

        if(fabs(dotp[0]) < DOTP_TOL && fabs(dotp[1]) < DOTP_TOL) {
            // This is actually an edge on face case, just that the face
            // is split into two pieces joining at our edge.
            *indir  = coinc;
            *outdir = coinc;
        } else if(fabs(dotp[0]) < DOTP_TOL && dotp[1] > DOTP_TOL) {
            if(edge_n_out.Dot(inter_edge_n[0]) > 0) {
                *indir  = coinc;
                *outdir = OUTSIDE;
            } else {
                *indir  = INSIDE;
                *outdir = coinc;
            }
        } else if(fabs(dotp[0]) < DOTP_TOL && dotp[1] < -DOTP_TOL) {
            if(edge_n_out.Dot(inter_edge_n[0]) > 0) {
                *indir  = coinc;
                *outdir = INSIDE;
            } else {
                *indir  = OUTSIDE;
                *outdir = coinc;
            }
        } else if(dotp[0] > DOTP_TOL && dotp[1] > DOTP_TOL) {
            *indir  = INSIDE;
            *outdir = OUTSIDE;
        } else if(dotp[0] < -DOTP_TOL && dotp[1] < -DOTP_TOL) {
            *indir  = OUTSIDE;
            *outdir = INSIDE;
        } else {
            // Edge is tangent to the shell at shell's edge, so can't be
            // a boundary of the surface.
            return false;
        }
        return true;
    }
//.........这里部分代码省略.........

void StepFileWriter::ExportSurfacesTo(char *file) {
    Group *g = SK.GetGroup(SS.GW.activeGroup);
    SShell *shell = &(g->runningShell);

    if(shell->surface.n == 0) {
        Error("The model does not contain any surfaces to export.%s",
            g->runningMesh.l.n > 0 ?
                "\n\nThe model does contain triangles from a mesh, but "
                "a triangle mesh cannot be exported as a STEP file. Try "
                "File -> Export Mesh... instead." : "");
        return;
    }

    f = fopen(file, "wb");
    if(!f) {
        Error("Couldn't write to '%s'", file);
        return;
    }

    WriteHeader();
	WriteProductHeader();

    ZERO(&advancedFaces);

    SSurface *ss;
    for(ss = shell->surface.First(); ss; ss = shell->surface.NextAfter(ss)) {
        if(ss->trim.n == 0) continue;

        // Get all of the loops of Beziers that trim our surface (with each
        // Bezier split so that we use the section as t goes from 0 to 1), and
        // the piecewise linearization of those loops in xyz space.
        SBezierList sbl;
        ZERO(&sbl);
        ss->MakeSectionEdgesInto(shell, NULL, &sbl);

        // Apply the export scale factor.
        ss->ScaleSelfBy(1.0/SS.exportScale);
        sbl.ScaleSelfBy(1.0/SS.exportScale);

        ExportSurface(ss, &sbl);

        sbl.Clear();
    }

    fprintf(f, "#%d=CLOSED_SHELL('',(", id);
    int *af;
    for(af = advancedFaces.First(); af; af = advancedFaces.NextAfter(af)) {
        fprintf(f, "#%d", *af);
        if(advancedFaces.NextAfter(af) != NULL) fprintf(f, ",");
    }
    fprintf(f, "));\n");
    fprintf(f, "#%d=MANIFOLD_SOLID_BREP('brep',#%d);\n", id+1, id);
    fprintf(f, "#%d=ADVANCED_BREP_SHAPE_REPRESENTATION('',(#%d,#170),#168);\n",
        id+2, id+1);
    fprintf(f, "#%d=SHAPE_REPRESENTATION_RELATIONSHIP($,$,#169,#%d);\n",
        id+3, id+2);

    WriteFooter();

    fclose(f);
    advancedFaces.Clear();
}