Replace the closed-form solutions for entity-entity splitting with

a method that works on the piecewise linear segments, and then
refines any intersections that it finds by Newton's method. So now
I support cubics too, and circle-circle intersections, and the code
is much simpler.

[git-p4: depot-paths = "//depot/solvespace/": change = 2012]

Makefile

@@ -59,7 +59,7 @@ LIBS = user32.lib gdi32.lib comctl32.lib advapi32.lib shell32.lib opengl32.lib g
 
 all: $(OBJDIR)/solvespace.exe
     @cp $(OBJDIR)/solvespace.exe .
-    solvespace block.slvs
+    solvespace t8.slvs
 
 clean:
 	rm -f obj/*

dsc.h

@@ -59,6 +59,9 @@ public:
     static Vector AtIntersectionOfPlanes(Vector na, double da,
                                          Vector nb, double db,
                                          Vector nc, double dc, bool *parallel);
+    static void ClosestPointBetweenLines(Vector pa, Vector da,
+                                         Vector pb, Vector db,
+                                         double *ta, double *tb);
 
     double Element(int i);
     bool Equals(Vector v, double tol=LENGTH_EPS);

entity.cpp

@@ -680,11 +680,27 @@ void EntityBase::GenerateEquations(IdList<Equation,hEquation> *l) {
             // If this is a copied entity, with its point already fixed
             // with respect to each other, then we don't want to generate
             // the distance constraint!
-            if(SK.GetEntity(point[0])->type == POINT_IN_2D) {
+            if(SK.GetEntity(point[0])->type != POINT_IN_2D) break;
+
+            // If the two endpoints of the arc are constrained coincident
+            // (to make a complete circle), then our distance constraint
+            // would be redundant and therefore overconstrain things.
+            Constraint *c;
+            for(c = SK.constraint.First(); c; c = SK.constraint.NextAfter(c)) {
+                if(c->group.v != group.v) continue;
+                if(c->type != Constraint::POINTS_COINCIDENT) continue;
+                
+                if((c->ptA.v == point[1].v && c->ptB.v == point[2].v) || 
+                   (c->ptA.v == point[2].v && c->ptB.v == point[1].v))
+                {
+                    break;
+                }
+            }
+            if(c) break;
+
             Expr *ra = Constraint::Distance(workplane, point[0], point[1]);
             Expr *rb = Constraint::Distance(workplane, point[0], point[2]);
             AddEq(l, ra->Minus(rb), 0);
-            }
             break;
         }
         default:;

modify.cpp

@@ -157,7 +157,7 @@ void GraphicsWindow::MakeTangentArc(void) {
     SS.later.generateAll = true;
 }
 
-void GraphicsWindow::SplitLine(hEntity he, Vector pinter) {
+hEntity GraphicsWindow::SplitLine(hEntity he, Vector pinter) {
     // Save the original endpoints, since we're about to delete this entity.
     Entity *e01 = SK.GetEntity(he);
     hEntity hep0 = e01->point[0], hep1 = e01->point[1];
@@ -181,26 +181,11 @@ void GraphicsWindow::SplitLine(hEntity he, Vector pinter) {
 
     ReplacePointInConstraints(hep0, e0i->point[0]);
     ReplacePointInConstraints(hep1, ei1->point[1]);
-
+    Constraint::ConstrainCoincident(e0i->point[1], ei1->point[0]);
-    // Finally, delete the original line
+    return e0i->point[1];
-    int i;
-    SK.request.ClearTags();
-    for(i = 0; i < SK.request.n; i++) {
-        Request *r = &(SK.request.elem[i]);
-        if(r->group.v != activeGroup.v) continue;
-        if(r->type != Request::LINE_SEGMENT) continue;
-    
-        // If the user wants to keep the old lines around, they can just
-        // mark it construction first.
-        if(he.v == r->h.entity(0).v && !r->construction) {
-            r->tag = 1;
-            break;
-        }
-    }
-    DeleteTaggedRequests();
 }
 
-void GraphicsWindow::SplitCircle(hEntity he, Vector pinter) {
+hEntity GraphicsWindow::SplitCircle(hEntity he, Vector pinter) {
     SS.UndoRemember();
 
     Entity *circle = SK.GetEntity(he);
@@ -216,11 +201,17 @@ void GraphicsWindow::SplitCircle(hEntity he, Vector pinter) {
         SK.GetEntity(arc->point[0])->PointForceTo(center);
         SK.GetEntity(arc->point[1])->PointForceTo(pinter);
         SK.GetEntity(arc->point[2])->PointForceTo(pinter);
+
+        Constraint::ConstrainCoincident(arc->point[1], arc->point[2]);
+        return arc->point[1];
     } else {
         // Start with an arc, break it in to two arcs
-        Vector center = SK.GetEntity(circle->point[0])->PointGetNum(),
+        hEntity hc = circle->point[0],
-               start  = SK.GetEntity(circle->point[1])->PointGetNum(),
+                hs = circle->point[1],
-               finish = SK.GetEntity(circle->point[2])->PointGetNum();
+                hf = circle->point[2];
+        Vector center = SK.GetEntity(hc)->PointGetNum(),
+               start  = SK.GetEntity(hs)->PointGetNum(),
+               finish = SK.GetEntity(hf)->PointGetNum();
 
         circle = NULL; // shortly invalid!
         hRequest hr0 = AddRequest(Request::ARC_OF_CIRCLE, false),
@@ -236,26 +227,99 @@ void GraphicsWindow::SplitCircle(hEntity he, Vector pinter) {
         SK.GetEntity(arc1->point[0])->PointForceTo(center);
         SK.GetEntity(arc1->point[1])->PointForceTo(pinter);
         SK.GetEntity(arc1->point[2])->PointForceTo(finish);
+
+        ReplacePointInConstraints(hs, arc0->point[1]);
+        ReplacePointInConstraints(hf, arc1->point[2]);
+        Constraint::ConstrainCoincident(arc0->point[2], arc1->point[1]);
+        return arc0->point[2];
+    }
 }
 
-    // Finally, delete the original circle or arc
+hEntity GraphicsWindow::SplitCubic(hEntity he, Vector pinter) {
+    // Save the original endpoints, since we're about to delete this entity.
+    Entity *e01 = SK.GetEntity(he);
+    hEntity hep0 = e01->point[0],
+            hep1 = e01->point[1],
+            hep2 = e01->point[2],
+            hep3 = e01->point[3];
+    Vector p0 = SK.GetEntity(hep0)->PointGetNum(),
+           p1 = SK.GetEntity(hep1)->PointGetNum(),
+           p2 = SK.GetEntity(hep2)->PointGetNum(),
+           p3 = SK.GetEntity(hep3)->PointGetNum();
+
+    SS.UndoRemember();
+
+    SBezier b0i, bi1, b01 = SBezier::From(p0, p1, p2, p3);
+    double t;
+    b01.ClosestPointTo(pinter, &t, true);
+    b01.SplitAt(t, &b0i, &bi1);
+
+    // Add the two line segments this one gets split into.
+    hRequest r0i = AddRequest(Request::CUBIC, false),
+             ri1 = AddRequest(Request::CUBIC, false);
+    // Don't get entities till after adding, realloc issues
+
+    Entity *e0i = SK.GetEntity(r0i.entity(0)),
+           *ei1 = SK.GetEntity(ri1.entity(0));
+
+    SK.GetEntity(e0i->point[0])->PointForceTo(b0i.ctrl[0]);
+    SK.GetEntity(e0i->point[1])->PointForceTo(b0i.ctrl[1]);
+    SK.GetEntity(e0i->point[2])->PointForceTo(b0i.ctrl[2]);
+    SK.GetEntity(e0i->point[3])->PointForceTo(b0i.ctrl[3]);
+
+    SK.GetEntity(ei1->point[0])->PointForceTo(bi1.ctrl[0]);
+    SK.GetEntity(ei1->point[1])->PointForceTo(bi1.ctrl[1]);
+    SK.GetEntity(ei1->point[2])->PointForceTo(bi1.ctrl[2]);
+    SK.GetEntity(ei1->point[3])->PointForceTo(bi1.ctrl[3]);
+
+    ReplacePointInConstraints(hep0, e0i->point[0]);
+    ReplacePointInConstraints(hep1, ei1->point[3]);
+    Constraint::ConstrainCoincident(e0i->point[3], ei1->point[0]);
+    return e0i->point[3];
+}
+
+hEntity GraphicsWindow::SplitEntity(hEntity he, Vector pinter) {
+    Entity *e = SK.GetEntity(he);
+    int entityType = e->type;
+
+    hEntity ret;
+    if(e->IsCircle()) {
+        ret = SplitCircle(he, pinter);
+    } else if(e->type == Entity::LINE_SEGMENT) {
+        ret = SplitLine(he, pinter);
+    } else if(e->type == Entity::CUBIC) {
+        ret = SplitCubic(he, pinter);
+    } else {
+        Error("Couldn't split this entity; lines, circles, or cubics only.");
+        return Entity::NO_ENTITY;
+    }
+
+    // Finally, delete the request that generated the original entity.
+    int reqType;
+    switch(entityType) {
+        case Entity::CIRCLE:        reqType = Request::CIRCLE; break;
+        case Entity::ARC_OF_CIRCLE: reqType = Request::ARC_OF_CIRCLE; break;
+        case Entity::LINE_SEGMENT:  reqType = Request::LINE_SEGMENT; break;
+        case Entity::CUBIC:         reqType = Request::CUBIC; break;
+        default: oops();
+    }
     int i;
     SK.request.ClearTags();
     for(i = 0; i < SK.request.n; i++) {
         Request *r = &(SK.request.elem[i]);
         if(r->group.v != activeGroup.v) continue;
-        if(r->type != Request::CIRCLE && r->type != Request::ARC_OF_CIRCLE) {
+        if(r->type != reqType) continue;
-            continue;
-        }
     
-        // If the user wants to keep the old lines around, they can just
+        // If the user wants to keep the old entities around, they can just
-        // mark it construction first.
+        // mark them construction first.
         if(he.v == r->h.entity(0).v && !r->construction) {
             r->tag = 1;
             break;
         }
     }
     DeleteTaggedRequests();
+
+    return ret;
 }
 
 void GraphicsWindow::SplitLinesOrCurves(void) {
@@ -265,144 +329,47 @@ void GraphicsWindow::SplitLinesOrCurves(void) {
     }
 
     GroupSelection();
-    if(gs.n == 2 && gs.lineSegments == 2) {
+    if(!(gs.n == 2 && (gs.lineSegments + gs.circlesOrArcs + gs.cubics) == 2)) {
-        Entity *la = SK.GetEntity(gs.entity[0]),
+        Error("Select two entities that intersect each other (e.g. two lines "
-               *lb = SK.GetEntity(gs.entity[1]);
+              "or two circles or a circle and a line).");
-        Vector a0 = SK.GetEntity(la->point[0])->PointGetNum(),
+        return;
-               a1 = SK.GetEntity(la->point[1])->PointGetNum(),
-               b0 = SK.GetEntity(lb->point[0])->PointGetNum(),
-               b1 = SK.GetEntity(lb->point[1])->PointGetNum();
-        Vector da = a1.Minus(a0), db = b1.Minus(b0);
-        
-        // First, check if the lines intersect.
-        bool skew;
-        Vector pinter = Vector::AtIntersectionOfLines(a0, a1, b0, b1, &skew);
-        if(skew) {
-            Error("Lines are parallel or skew; no intersection to split.");
-            goto done;
     }
 
-        double ta = (pinter.Minus(a0)).DivPivoting(da),
+    hEntity ha = gs.entity[0],
-               tb = (pinter.Minus(b0)).DivPivoting(db);
+            hb = gs.entity[1];
+    Entity *ea = SK.GetEntity(ha),
+           *eb = SK.GetEntity(hb);
    
-        double tola = LENGTH_EPS/da.Magnitude(),
+    // Compute the possibly-rational Bezier curves for each of these entities
-               tolb = LENGTH_EPS/db.Magnitude();
+    SBezierList sbla, sblb;
+    ZERO(&sbla);
+    ZERO(&sblb);
+    ea->GenerateBezierCurves(&sbla);
+    eb->GenerateBezierCurves(&sblb);
+    // and then compute the points where they intersect, based on those curves.
+    SPointList inters;
+    ZERO(&inters);
+    sbla.AllIntersectionsWith(&sblb, &inters);
 
-        hEntity ha = la->h, hb = lb->h;
+    // If there's multiple points, then just take the first one.
-        la = NULL; lb = NULL;
+    if(inters.l.n > 0) {
-        // Following adds will cause a realloc, break the pointers
+        Vector pi = inters.l.elem[0].p;
-
+        hEntity hia = SplitEntity(ha, pi),
-        bool didSomething = false;
+                hib = SplitEntity(hb, pi);
-        if(ta > tola && ta < (1 - tola)) {
+        // SplitEntity adds the coincident constraints to join the split halves
-            SplitLine(ha, pinter);
+        // of each original entity; and then we add the constraint to join
-            didSomething = true;
+        // the two entities together at the split point.
-        }
+        if(hia.v && hib.v) {
-        if(tb > tolb && tb < (1 - tolb)) {
+            Constraint::ConstrainCoincident(hia, hib);
-            SplitLine(hb, pinter);
-            didSomething = true;
-        }
-        if(!didSomething) {
-            Error(
-                "Nothing to split; intersection does not lie on either line.");
-        }
-    } else if(gs.n == 2 && gs.lineSegments == 1 && gs.circlesOrArcs == 1) {
-        Entity *line   = SK.GetEntity(gs.entity[0]),
-               *circle = SK.GetEntity(gs.entity[1]);
-        if(line->type != Entity::LINE_SEGMENT) {
-            SWAP(Entity *, line, circle);
-        }
-        hEntity hline = line->h, hcircle = circle->h;
-
-        Vector l0 = SK.GetEntity(line->point[0])->PointGetNum(),
-               l1 = SK.GetEntity(line->point[1])->PointGetNum();
-        Vector dl = l1.Minus(l0);
-
-        Quaternion q = SK.GetEntity(circle->normal)->NormalGetNum();
-        Vector cn = q.RotationN();
-        Vector cc = SK.GetEntity(circle->point[0])->PointGetNum();
-        double cd = cc.Dot(cn);
-        double cr = circle->CircleGetRadiusNum();
-
-        if(fabs(l0.Dot(cn) - cd) > LENGTH_EPS ||
-           fabs(l1.Dot(cn) - cd) > LENGTH_EPS)
-        {
-            Error("Lines does not lie in same plane as circle.");
-            goto done;
-        }
-        // Now let's see if they intersect; transform everything into a csys
-        // with origin at the center of the circle, and where the line is
-        // horizontal.
-        Vector n = cn.WithMagnitude(1);
-        Vector u = dl.WithMagnitude(1);
-        Vector v = n.Cross(u);
-        
-        Vector nl0 = (l0.Minus(cc)).DotInToCsys(u, v, n),
-               nl1 = (l1.Minus(cc)).DotInToCsys(u, v, n);
-
-        double yint = nl0.y;
-        if(fabs(yint) > (cr - LENGTH_EPS)) {
-            Error("Line does not intersect (or is tangent to) circle.");
-            goto done;
-        }
-        double xint = sqrt(cr*cr - yint*yint);
-
-        Vector inter0 = Vector::From( xint, yint, 0),
-               inter1 = Vector::From(-xint, yint, 0);
-
-        // While we're here, let's calculate the angles at which the
-        // intersections (and the endpoints of the arc) occur.
-        double theta0 = atan2(yint, xint), theta1 = atan2(yint, -xint);
-        double thetamin, thetamax;
-        if(circle->type == Entity::CIRCLE) {
-            thetamin = 0;
-            thetamax = 2.1*PI; // fudge, make sure it's a good complete circle
-        } else {
-            Vector start  = SK.GetEntity(circle->point[1])->PointGetNum(),
-                   finish = SK.GetEntity(circle->point[2])->PointGetNum();
-            start  = (start .Minus(cc)).DotInToCsys(u, v, n);
-            finish = (finish.Minus(cc)).DotInToCsys(u, v, n);
-            thetamin = atan2(start.y, start.x);
-            thetamax = atan2(finish.y, finish.x);
-
-            // Normalize; arc is drawn with increasing theta from start,
-            // so subtract that off and make all angles in (0, 2*pi]
-            theta0   = WRAP_NOT_0(theta0   - thetamin, 2*PI);
-            theta1   = WRAP_NOT_0(theta1   - thetamin, 2*PI);
-            thetamax = WRAP_NOT_0(thetamax - thetamin, 2*PI);
-        }
-
-        // And move our intersections back out to the base frame.
-        inter0 = inter0.ScaleOutOfCsys(u, v, n).Plus(cc);
-        inter1 = inter1.ScaleOutOfCsys(u, v, n).Plus(cc);
-
-        // So now we have our intersection points. Let's see where they are
-        // on the line.
-        double t0 = (inter0.Minus(l0)).DivPivoting(dl),
-               t1 = (inter1.Minus(l0)).DivPivoting(dl);
-        double tol = LENGTH_EPS/dl.Magnitude();
-   
-        bool didSomething = false;
-        // Split only once, even if it crosses multiple times; just pick
-        // arbitrarily which.
-        if(t0 > tol && t0 < (1 - tol) && theta0 < thetamax) {
-            SplitLine(hline, inter0);
-            SplitCircle(hcircle, inter0);
-            didSomething = true;
-        } else if(t1 > tol && t1 < (1 - tol) && theta1 < thetamax) {
-            SplitLine(hline, inter1);
-            SplitCircle(hcircle, inter1);
-            didSomething = true;
-        }
-        if(!didSomething) {
-            Error("Nothing to split; neither intersection lies on both the "
-                  "line and the circle.");
         }
     } else {
-        Error("Can't split these entities; select two lines, a line and "
+        Error("Can't split; no intersection found.");
-              "a circle, or a line and an arc.");
     }
 
-done:
+    // All done, clean up and regenerate.
+    inters.Clear();
+    sbla.Clear();
+    sblb.Clear();
     ClearSelection();
     SS.later.generateAll = true;
 }

polygon.cpp

@@ -148,7 +148,9 @@ bool SEdgeList::AssemblePolygon(SPolygon *dest, SEdge *errorAt, bool keepDir) {
 // but they are considered to cross if they are coincident and overlapping.
 // If pi is not NULL, then a crossing is returned in that.
 //-----------------------------------------------------------------------------
-int SEdgeList::AnyEdgeCrossings(Vector a, Vector b, Vector *ppi) {
+int SEdgeList::AnyEdgeCrossings(Vector a, Vector b, 
+                                Vector *ppi, SPointList *spl)
+{
     Vector d = b.Minus(a);
     double t_eps = LENGTH_EPS/d.Magnitude();
 
@@ -214,6 +216,7 @@ int SEdgeList::AnyEdgeCrossings(Vector a, Vector b, Vector *ppi) {
             // inside of the other (or if they cross away from either's
             // vertex).
             if(ppi) *ppi = pi;
+            if(spl) spl->Add(pi);
             goto intersects;
         }
         continue;

polygon.h

@@ -2,6 +2,7 @@
 #ifndef __POLYGON_H
 #define __POLYGON_H
 
+class SPointList;
 class SPolygon;
 class SContour;
 class SMesh;
@@ -25,7 +26,8 @@ public:
     bool AssemblePolygon(SPolygon *dest, SEdge *errorAt, bool keepDir=false);
     bool AssembleContour(Vector first, Vector last, SContour *dest,
                             SEdge *errorAt, bool keepDir);
-    int AnyEdgeCrossings(Vector a, Vector b, Vector *pi=NULL);
+    int AnyEdgeCrossings(Vector a, Vector b,
+        Vector *pi=NULL, SPointList *spl=NULL);
     bool ContainsEdgeFrom(SEdgeList *sel);
     bool ContainsEdge(SEdge *se);
     void CullExtraneousEdges(void);

srf/curve.cpp

@@ -220,6 +220,48 @@ void SBezierList::CullIdenticalBeziers(void) {
     l.RemoveTagged();
 }
 
+//-----------------------------------------------------------------------------
+// Find all the points where a list of Bezier curves intersects another list
+// of Bezier curves. We do this by intersecting their piecewise linearizations,
+// and then refining any intersections that we find to lie exactly on the
+// curves. So this will screw up on tangencies and stuff, but otherwise should
+// be fine.
+//-----------------------------------------------------------------------------
+void SBezierList::AllIntersectionsWith(SBezierList *sblb, SPointList *spl) {
+    SBezier *sba, *sbb;
+    for(sba = l.First(); sba; sba = l.NextAfter(sba)) {
+        for(sbb = sblb->l.First(); sbb; sbb = sblb->l.NextAfter(sbb)) {
+            sbb->AllIntersectionsWith(sba, spl);
+        }
+    }
+}
+void SBezier::AllIntersectionsWith(SBezier *sbb, SPointList *spl) {
+    SPointList splRaw;
+    ZERO(&splRaw);
+    SEdgeList sea, seb;
+    ZERO(&sea);
+    ZERO(&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();
+}
 
 SBezierLoop SBezierLoop::FromCurves(SBezierList *sbl,
                                     bool *allClosed, SEdge *errorAt)

srf/ratpoly.cpp

@@ -162,6 +162,31 @@ void SBezier::ClosestPointTo(Vector p, double *t, bool converge) {
     }
 }
 
+bool SBezier::PointOnThisAndCurve(SBezier *sbb, Vector *p) {
+    double ta, tb;
+    this->ClosestPointTo(*p, &ta, false);
+    sbb ->ClosestPointTo(*p, &tb, false);
+
+    int i;
+    for(i = 0; i < 20; i++) {   
+        Vector pa = this->PointAt(ta),
+               pb = sbb ->PointAt(tb),
+               da = this->TangentAt(ta),
+               db = sbb ->TangentAt(tb);
+
+        if(pa.Equals(pb, RATPOLY_EPS)) {
+            *p = pa;
+            return true;
+        }
+      
+        double tta, ttb;
+        Vector::ClosestPointBetweenLines(pa, da, pb, db, &tta, &ttb);
+        ta += tta;
+        tb += ttb;
+    }
+    return false;
+}
+
 void SBezier::SplitAt(double t, SBezier *bef, SBezier *aft) {
     Vector4 ct[4];
     int i;
@@ -201,6 +226,16 @@ void SBezier::SplitAt(double t, SBezier *bef, SBezier *aft) {
     }
 }
 
+void SBezier::MakePwlInto(SEdgeList *sel) {
+    List<Vector> lv;
+    ZERO(&lv);
+    MakePwlInto(&lv);
+    int i;
+    for(i = 1; i < lv.n; i++) {
+        sel->AddEdge(lv.elem[i-1], lv.elem[i]);
+    }
+    lv.Clear();
+}
 void SBezier::MakePwlInto(List<SCurvePt> *l) {
     List<Vector> lv;
     ZERO(&lv);

srf/surface.h

@@ -69,14 +69,17 @@ public:
     Vector TangentAt(double t);
     void ClosestPointTo(Vector p, double *t, bool converge=true);
     void SplitAt(double t, SBezier *bef, SBezier *aft);
+    bool PointOnThisAndCurve(SBezier *sbb, Vector *p);
 
     Vector Start(void);
     Vector Finish(void);
     bool Equals(SBezier *b);
+    void MakePwlInto(SEdgeList *sel);
     void MakePwlInto(List<SCurvePt> *l);
     void MakePwlInto(List<Vector> *l);
     void MakePwlWorker(List<Vector> *l, double ta, double tb);
 
+    void AllIntersectionsWith(SBezier *sbb, SPointList *spl);
     void GetBoundingProjd(Vector u, Vector orig, double *umin, double *umax);
     void Reverse(void);
 
@@ -101,6 +104,7 @@ public:
 
     void Clear(void);
     void CullIdenticalBeziers(void);
+    void AllIntersectionsWith(SBezierList *sblb, SPointList *spl);
 };
 
 class SBezierLoop {

ui.h

@@ -347,8 +347,10 @@ public:
 
     void MakeTangentArc(void);
     void SplitLinesOrCurves(void);
-    void SplitLine(hEntity he, Vector pinter);
+    hEntity SplitEntity(hEntity he, Vector pinter);
-    void SplitCircle(hEntity he, Vector pinter);
+    hEntity SplitLine(hEntity he, Vector pinter);
+    hEntity SplitCircle(hEntity he, Vector pinter);
+    hEntity SplitCubic(hEntity he, Vector pinter);
     void ReplacePointInConstraints(hEntity oldpt, hEntity newpt);
     
     // The current selection.

util.cpp

@@ -660,14 +660,15 @@ Vector Vector::AtIntersectionOfPlanes(Vector n1, double d1,
     return (n1.ScaledBy(c1)).Plus(n2.ScaledBy(c2));
 }
 
-Vector Vector::AtIntersectionOfLines(Vector a0, Vector a1,
+void Vector::ClosestPointBetweenLines(Vector a0, Vector da,
-                                     Vector b0, Vector b1,
+                                      Vector b0, Vector db,
-                                     bool *skew,
+                                      double *ta, double *tb)
-                                     double *parama, double *paramb)
 {
-    Vector da = a1.Minus(a0), db = b1.Minus(b0);
+    Vector a1 = a0.Plus(da),
+           b1 = a1.Plus(db);
 
-    // Make an orthogonal coordinate system from those directions
+    // Make a semi-orthogonal coordinate system from those directions;
+    // note that dna and dnb need not be perpendicular.
     Vector dn = da.Cross(db); // normal to both
     Vector dna = dn.Cross(da); // normal to da
     Vector dnb = dn.Cross(db); // normal to db
@@ -676,8 +677,20 @@ Vector Vector::AtIntersectionOfLines(Vector a0, Vector a1,
     //    a0 + pa*da = b0 + pb*db (where pa, pb are scalar params)
     // So dot this equation against dna and dnb to get two equations
     // to solve for da and db
-    double pb =  ((a0.Minus(b0)).Dot(dna))/(db.Dot(dna));
+    *tb =  ((a0.Minus(b0)).Dot(dna))/(db.Dot(dna));
-    double pa = -((a0.Minus(b0)).Dot(dnb))/(da.Dot(dnb));
+    *ta = -((a0.Minus(b0)).Dot(dnb))/(da.Dot(dnb));   
+}
+
+Vector Vector::AtIntersectionOfLines(Vector a0, Vector a1,
+                                     Vector b0, Vector b1,
+                                     bool *skew,
+                                     double *parama, double *paramb)
+{
+    Vector da = a1.Minus(a0), db = b1.Minus(b0);
+
+    double pa, pb;
+    Vector::ClosestPointBetweenLines(a0, da, b0, db, &pa, &pb);
+
     if(parama) *parama = pa;
     if(paramb) *paramb = pb;