change delete points for ink to try to preserve the shape as much as possible. Shift + backspace deletes the point without preserving geometry

4 files changed, 1478 insertions, 13 deletions

diff --git a/src/client/util/bezierFit.ts b/src/client/util/bezierFit.ts
new file mode 100644
index 000000000..57c6dbbde
--- /dev/null
+++ b/src/client/util/bezierFit.ts
@@ -0,0 +1,1446 @@
+import { Point } from "../../pen-gestures/ndollar";
+import { max } from "lodash";
+
+class SmartRect {
+    minx: number = 0;
+    miny: number = 0;
+    maxx: number = 0;
+    maxy: number = 0;
+
+    constructor(mix: number = 0, miy: number = 0, max: number = 0, may: number = 0) { this.minx = mix; this.miny = miy; this.maxx = max; this.maxy = may; }
+
+    public get Center() { return new Point((this.maxx + this.minx) / 2.0, (this.maxy + this.miny) / 2.0); }
+    public get TopLeft() { return new Point(this.minx, this.miny); }
+    public get TopRight() { return new Point(this.maxx, this.miny); }
+    public get BotLeft() { return new Point(this.minx, this.maxy); }
+    public get BotRight() { return new Point(this.maxx, this.maxy); }
+    public get Width() { return this.maxx - this.minx; }
+    public get Height() { return this.maxy - this.miny; }
+    public static Intersect(a: SmartRect, b: SmartRect) { return a.Intersect(b); }
+    public Intersect(b: SmartRect) { return !((this.minx > b.maxx) || (this.miny > b.maxy) || (b.minx > this.maxx) || (b.miny > this.maxy)); }
+
+    public ContainsPercentage(other: SmartRect, axis: Point) {
+        var ret = 0;
+        var minx = Math.max(other.TopLeft.X * axis.X + other.TopLeft.Y * axis.Y, this.TopLeft.X * axis.X + this.TopLeft.Y * axis.Y);
+        var maxx = Math.max(other.BotRight.X * axis.X + other.BotRight.Y * axis.Y, this.BotRight.X * axis.X + this.BotRight.Y * axis.Y);
+        ret = maxx > minx ? (maxx - minx) / (axis == new Point(1, 0) ? other.Width : other.Height) : 0;
+        return ret;
+    }
+    public static Bounds(p: Point[]) {
+        var r = new SmartRect();
+        if (p.length > 0) {
+            r.minx = p[0].X; // These are the most likely to be extremal
+            r.maxx = p.lastElement().X;
+            r.miny = p[0].Y;
+            r.maxy = p.lastElement().Y;
+
+            if (r.minx > r.maxx) {
+                var tmp = r.minx;
+                r.minx = r.maxx;
+                r.maxx = tmp;
+            }
+            if (r.miny > r.maxy) {
+                var tmp = r.miny;
+                r.miny = r.maxy;
+                r.maxy = tmp;
+            }
+
+            for (var pt of p) {
+                if (pt.X < r.minx)
+                    r.minx = pt.X;
+                else if (pt.X > r.maxx)
+                    r.maxx = pt.X;
+
+                if (pt.Y < r.miny)
+                    r.miny = pt.Y;
+                else if (pt.Y > r.maxy)
+                    r.maxy = pt.Y;
+            }
+        }
+        return r;
+    }
+};
+
+function Normalize(p: Point) {
+    const len = Math.sqrt(p.X * p.X + p.Y * p.Y);
+    return new Point(p.X / len, p.Y / len);
+}
+
+function ReparameterizeBezier(d: Point[], first: number, last: number, u: number[], bezCurve: Point[]) {
+    var uPrime = new Array<number>(last - first + 1); //  New parameter values
+
+    for (var i = first; i <= last; i++) {
+        uPrime[i - first] = NewtonRaphsonRootFind(bezCurve, d[i], u[i - first]);
+    }
+    return uPrime;
+}
+function ComputeMaxError(d: Point[], first: number, last: number, bezCurve: Point[], u: number[]) {
+    var maxError = 0; //  Maximum error
+
+    var splitPoint2D = (last - first + 1) / 2;
+    for (var i = first + 1; i < last; i++) {
+        var P = [0, 0];  //  point on curve
+        EvalBezierFast(bezCurve, u[i - first], P);
+        var dx = P[0] - d[i].X;//  offset from point to curve
+        var dy = P[1] - d[i].Y;
+        var dist = Math.sqrt(dx * dx + dy * dy); //  Current error
+        if (dist >= maxError) {
+            maxError = dist;
+            if (splitPoint2D)
+                splitPoint2D = i;
+        }
+    }
+    return { maxError, splitPoint2D };
+}
+function ChordLengthParameterize(d: Point[], first: number, last: number) {
+    var u = new Array<number>(last - first + 1);//  Parameterization
+
+    var prev = 0.0;
+    u[0] = prev;
+    for (var i = first + 1; i <= last; i++) {
+        var lastd = d[i - 1];
+        var curd = d[i];
+        var dx = lastd.X - curd.X;
+        var dy = lastd.Y - curd.Y;
+        prev = u[i - first] = prev + Math.sqrt(dx * dx + dy * dy);
+    }
+
+    var ulastfirst = u[last - first];
+    for (var i = first + 1; i <= last; i++) {
+        u[i - first] /= ulastfirst;
+    }
+
+    return u;
+}
+/*
+*  B0, B1, B2, B3 :
+*  Bezier multipliers
+*/
+function B0(u: number) { var tmp = 1.0 - u; return tmp * tmp * tmp; }
+function B1(u: number) { var tmp = 1.0 - u; return 3 * u * tmp * tmp; }
+function B2(u: number) { var tmp = 1.0 - u; return 3 * u * u * tmp; }
+function B3(u: number) { return u * u * u; }
+function bounds(p: Point[]) {
+    var r = new SmartRect(p[0].X, p[0].Y, p[3].X, p[3].Y); // These are the most likely to be extremal
+
+    if (r.minx > r.maxx) (r.minx, r.maxx);
+    if (r.miny > r.maxy) [r.miny, r.maxy] = [r.maxy, r.miny]; // swap min & max
+
+    for (var i = 1; i < 3; i++) {
+        if (p[i].X < r.minx) r.minx = p[i].X;
+        else if (p[i].X > r.maxx) r.maxx = p[i].X;
+
+        if (p[i].Y < r.miny) r.miny = p[i].Y;
+        else if (p[i].Y > r.maxy) r.maxy = p[i].Y;
+    }
+    return r;
+}
+
+
+
+function splitCubic(p: Point[], t: number, left: Point[], right: Point[]) {
+    var sz = 4;
+    var Vtemp = new Array<Array<Point>>(4);
+    for (var i = 0; i < 4; i++) Vtemp[i] = new Array<Point>(4);
+
+    /* Copy control points */
+    // std::copy(p.begin(), p.end(), Vtemp[0]);
+    for (var i = 0; i < sz; i++) {
+        Vtemp[0][i].X = p[i].X;
+        Vtemp[0][i].Y = p[i].Y;
+    }
+
+    /* Triangle computation */
+    for (var i = 1; i < sz; i++) {
+        for (var j = 0; j < sz - i; j++) {
+            var a = Vtemp[i - 1][j];
+            var b = Vtemp[i - 1][j + 1];
+            Vtemp[i][j].X = b.X * t + a.X * (1 - t);
+            Vtemp[i][j].Y = b.Y * t + a.Y * (1 - t);  // Vtemp[i][j] = Point2D::Lerp(Vtemp[i - 1][j], Vtemp[i - 1][j + 1], t);
+        }
+    }
+
+    if (left) {
+        for (var j = 0; j < sz; j++) {
+            left[j].X = Vtemp[j][0].X;
+            left[j].Y = Vtemp[j][0].Y;
+        }
+    }
+    if (right) {
+        for (var j = 0; j < sz; j++) {
+            right[j].X = Vtemp[sz - 1 - j][j].X;
+            right[j].Y = Vtemp[sz - 1 - j][j].Y;
+        }
+    }
+}
+
+/*
+* Recursively intersect two curves keeping track of their real parameters
+* and depths of intersection.
+* The results are returned in a 2-D array of doubles indicating the parameters
+* for which intersections are found.  The parameters are in the order the
+* intersections were found, which is probably not in sorted order.
+* When an intersection is found, the parameter value for each of the two
+* is stored in the index elements array, and the index is incremented.
+*
+* If either of the curves has subdivisions left before it is straight
+* (depth > 0)
+* that curve (possibly both) is (are) subdivided at its (their) midpoint(s).
+* the depth(s) is (are) decremented, and the parameter value(s) corresponding
+* to the midpoints(s) is (are) computed.
+* Then each of the subcurves of one curve is intersected with each of the
+* subcurves of the other curve, first by testing the bounding boxes for
+* interference.  If there is any bounding box interference, the corresponding
+* subcurves are recursively intersected.
+*
+* If neither curve has subdivisions left, the line segments from the first
+* to last control point of each segment are intersected.  (Actually the
+* only the parameter value corresponding to the intersection point is found).
+*
+* The apriori flatness test is probably more efficient than testing at each
+* level of recursion, although a test after three or four levels would
+* probably be worthwhile, since many curves become flat faster than their
+* asymptotic rate for the first few levels of recursion.
+*
+* The bounding box test fails much more frequently than it succeeds, providing
+* substantial pruning of the search space.
+*
+* Each (sub)curve is subdivided only once, hence it is not possible that for
+* one final line intersection test the subdivision was at one level, while
+* for another final line intersection test the subdivision (of the same curve)
+* was at another.  Since the line segments share endpoints, the intersection
+* is robust: a near-tangential intersection will yield zero or two
+* intersections.
+*/
+function recursively_intersect(a: Point[], t0: number, t1: number, deptha: number, b: Point[], u0: number, u1: number, depthb: number, parameters: number[][]) {
+    if (deptha > 0) {
+        var a1 = new Array<Point>(4), a2 = new Array<Point>(4);
+        splitCubic(a, 0.5, a1, a2);
+        var tmid = (t0 + t1) * 0.5;
+        deptha--;
+        if (depthb > 0) {
+            var b1 = new Array<Point>(4), b2 = new Array<Point>(4);
+            splitCubic(b, 0.5, b1, b2);
+            var umid = (u0 + u1) * 0.5;
+            depthb--;
+            if (SmartRect.Intersect(bounds(a1), bounds(b1))) {
+                recursively_intersect(a1, t0, tmid, deptha, b1, u0, umid, depthb, parameters);
+            }
+            if (SmartRect.Intersect(bounds(a2), bounds(b1))) {
+                recursively_intersect(a2, tmid, t1, deptha, b1, u0, umid, depthb, parameters);
+            }
+            if (SmartRect.Intersect(bounds(a1), bounds(b2))) {
+                recursively_intersect(a1, t0, tmid, deptha, b2, umid, u1, depthb, parameters);
+            }
+            if (SmartRect.Intersect(bounds(a2), bounds(b2))) {
+                recursively_intersect(a2, tmid, t1, deptha, b2, umid, u1, depthb, parameters);
+            }
+        }
+        else {
+            if (SmartRect.Intersect(bounds(a1), bounds(b))) {
+                recursively_intersect(a1, t0, tmid, deptha, b, u0, u1, depthb, parameters);
+            }
+            if (SmartRect.Intersect(bounds(a2), bounds(b))) {
+                recursively_intersect(a2, tmid, t1, deptha, b, u0, u1, depthb, parameters);
+            }
+        }
+    }
+    else {
+        if (depthb > 0) {
+            var b1 = new Array<Point>(4), b2 = new Array<Point>(4);
+            splitCubic(b, 0.5, b1, b2);
+            var umid = (u0 + u1) * 0.5;
+            depthb--;
+            if (SmartRect.Intersect(bounds(a), bounds(b1))) {
+                recursively_intersect(a, t0, t1, deptha, b1, u0, umid, depthb, parameters);
+            }
+            if (SmartRect.Intersect(bounds(a), bounds(b2))) {
+                recursively_intersect(a, t0, t1, deptha, b2, umid, u1, depthb, parameters);
+            }
+        }
+        else // Both segments are fully subdivided; now do line segments
+        {
+            var xlk = a[3].X - a[0].X;
+            var ylk = a[3].Y - a[0].Y;
+            var xnm = b[3].X - b[0].X;
+            var ynm = b[3].Y - b[0].Y;
+            var xmk = b[0].X - a[0].X;
+            var ymk = b[0].Y - a[0].Y;
+            var det = xnm * ylk - ynm * xlk;
+            if (1.0 + det == 1.0) {
+                return;
+            }
+            else {
+                var detinv = 1.0 / det;
+                var s = (xnm * ymk - ynm * xmk) * detinv;
+                var t = (xlk * ymk - ylk * xmk) * detinv;
+                if ((s < 0.0) || (s > 1.0) || (t < 0.0) || (t > 1.0) || Number.isNaN(s) || Number.isNaN(t)) {
+                    return;
+                }
+                parameters.push([t0 + s * (t1 - t0), u0 + t * (u1 - u0)]);
+            }
+        }
+    }
+}
+
+/*
+*  EvalBezier :
+*      Evaluate a Bezier curve at a particular parameter value
+*
+*/
+const MAX_DEGREE = 5;
+function EvalBezier(V: Point[], degree: number, t: number, result: number[]) {
+    if (degree + 1 > MAX_DEGREE) {
+        result[0] = V[0].X;
+        result[1] = V[0].Y;
+        return;
+    }
+
+    var Vtemp = [new Point(0, 0), new Point(0, 0), new Point(0, 0), new Point(0, 0)]; // Local copy of control points
+
+    /* Copy array   */
+    for (var i = 0; i <= degree; i++) {
+        Vtemp[i].X = V[i].X;
+        Vtemp[i].Y = V[i].Y;
+    }
+
+    /* Triangle computation */
+    for (var i = 1; i <= degree; i++) {
+        for (var j = 0; j <= degree - i; j++) {
+            Vtemp[j].X = (1.0 - t) * Vtemp[j].X + t * Vtemp[j + 1].X;
+            Vtemp[j].Y = (1.0 - t) * Vtemp[j].Y + t * Vtemp[j + 1].Y;
+        }
+    }
+
+    result[0] = Vtemp[0].X;
+    result[1] = Vtemp[0].Y;// Point on curve at parameter t
+}
+
+function EvalBezierFast(p: Point[], t: number, result: number[]) {
+    var n = 3;
+    var u: number, bc: number, tn: number, tmpX: number, tmpY: number;
+    u = 1.0 - t;
+    bc = 1;
+    tn = 1;
+
+    tmpX = p[0].X * u;
+    tmpY = p[0].Y * u;
+    tn = tn * t;
+    bc = bc * (n - 1 + 1) / 1;
+    tmpX = (tmpX + tn * bc * p[1].X) * u;
+    tmpY = (tmpY + tn * bc * p[1].Y) * u;
+    tn = tn * t;
+    bc = bc * (n - 2 + 1) / 2;
+    tmpX = (tmpX + tn * bc * p[2].X) * u;
+    tmpY = (tmpY + tn * bc * p[2].Y) * u;
+
+    result[0] = tmpX + tn * t * p[3].X;
+    result[1] = tmpY + tn * t * p[3].Y;
+}
+/*
+* ComputeLeftTangent, ComputeRightTangent, ComputeCenterTangent :
+*Approximate unit tangents at endpoints and "center" of digitized curve
+*/
+function ComputeLeftTangent(d: Point[], end: number) {
+    var use = 1;
+    var tHat1 = new Point(d[end + use].X - d[end].X, d[end + use].Y - d[end].Y);
+    return Normalize(tHat1);
+}
+function ComputeRightTangent(d: Point[], end: number) {
+    var available = end;
+    var use = 1;
+    var tHat2 = new Point(d[end - use].X - d[end].X, d[end - use].Y - d[end].Y);
+    return Normalize(tHat2);
+}
+function ComputeCenterTangent(d: Point[], center: number) {
+    if (center == 0) {
+        return ComputeLeftTangent(d, center);
+    }
+    var V1 = ComputeLeftTangent(d, center); // d[center] - d[center-1];
+    var V2 = ComputeRightTangent(d, center); // d[center] - d[center + 1];
+    var tHatCenter = new Point((-V1.X + V2.X) / 2.0, (-V1.Y + V2.Y) / 2.0);
+    if (tHatCenter === new Point(0, 0)) {
+        tHatCenter = new Point(-V1.Y, -V1.X);// V1.Perp();
+    }
+    return Normalize(tHatCenter);
+}
+function GenerateBezier(d: Point[], first: number, last: number, uPrime: number[], tHat1: Point, tHat2: Point, result: Point[] /* must be prealloacted to size 4 */) {
+    var nPts = last - first + 1; // Number of pts in sub-curve
+    var Ax = new Array<number>(nPts * 2);// Precomputed rhs for eqn  //std::vector<Vector2D> A(nPts * 2);
+    var Ay = new Array<number>(nPts * 2);// Precomputed rhs for eqn  //std::vector<Vector2D> A(nPts * 2);
+
+    /* Compute the A's  */
+    for (var i = 0; i < nPts; i++) {
+        var uprime = uPrime[i];
+        var b1 = B1(uprime);
+        var b2 = B2(uprime);
+        Ax[i] = tHat1.X * b1;
+        Ay[i] = tHat1.Y * b1;
+        Ax[i + 1 * nPts] = tHat2.X * b2;
+        Ay[i + 1 * nPts] = tHat2.Y * b2;
+    }
+
+    /* Create the C and X matrices  */
+    var C = [[0, 0], [0, 0]];
+    var df = d[first];
+    var dl = d[last];
+
+    var X = [0, 0];    // Matrix X
+    for (var i = 0; i < nPts; i++) {
+        C[0][0] += Ax[i] * Ax[i] + Ay[i] * Ay[i];     //A[i+0*nPts].Dot(A[i+0*nPts]);
+        C[0][1] += Ax[i] * Ax[i + nPts] + Ay[i] * Ay[i + nPts];//A[i+0*nPts].Dot(A[i+1*nPts]);
+        C[1][0] = C[0][1];
+        C[1][1] += Ax[i + nPts] * Ax[i + nPts] + Ay[i + nPts] * Ay[i + nPts];// A[i+1*nPts].Dot(A[i+1*nPts]);
+        var uprime = uPrime[i];
+        var b0plb1 = B0(uprime) + B1(uprime);
+        var b2plb3 = B2(uprime) + B3(uprime);
+        var df1 = d[first + i];
+        var tmpX = df1.X - (df.X * b0plb1 + (dl.X * b2plb3));
+        var tmpY = df1.Y - (df.Y * b0plb1 + (dl.Y * b2plb3));
+
+        X[0] += Ax[i] * tmpX + Ay[i] * tmpY; // A[i+0*nPts].Dot(tmp)
+        X[1] += Ax[i + nPts] * tmpX + Ay[i + nPts] * tmpY; //A[i+1*nPts].Dot(tmp)
+    }
+
+    /* Compute the determinants of C and X  */
+    var det_C0_C1 = C[0][0] * C[1][1] - C[1][0] * C[0][1];
+    var det_C0_X = C[0][0] * X[1] - C[0][1] * X[0];
+    var det_X_C1 = X[0] * C[1][1] - X[1] * C[0][1];
+
+    /* Finally, derive alpha values */
+    if (det_C0_C1 == 0.0) {
+        det_C0_C1 = (C[0][0] * C[1][1]) * 10e-12;
+    }
+    var alpha_l = (det_C0_C1 == 0) ? 0.0 : det_X_C1 / det_C0_C1;
+    var alpha_r = (det_C0_C1 == 0) ? 0.0 : det_C0_X / det_C0_C1;
+
+    /* If alpha negative, use the Wu/Barsky heuristic (see text) */
+    /* (if alpha is 0, you get coincident control points that lead to
+    * divide by zero in any subsequent NewtonRaphsonRootFind() call. */
+    var segLength = Math.sqrt((df.X - dl.X) * (df.X - dl.X) + (df.Y - dl.Y) * (df.Y - dl.Y));
+    var epsilon = 1.0e-6 * segLength;
+    if (alpha_l < epsilon || alpha_r < epsilon) {
+        /* fall back on standard (probably inaccurate) formula, and subdivide further if needed. */
+        alpha_l = alpha_r = segLength / 3.0;
+    }
+
+    /*  First and last control points of the Bezier curve are */
+    /*  positioned exactly at the first and last data points */
+    /*  Control points 1 and 2 are positioned an alpha distance out */
+    /*  on the tangent vectors, left and right, respectively */
+    result[0] = df;// RETURN bezier curve ctl pts
+    result[3] = dl;
+    result[1] = new Point(df.X + (tHat1.X * alpha_l), df.Y + (tHat1.Y * alpha_l));
+    result[2] = new Point(dl.X + (tHat2.X * alpha_r), dl.Y + (tHat2.Y * alpha_r));
+
+}
+/*
+ *  NewtonRaphsonRootFind :
+ *  Use Newton-Raphson iteration to find better root.
+ */
+function NewtonRaphsonRootFind(Q: Point[], P: Point, u: number) {
+    var Q1 = [new Point(0, 0), new Point(0, 0), new Point(0, 0)], Q2 = [new Point(0, 0), new Point(0, 0)]; //  Q' and Q''
+    var Q_u = [0, 0], Q1_u = [0, 0], Q2_u = [0, 0]; //u evaluated at Q, Q', & Q''
+
+    /* Compute Q(u) */
+    var uPrime: number; //  Improved u
+    EvalBezierFast(Q, u, Q_u);
+
+    /* Generate control vertices for Q' */
+    for (var i = 0; i <= 2; i++) {
+        Q1[i].X = (Q[i + 1].X - Q[i].X) * 3.0;
+        Q1[i].Y = (Q[i + 1].Y - Q[i].Y) * 3.0;
+    }
+
+    /* Generate control vertices for Q'' */
+    for (var i = 0; i <= 1; i++) {
+        Q2[i].X = (Q1[i + 1].X - Q1[i].X) * 2.0;
+        Q2[i].Y = (Q1[i + 1].Y - Q1[i].Y) * 2.0;
+    }
+
+    /* Compute Q'(u) and Q''(u) */
+    EvalBezier(Q1, 2, u, Q1_u);
+    EvalBezier(Q2, 1, u, Q2_u);
+
+    /* Compute f(u)/f'(u) */
+    var numerator = (Q_u[0] - P.X) * (Q1_u[0]) + (Q_u[1] - P.Y) * (Q1_u[1]);
+    var denominator = (Q1_u[0]) * (Q1_u[0]) + (Q1_u[1]) * (Q1_u[1]) + (Q_u[0] - P.X) * (Q2_u[0]) + (Q_u[1] - P.Y) * (Q2_u[1]);
+    if (denominator == 0.0)
+        uPrime = u;
+    else uPrime = u - (numerator / denominator);/* u = u - f(u)/f'(u) */
+
+    return uPrime;
+}
+function FitCubic(d: Point[], first: number, last: number, tHat1: Point, tHat2: Point, error: number, result: Point[]) {
+    var bezCurve = new Array<Point>(4);          //  Control points of fitted Bezier curve
+    var splitPoint2D: number;           //  Point2D to split point set at
+    var maxIterations = 4;      //  Max times to try iterating
+
+    var iterationError = error * error; //  Error below which you try iterating
+    var nPts = last - first + 1;        //  Number of points in subset
+
+    /*  Use heuristic if region only has two points in it */
+    if (nPts == 2) {
+        var dist = Math.sqrt((d[first].X - d[last].X) * (d[first].X - d[last].X) + (d[first].Y - d[last].Y) * (d[first].Y - d[last].Y)) / 3;
+
+        bezCurve[0] = d[first];
+        bezCurve[3] = d[last];
+        bezCurve[1] = new Point(bezCurve[0].X + (tHat1.X * dist), bezCurve[0].Y + (tHat1.Y * dist));
+        bezCurve[2] = new Point(bezCurve[3].X + (tHat2.X * dist), bezCurve[3].Y + (tHat2.Y * dist));
+
+        result.push(bezCurve[1]);
+        result.push(bezCurve[2]);
+        result.push(bezCurve[3]);
+        return;
+    }
+
+    /*  Parameterize points, and attempt to fit curve */
+    var u = ChordLengthParameterize(d, first, last);
+    GenerateBezier(d, first, last, u, tHat1, tHat2, bezCurve);
+
+    /*  Find max deviation of points to fitted curve */
+    var { maxError, splitPoint2D } = ComputeMaxError(d, first, last, bezCurve, u);  //  Maximum fitting error
+    if (maxError < Math.abs(error)) {
+        result.push(bezCurve[1]);
+        result.push(bezCurve[2]);
+        result.push(bezCurve[3]);
+        return;
+    }
+
+    /*  If error not too large, try some reparameterization  */
+    /*  and iteration */
+    if (maxError < iterationError) {
+        for (var i = 0; i < maxIterations; i++) {
+            var uPrime = ReparameterizeBezier(d, first, last, u, bezCurve);   //  Improved parameter values
+            GenerateBezier(d, first, last, uPrime, tHat1, tHat2, bezCurve);
+            var { maxError, splitPoint2D } = ComputeMaxError(d, first, last, bezCurve, uPrime);
+            if (maxError < error) {
+                result.push(bezCurve[1]);
+                result.push(bezCurve[2]);
+                result.push(bezCurve[3]);
+                return;
+            }
+            u = uPrime;
+        }
+    }
+
+    /* Fitting failed -- split at max error point and fit recursively */
+    var tHatCenter = splitPoint2D >= last - 1 ? ComputeRightTangent(d, splitPoint2D) : ComputeCenterTangent(d, splitPoint2D);
+    FitCubic(d, first, splitPoint2D, tHat1, tHatCenter, error, result);
+    var negThatCenter = new Point(-tHatCenter.X, -tHatCenter.Y);
+    FitCubic(d, splitPoint2D, last, negThatCenter, tHat2, error, result);
+}
+export function FitCurve(d: Point[], error: number) {
+    var tHat1 = ComputeLeftTangent(d, 0); //  Unit tangent vectors at endpoints
+    var tHat2 = ComputeRightTangent(d, d.length - 1);
+    var result: Point[] = [];
+    result.push(d[0]);
+    FitCubic(d, 0, d.length - 1, tHat1, tHat2, error, result);
+    return result;
+}
+export function FitOneCurve(d: Point[], tHat1?: Point, tHat2?: Point) {
+    tHat1 = tHat1 ?? Normalize(ComputeLeftTangent(d, 0));
+    tHat2 = tHat2 ?? Normalize(ComputeRightTangent(d, d.length - 1));
+    tHat2 = new Point(-tHat2.X, -tHat2.Y);
+    var u = ChordLengthParameterize(d, 0, d.length - 1);
+    var bezCurveCtrls = [new Point(0, 0), new Point(0, 0), new Point(0, 0), new Point(0, 0)];
+    GenerateBezier(d, 0, d.length - 1, u, tHat1, tHat2, bezCurveCtrls);    /*  Find max deviation of points to fitted curve */
+    var finalCtrls = bezCurveCtrls.slice();
+    var { maxError: error } = ComputeMaxError(d, 0, d.length - 1, bezCurveCtrls, u);
+    for (var i = 0; i < 10; i++) {
+        var uPrime = ReparameterizeBezier(d, 0, d.length - 1, u, bezCurveCtrls);   //  Improved parameter values
+        GenerateBezier(d, 0, d.length - 1, uPrime, tHat1, tHat2, bezCurveCtrls);
+        var { maxError } = ComputeMaxError(d, 0, d.length - 1, bezCurveCtrls, uPrime);
+        if (maxError < error) {
+            error = maxError;
+            finalCtrls = bezCurveCtrls.slice();
+        }
+        u = uPrime;
+    }
+    return { finalCtrls, error };
+}
+
+/*
+static double                GetTValueFromSValue  (const BezierRep &parent, double t, double endT, bool left, double influenceDistance, double &excess) {
+double dist = 0;
+double step = 0.01;
+double spliceT = t;
+for (spliceT = t+(left?-1:1)*step; dist < influenceDistance && (left ? (spliceT > endT) : (spliceT < endT)); spliceT += step * (left ? -1 : 1)) {
+dist +=  (parent[spliceT]-parent[spliceT-step*(left ? -1:1)]).Length();
+}
+if ((left && spliceT < endT) || (!left && spliceT > endT))
+spliceT = endT;
+excess = influenceDistance - dist;
+return spliceT;
+}
+static BezierRep::BezierLock FindSplitIndex       (const BezierRep &parent, double t, bool left, std::vector<BezierRep::BezierLock> &locked)
+{
+BezierRep::BezierLock  cuspIndex = { left ? 0.0 : 1.0*parent.MaxIndex(), true};
+double                 tprev = t;
+for (int tstep = (left ? std::floor(t) : std::ceil(t)) * 3; left ? (tstep >= 0) : (tstep < parent.p.size()); tstep += (left ? -1 : 1))
+{
+double near = HUGE_VAL;
+for (auto &l : locked) {
+if ((( left && tprev > l.T && tstep <= l.T) ||
+(!left && tprev < l.T && tstep >= l.T)) && std::abs(tprev-l.T) < near) {
+near = std::abs(tprev-l.T);
+cuspIndex = l;
+}
+}
+if (near != HUGE_VAL)
+break;
+}
+return cuspIndex;
+}
+size_t                       SampleBezier         (const BezierRep &bez, Point2D *&multiSegmentSamplePts, size_t numMultiSegmentSamples, size_t samplesPerSegment)
+{
+auto numSamples = bez.MaxIndex() * samplesPerSegment + 1;
+if (numSamples > numMultiSegmentSamples) {
+if (numMultiSegmentSamples)
+delete [] multiSegmentSamplePts;
+multiSegmentSamplePts  = new Point2D[numSamples];
+}
+for (auto seg = 0; seg < bez.MaxIndex(); seg++)
+{
+   double  result[2];
+Point2D tmp[4] = { bez.p[seg * 3], bez.p[seg * 3 +1], bez.p[seg * 3 +2], bez.p[seg * 3 +3] };
+for (auto index = 0; index < samplesPerSegment; index++) {
+EvalBezierFast(tmp, 1.0 * index / samplesPerSegment, result);
+multiSegmentSamplePts[seg * samplesPerSegment + index].X = result[0];
+multiSegmentSamplePts[seg * samplesPerSegment + index].Y = result[1];
+}
+}
+multiSegmentSamplePts[numSamples-1] = bez.p.back();
+return numSamples;
+}
+static double                GetSpliceCurve       (const BezierRep &parent, double t, BezierRep::BezierLock * isCusp, BezierRep::BezierLock tEnd, std::vector<BezierRep::BezierLock> &locked, const Vector2D &v, Point2D singleSegmentSpliceCurve[4], double errorTolerance, double influenceDistance, double &excess)
+{
+Point2D  *multiSegmentSamplePts  = NULL;
+size_t    numMultiSegmentSamples = 0;
+double    spliceT    = tEnd.T;
+bool      left       = tEnd.T < t;
+auto      parTangent = parent.Tangent(t + (left ? -1e-7:1e-7));
+if (_isnan(parTangent.X))
+parTangent = Vector2D();
+for  (auto &l : locked) {
+if (l.T == t && l.Cusp) {
+parTangent = Vector2D();
+if (left && (l.Side == 2) && t<= parent.MaxIndex())
+parTangent = (parent[t+1]-(parent[t]+v)).Normal();
+else if (!left && (l.Side == 1) && t >= 1)
+parTangent = (parent[t]+v - parent[t-1]).Normal();
+}
+}
+
+if (_isnan(influenceDistance) && isCusp && abs(tEnd.T - t) <= 1 && tEnd.Cusp && (((tEnd.Side & 2)  && left) || ((tEnd.Side & 1) && !left))) {
+singleSegmentSpliceCurve[0] = parent[ left ? tEnd.T : t];
+singleSegmentSpliceCurve[2] = parent[!left ? tEnd.T : t];
+singleSegmentSpliceCurve[1] = singleSegmentSpliceCurve[0];
+singleSegmentSpliceCurve[3] = singleSegmentSpliceCurve[2];
+return spliceT;
+}
+
+for (auto startSample = t, endSample = tEnd.T; !((left && startSample < endSample+1e-5) || (!left && startSample > endSample-1e-5)); spliceT = (endSample + startSample)/2)
+{
+if (!_isnan(influenceDistance)) // if influenceDistance has been set, we just use it without subdividing.
+endSample = startSample = spliceT = GetTValueFromSValue(parent, t, tEnd.T, left, influenceDistance, excess);
+
+bool endCusp = spliceT == tEnd.T && tEnd.Cusp && tEnd.Side == 3;
+auto multiSegmentSplitCurve = BezierRep(parent.Split(left ? spliceT : t, left ? t : spliceT) );
+double singleToMultiSegmentError = 0;
+if (multiSegmentSplitCurve.p.size() == 4) { // if split curve is a single-segment bezier, then we it should be 100% accurate
+singleSegmentSpliceCurve[0] = multiSegmentSplitCurve.p[0];
+singleSegmentSpliceCurve[3] = multiSegmentSplitCurve.p[3];
+singleSegmentSpliceCurve[1] = !left && isCusp && isCusp->Side == 3 ? singleSegmentSpliceCurve[0] : multiSegmentSplitCurve.p[1];
+singleSegmentSpliceCurve[2] =  left && isCusp && isCusp->Side == 3 ? singleSegmentSpliceCurve[3] : multiSegmentSplitCurve.p[2];
+if (spliceT == endSample)
+break;
+} else {
+const size_t SAMPLES_PER_SEGMENT = 20;
+numMultiSegmentSamples = SampleBezier(multiSegmentSplitCurve, multiSegmentSamplePts, numMultiSegmentSamples, SAMPLES_PER_SEGMENT);
+
+auto   endTan = (endSample == tEnd.T && tEnd.Cusp && tEnd.Side == (left ? 1 : 2)) ? parent.Tangent(endSample + (left ? -0.001 : 0.001)) : multiSegmentSplitCurve.Tangent(left ? 0.0 : 1.0*multiSegmentSplitCurve.MaxIndex());
+auto   tHat1  = endCusp &&  left ? Vector2D() : !left ?  parTangent :  endTan;
+auto   tHat2  = endCusp && !left ? Vector2D() :  left ? -parTangent : -endTan;
+auto   u      = BezierRep::ChordLengthParameterize(multiSegmentSamplePts, 0, numMultiSegmentSamples-1);
+GenerateBezier(multiSegmentSamplePts, 0, numMultiSegmentSamples-1, u, tHat1, tHat2, singleSegmentSpliceCurve);
+
+singleToMultiSegmentError = BezierRep::ComputeMaxError(multiSegmentSamplePts, 0, numMultiSegmentSamples-1, singleSegmentSpliceCurve, u);
+}
+if (singleToMultiSegmentError > (endCusp ? 5 : 1) * errorTolerance)
+endSample   = spliceT;
+else startSample = spliceT;
+}
+
+if (numMultiSegmentSamples)
+delete [] multiSegmentSamplePts;
+return spliceT;
+}
+
+static void MoveCurveSplice(double t, Point2D splice[4], BezierRep::BezierLock &stepLock, double &extra, bool left, BezierRep::BezierLock *moveLock, const Vector2D &v, double influenceDistance, const Vector2D &smoothParTangent, double ctrlPtScale, double ctrlPtRotate)
+{
+if (!_isnan(influenceDistance) && influenceDistance < 0) {
+splice[left ? 2 : 1] = (splice[left ? 3 : 0] += v);
+if (moveLock) {
+moveLock->Side |= (left ? 1 : 2);
+moveLock->Cusp = true;
+}
+}
+else {
+auto tan = (splice[left?2:1]-splice[left?3:0]);
+splice[left?3:0] += v;
+splice[left?2:1]  = splice[left?3:0] + Mat::Rotate(ctrlPtRotate) * tan * ctrlPtScale;
+if (influenceDistance > 0 && t <= stepLock.T ) {
+LnSeg   tangent(splice[left?3:0], tan == Vector2D() ? (splice[left?3:0]+smoothParTangent):splice[left?2:1]), otherTangent(splice[left?0:3], splice[left?1:2]);
+auto    inter = otherTangent.LnIntersection(tangent);
+auto    seglen = (splice[0] - splice[3]).Length();
+if (inter == Point2D::Null()) {
+if (otherTangent.Length() == 0) {
+auto ang = tangent.Direction().UnsignedAngle(splice[left?0:3]-splice[left?3:0]) / M_PI;
+auto target = splice[left?3:0] +  tangent.Direction() * .5519 * (ang < 0.01 ? 0 : 1) * seglen;
+
+splice[left?2:1] = (target * std::min(1.0, extra/25) + splice[left?2:1] * std::max(0.0, 1-extra/25));
+}
+} else {
+bool behind = tangent.ClosestFraction(inter) <= 0 && otherTangent.ClosestFraction(inter) <= 0;
+auto tandir   = tangent.ClosestFraction(inter) <=0 ? -tangent.Direction() : tangent.Direction();
+auto leglen   = std::max(seglen/4, (splice[left?3:0]-inter).Length());
+auto aspect   = std::sqrt(leglen / seglen / .7071);
+auto modinter =  splice[left?3:0] + tandir * leglen*std::min(1.0,.5519/aspect);
+if (leglen / seglen > 2) {
+if (tangent.Direction().Dot(otherTangent.Direction()) < 0)
+modinter = splice[left?3:0] + tandir * seglen*(.5519);
+else modinter = splice[left?3:0] + tandir * seglen*.7071;
+}
+if (behind && tangent.Direction().Dot(otherTangent.Direction()) < 0)
+modinter = (splice[0] + splice[3])/2;
+auto targetFrac  = (modinter-splice[left?3:0]).Length();
+auto target      = splice[left?3:0] + targetFrac * tangent.Direction();
+splice[left?2:1] = (target * std::min(1.0, extra/25) + splice[left?2:1] * std::max(0.0, 1-extra/25));
+if (extra> 25) {
+//LnSeg   tangent(splice[left?3:0], splice[left?2:1]), otherTangent(splice[left?0:3], splice[left?1:2]);
+auto oextra = extra - 25;
+auto otandir   = otherTangent.ClosestFraction(inter) <=0 ? -otherTangent.Direction() : otherTangent.Direction();
+auto oleglen   = std::max(seglen/4, (splice[!left?3:0]-inter).Length());
+auto oaspect   = std::sqrt(oleglen / seglen / .7071);
+auto omodinter =  splice[!left?3:0] + otandir * oleglen*std::min(1.0,.5519/oaspect);
+if (oleglen/ seglen > 2) {
+if (tangent.Direction().Dot(otherTangent.Direction()) < 0)
+omodinter = splice[!left?3:0] + tandir * seglen*(.5519);
+else omodinter = splice[!left?3:0] + otandir * seglen *.7071;
+}
+if (behind && tangent.Direction().Dot(otherTangent.Direction()) < 0)
+omodinter = (splice[0] + splice[3])/2;
+auto otargetFrac  = (omodinter-splice[!left?3:0]).Length();
+auto otarget      = splice[!left?3:0] + otargetFrac * otherTangent.Direction();
+splice[!left?2:1] = (otarget * std::min(1.0, oextra/25) + splice[!left?2:1] * std::max(0.0, 1-oextra/25));
+}
+}
+}
+}
+}
+static void                  MoveTAux             (BezierRep &curve, double tMove, const Vector2D &v, bool moveEnds)
+{
+auto &p = curve.p;
+auto tstart = static_cast<int>(tMove);
+auto tend   = static_cast<int>(ceil(tMove));
+if (tend == tstart)
+{
+if (tend == 0)
+{
+tend = 1;
+}
+else
+{
+tstart = tend - 1;
+}
+}
+auto t = tMove - tstart;
+
+auto b0 = pow(1 - t, 3);
+auto b1 = 3 * t * pow(1 - t, 2);
+auto b2 = 3 * t * t * (1 - t);
+auto b3 = t * t * t;
+
+auto ind = t < 0.4 ? 1 : t > 0.6 ? -1 : 0;
+if (ind == 0) {
+auto norm = (b1 + b2);
+p[tstart * 3 + 1] +=  (b1/norm * v)/b1;
+p[tend * 3   - 1] +=  (b2/norm * v)/b2;
+}
+else if (ind ==  1 && b1 != 0) {
+auto pt = curve[tMove] + v;
+p[tstart * 3 +1].X = (pt.X - b0 * p[tstart*3].X - b3 * p[tend*3].X - b2 * p[tend*3-1].X) / b1;
+p[tstart * 3 +1].Y = (pt.Y - b0 * p[tstart*3].Y - b3 * p[tend*3].Y - b2 * p[tend*3-1].Y) / b1;
+}
+else if (ind == -1 && b2 != 0) {
+auto pt = curve[tMove] + v;
+p[tend * 3 -1].X = (pt.X - b0 * p[tstart*3].X - b3 * p[tend*3].X - b1 * p[tstart*3+1].X) / b2;
+p[tend * 3 -1].Y = (pt.Y - b0 * p[tstart*3].Y - b3 * p[tend*3].Y - b1 * p[tstart*3+1].Y) / b2;
+}
+
+if (moveEnds) {
+p[tstart * 3] += b0 * v;
+p[tend   * 3] += b3 * v;
+}
+
+//p[tstart * 3 + 1] += b1 * v;
+//p[tend * 3 - 1] += b2 * v;
+//p[tstart*3] += v * b0;
+//p[tend*3] += v * b3;
+
+//            fx(t):=(1−t)3p1x+3t(1−t)2p2x+3t2(1−t)p3x+t3p4x
+//fy(t):=(1−t)3p1y+3t(1−t)2p2y+3t2(1−t)p3y+t3p4y
+
+//Call the curve C(t) = b0(t) P0 + b1(t) P1 + b2(t) P2 + b3(t) P3. The user clicks at some point Q and drags to a new point R.
+//            3. Compute c0 = b0(s); c1 = b1(s), c2 = b2(s), and c3 = b3(s), the coefficients of the control points at parameter s.
+
+//4. Adjust the Ps like this:
+
+//P0 += c0 * v
+//P1 += c1 * v;
+//P2 += c2 * v;
+//P3 += c3 * v.
+}
+static void                  MoveTAdaptive        (BezierRep &curve, double tMove, const Vector2D &v, std::vector<BezierRep::BezierLock> &locked, bool rangeDrag, double errorTolerance, double influenceLDistance, double influenceRDistance, double ctrlPtLRotate, double ctrlPtRRotate, double ctrlPtScale, bool moveEnds)
+{
+auto tleftMove = rangeDrag ? (static_cast<int>(tMove) == tMove ? tMove-1 : static_cast<int>(tMove)) : tMove;
+auto trightMove = rangeDrag ? (static_cast<int>(tMove) == tMove ? tMove+1 : static_cast<int>(tMove)+1) : tMove;
+auto leftStep  = FindSplitIndex(curve,tleftMove, true,  locked);
+auto rightStep = FindSplitIndex(curve, trightMove, false, locked);
+auto leftTan   = curve.Tangent(std::max(0.0, tleftMove-1e-5));
+auto rightTan  = curve.Tangent(std::min(curve.MaxIndex() * 1.0, trightMove + 1e-5));
+auto smoothParTangent = (leftTan + rightTan)/2;
+auto smoothParDist    = (curve[std::max(0.0, tMove-1)] - curve[std::min(curve.MaxIndex() * 1.0, tMove+1)]).Length()/4;
+
+BezierRep::BezierLock *isCusp = NULL, *moveLock = NULL;
+for (auto &lck : locked) {
+if (lck.T == tMove) {
+moveLock = &lck;
+if (influenceLDistance > 0 && moveLock && moveLock->Cusp) {
+moveLock->Cusp = false;
+if (moveLock->T*3 - 1 >= 0)
+curve.p[moveLock->T * 3 - 1] = curve.p[moveLock->T * 3] - smoothParTangent*smoothParDist;
+if (moveLock->T*3 + 1 < curve.p.size())
+curve.p[moveLock->T * 3 + 1] = curve.p[moveLock->T * 3] + smoothParTangent*smoothParDist;
+}
+if (lck.Cusp)
+isCusp = &lck;
+break;
+}
+}
+if (moveLock && moveLock->T * 3 -1 >= 0 && moveLock->T*3+1 < curve.p.size() &&
+curve.p[moveLock->T*3-1] == curve.p[moveLock->T*3+1])
+leftTan = rightTan = smoothParTangent;
+// splice the left side of the point that is moved
+Point2D spliceL[4], spliceR[4];
+double lextra=0, rextra= 0;
+auto l = GetSpliceCurve(curve, tleftMove,  isCusp, leftStep,  locked, v, spliceL, errorTolerance, abs(influenceLDistance), lextra);
+auto r = GetSpliceCurve(curve, trightMove, isCusp, rightStep, locked, v, spliceR, errorTolerance, abs(influenceRDistance), rextra);
+
+BezierRep   splicedCurve;
+if (tMove != 0) {
+if (l == -1)
+return;
+
+MoveCurveSplice(l, spliceL, leftStep, lextra, true, moveLock, v, influenceLDistance, -rightTan, ctrlPtScale, ctrlPtLRotate);
+
+// add on the remaining left side of the curve
+if (l != 0)
+splicedCurve = curve.Split(0,l);
+
+// add the spliced left side of the curve
+for (auto i = l == 0 ? 0 : 1; i < 4; i++)
+splicedCurve.p.push_back(spliceL[i]);
+
+if (tleftMove != tMove)
+{
+auto fixedL = curve.Split(tleftMove, tMove);
+for (auto i = 1; i < 4; i++)
+splicedCurve.p.push_back(fixedL[i] + v);
+}
+}
+
+auto moveIndex = splicedCurve.p.size();
+auto insertEnd = moveIndex;
+
+// splice the right side of the point that is moved
+if (tMove != curve.MaxIndex()) {
+if (r == -1)
+return;
+
+if (trightMove != tMove)
+{
+auto fixedL = curve.Split(tMove, trightMove);
+for (auto i = 1; i < 4; i++)
+splicedCurve.p.push_back(fixedL[i] + v);
+}
+MoveCurveSplice(r, spliceR, rightStep, rextra, false, moveLock, v, influenceRDistance, leftTan, ctrlPtScale, ctrlPtRRotate);
+
+// add the spliced right side of the curve
+for (auto i = splicedCurve.p.size() == 0 ? 0 : 1; i < (r != curve.MaxIndex() ? 3 :4); i++) {
+insertEnd++;
+splicedCurve.p.push_back(spliceR[i]);
+}
+if (r != curve.MaxIndex()) {
+insertEnd++;
+for (auto & p : curve.Split(r, 1.0* curve.MaxIndex()))  // add on the remaining right side of the curve
+splicedCurve.p.push_back(p);
+}
+}
+
+for (auto & pt : splicedCurve.p) {
+if (_isnan(pt.X))
+break;
+}
+
+// adjust all lock t-values based on the size of the inserted splice segments
+for (auto i = 0; i < locked.size(); i++) {
+if (locked[i].T == tMove)
+locked[i].T = moveIndex ==0  ? 0.0 : (moveIndex*1.0-1)/3;
+else if (locked[i].T == l)
+locked[i].T = std::ceil(l);
+else if (locked[i].T == r)
+locked[i].T = (insertEnd*1.0-1)/3;
+else
+locked[i].T = splicedCurve.NearestT(curve[locked[i].T]);
+}
+curve.p = splicedCurve.p;
+}
+
+
+BezierRep                               BezierRep::Rotate(const BezierRep &bez, const double angle, const Point2D &center)
+{
+auto rot = Mat::Rotate(angle);
+auto tri = Mat::Translate(-center);
+auto tr  = Mat::Translate( center);
+BezierRep moved;
+for (auto &p : bez.p) {
+moved.p.push_back(tr * (rot * (tri *p)));
+}
+return moved;
+}
+BezierRep                               BezierRep::Move(const BezierRep &bez, const Vector2D &move)
+{
+BezierRep moved;
+for (auto &p : bez.p) {
+moved.p.push_back(p+move);
+}
+return moved;
+}
+BezierRep                               BezierRep::Interpolate(const BezierRep &start, const BezierRep &end, double t)
+{
+BezierRep interpolated;
+for (auto p=0; p < start.p.size() && p < end.p.size(); p++) {
+interpolated.p.push_back(start.p[p] + (end.p[p]-start.p[p])*t);
+}
+return interpolated;
+}
+std::vector<std::tuple<double, double>> BezierRep::Find_intersections(const BezierRep & a, const BezierRep & b, size_t t_a_off, size_t t_b_off)
+{
+auto ints = std::vector<std::tuple<double, double>>();
+if (a.p.size() == 0 || b.p.size() == 0)
+return ints;
+if (a.p.size() == 4 && b.p.size() == 4)
+{
+std::vector<std::tuple<double, double>> parameters;
+if (SmartRect::Intersect(a.Bounds(), b.Bounds()))
+{
+const int depth = 6;
+Point2D ap[4], bp[4];
+ap[0] = a.p[0];
+ap[1] = a.p[1];
+ap[2] = a.p[2];
+ap[3] = a.p[3];
+bp[0] = b.p[0];
+bp[1] = b.p[1];
+bp[2] = b.p[2];
+bp[3] = b.p[3];
+recursively_intersect(ap, 0, 1, depth, bp, 0, 1, depth, parameters);
+}
+
+std::vector<std::tuple<double, double>> modParameters;
+for (size_t i = 0; i < parameters.size(); i++) {
+modParameters.push_back(std::tuple<double,double>(std::get<0>(parameters[i]) + t_a_off, std::get<1>(parameters[i]) + t_b_off));
+}
+return modParameters;
+}
+for (size_t i = 0; i <= a.p.size() - 4; i += 3)
+{
+for (size_t j = 0; j <= b.p.size() - 4; j += 3)
+{
+std::vector<Point2D> tempVector2(4);
+tempVector2[0] = a.p[i];
+tempVector2[1] = a.p[i + 1];
+tempVector2[2] = a.p[i + 2];
+tempVector2[3] = a.p[i + 3];
+std::vector<Point2D> tempVector3(4);
+tempVector3[0] = b.p[j];
+tempVector3[1] = b.p[j + 1];
+tempVector3[2] = b.p[j + 2];
+tempVector3[3] = b.p[j + 3];
+auto fints = Find_intersections(BezierRep(tempVector2), BezierRep(tempVector3), t_a_off + i / 3, t_b_off + j / 3);
+for (auto inter = 0; inter < fints.size(); inter++) {
+bool newinter = true;
+for (auto & oint : ints)
+if (std::get<0>(oint) == std::get<0>(fints[inter]) &&
+std::get<1>(oint) == std::get<1>(fints[inter])) {
+newinter = false;
+break;
+}
+if (newinter)
+ints.push_back(fints[inter]);
+}
+}
+}
+return ints;
+}
+std::vector<std::vector<Point2D> >      BezierRep::FitCurveSet(    const Point2D d[], size_t dSize, double error, bool & isLoop) {
+std::vector<std::vector<Point2D>> fitSet;
+fitSet.push_back(::FitCurve(d, dSize, error));
+return fitSet;
+}
+std::vector<Point2D>                    BezierRep::FitCurve(       const std::vector<Point2D> &d, double error)
+{
+return ::FitCurve(d.data(), d.size(), error);
+}
+std::vector<Point2D>                    BezierRep::FitOneCurve(const std::vector<Point2D> &d)
+{
+return::FitOneCurve(d.data(), d.size());
+}
+
+std::vector<double>                     BezierRep::Reparameterize( const Point2D d[], size_t first, size_t last, const std::vector<double> &u, const Point2D bezCurve[4])
+{
+std::vector<double> uPrime(last - first + 1); //  New parameter values
+
+for (auto i = first; i <= last; i++)
+{
+uPrime[i - first] = NewtonRaphsonRootFind(bezCurve, d[i], u[i - first]);
+}
+return uPrime;
+}
+double                                  BezierRep::ComputeMaxError(const Point2D d[], size_t first, size_t last, const Point2D bezCurve[4], const std::vector<double> &u, size_t *splitPoint2D)
+{
+double    maxDist; //  Maximum error
+
+if (splitPoint2D)
+*splitPoint2D = (last - first + 1) / 2;
+maxDist      = 0.0;
+for (auto i = first + 1; i < last; i++)
+{
+   double P[2];  //  point on curve
+EvalBezierFast(bezCurve, u[i-first], P);
+double dx = P[0] - d[i].X;//  offset from point to curve
+double dy = P[1] - d[i].Y;
+auto dist = sqrt(dx*dx+dy*dy); //  Current error
+if (dist >= maxDist)
+{
+maxDist = dist;
+if (splitPoint2D)
+*splitPoint2D = i;
+}
+}
+return maxDist;
+}
+std::vector<double>                     BezierRep::ChordLengthParameterize(const Point2D d[], size_t first, size_t last)
+{
+std::vector<double> u(last-first+1);//  Parameterization
+
+double prev = 0.0;
+u[0] = prev;
+for (auto i = first + 1; i <= last; i++)
+{
+auto & lastd = d[i-1];
+auto & curd  = d[i];
+auto dx = lastd.X - curd.X;
+auto dy = lastd.Y - curd.Y;
+prev = u[i - first] = prev + sqrt(dx*dx+dy*dy);
+}
+
+double ulastfirst = u[last-first];
+for (auto i = first + 1; i <= last; i++)
+{
+u[i - first] /= ulastfirst;
+}
+
+return u;
+}
+
+void                                    BezierRep::InsertCpt(double tstart)
+{   auto &allPts  = p;
+    auto t_start_base = (size_t)tstart;
+    if (t_start_base >=  MaxIndex())
+        t_start_base =  MaxIndex() - 1;
+
+    Point2D   left[4], right[4];
+    splitCubic(&allPts[t_start_base*3], tstart - t_start_base, left, right);
+std::vector<Point2D> newP;
+for (size_t i = 0; i < t_start_base*3; i++)
+newP.push_back(allPts[i]);
+for (size_t i = 0; i < 4; i++)
+newP.push_back(left[i]);
+for (size_t i = 1; i < 4; i++)
+newP.push_back(right[i]);
+for (size_t i = t_start_base*3+4; i < allPts.size(); i++)
+newP.push_back(allPts[i]);
+p = newP;
+}
+std::vector<Point2D>                    BezierRep::Split(double tstart, double tend) const
+{
+    auto t_start_base = static_cast<size_t>(tstart);
+    auto t_end_base   = static_cast<size_t>(tend);
+    auto maxIndex      = MaxIndex();
+    if (t_start_base >= maxIndex)
+        t_start_base = maxIndex - 1;
+    if (t_end_base >= maxIndex)
+        t_end_base = maxIndex - 1;
+
+    Point2D              split[4];
+    std::vector<Point2D> splitPts(4);
+    if (t_start_base != t_end_base)
+    {
+bool used4 = true;
+if (tstart - t_start_base == 0) {
+splitPts[0] = p[t_start_base*3];
+splitPts[1] = p[t_start_base*3+1];
+splitPts[2] = p[t_start_base*3+2];
+splitPts[3] = p[t_start_base*3+3];
+} else {
+splitCubic(&(p[t_start_base*3]), tstart - t_start_base, NULL, split);
+if (!(split[0].X == split[1].X && split[1].X == split[2].X && split[2].X == split[3].X &&
+ split[0].Y == split[1].Y && split[1].Y == split[2].Y && split[2].Y == split[3].Y))
+for (size_t i = 0; i < 4; i++)
+splitPts[i] = split[i];
+else {
+splitPts[0] = split[0];
+used4 = false;
+}
+}
+        for (auto i = (t_start_base + 1) * 3; i < t_end_base * 3; i += 3) {
+if (!used4) {
+used4 = true;
+splitPts[1] = p[i+1];
+splitPts[2] = p[i+2];
+splitPts[3] = p[i+3];
+} else {
+splitPts.push_back(p[i+1]);
+splitPts.push_back(p[i+2]);
+splitPts.push_back(p[i+3]);
+}
+}
+        if (t_end_base * 3 < p.size() - 1 && tend - t_end_base != 0)
+        {
+splitCubic(&(p[t_end_base *3]), tend - t_end_base, split, NULL);
+if (!(split[0].X == split[1].X && split[1].X == split[2].X && split[2].X == split[3].X &&
+split[0].Y == split[1].Y && split[1].Y == split[2].Y && split[2].Y == split[3].Y))
+{
+if (!used4) {
+splitPts[1] = split[1];
+splitPts[2] = split[2];
+splitPts[3] = split[3];
+} else {
+splitPts.push_back(split[1]);
+splitPts.push_back(split[2]);
+splitPts.push_back(split[3]);
+}
+
+}
+        }
+    }
+    else
+    {
+Point2D   tmp[4];
+splitCubic(&(p[t_end_base *3]), tend-t_end_base, tmp, NULL);
+splitCubic(tmp, tstart==tend ? 0 : (tstart-t_end_base) / (tend-t_end_base), NULL, split);
+for (auto i = 0; i < 4; i++)
+splitPts[i] = split[i];
+    }
+return splitPts;
+}
+void                                    BezierRep::MoveT(double tMove, const Vector2D & v, bool moveEnds, std::vector<BezierLock> &locked, bool adaptive, bool rangeDrag, double errorTolerance, double influenceLDistance, double influenceRDistance, double ctrlPtLRotate, double ctrlPtRRotate, double ctrlPtScale) {
+if (adaptive)
+MoveTAdaptive(*this, tMove, v, locked, rangeDrag, errorTolerance, influenceLDistance, influenceRDistance, ctrlPtLRotate, ctrlPtRRotate, ctrlPtScale, moveEnds);
+else MoveTAux(*this, tMove, v, moveEnds);
+}
+void                                    BezierRep::GetPoint(const std::vector<Point2D> &p, double t, Point2D &result)
+{
+while (t < 0) {
+t += (p.size()-1)/3;
+}
+while (t > (p.size()-1)/3) {
+t -= (p.size()-1)/3;
+}
+if (p.size() == 0)
+return;
+size_t t_base = 0;
+if (p.size() > 4)
+{
+t_base = static_cast<size_t>(t);
+if (t_base * 3 + 1 >= p.size() - 2) {
+result.X = p.back().X;
+result.Y = p.back().Y;
+return;
+}
+t = t- t_base;
+}
+
+Point2D bez[4] = { p[t_base * 3 + 0], p[t_base * 3 + 1], p[t_base * 3 + 2], p[t_base * 3 + 3] };
+double res[2];
+EvalBezierFast(bez, t, res);
+result.X = res[0];
+result.Y = res[1];
+}
+double                                  BezierRep::NearestT(const Point2D &Pt) const
+{
+if (p.size() < 1)
+return 0;
+
+double closest = DBL_MAX;
+double tclosest = -1;
+for (size_t i = 0; i< MaxIndex(); i++) {
+std::vector<Point2D> tmppts;
+tmppts.push_back(p[i*3]);
+tmppts.push_back(p[i*3+1]);
+tmppts.push_back(p[i*3+2]);
+tmppts.push_back(p[i*3+3]);
+double tc;
+auto nrst = NearestPointOnCurve(Pt, tmppts, &tc);
+if ((nrst-Pt).Length() < closest) {
+closest = (nrst-Pt).Length();
+tclosest = tc+i;
+}
+}
+return tclosest;
+}
+Vector2D                                BezierRep::Tangent(double T) const
+{
+while (T < 0) {
+T += (p.size()-1)/3;
+}
+while (T > (p.size()-1)/3) {
+T -= (p.size()-1)/3;
+}
+if (T == 0)
+{
+for (auto i = 1; i < p.size(); i++)
+if (p[i] != p[0])
+return (p[i]-p[0]).Normal();
+//else return Vector2D();
+return Vector2D();
+}
+if (T == MaxIndex())
+{
+for (int i = static_cast<int>(p.size())-2; i >= 0; i--)
+if (p[i] != p.back())
+return (p.back()-p[i]).Normal();
+//else return Vector2D();
+return Vector2D();
+}
+
+
+int segStart = 3 * (static_cast<int>(T));
+auto t = T - static_cast<int>(T);
+auto A = p[segStart] - p[segStart];
+auto B = p[segStart + 1] - p[segStart];
+// if (B == Vector2D() && segStart > 0 && (p[segStart-1] - p[segStart]) == Vector2D())
+// return Vector2D();
+auto C = p[segStart + 2] - p[segStart];
+auto D = p[segStart + 3] - p[segStart];
+// note that abcd are aka x0 x1 x2 x3
+
+auto tan = -3*A*(1-t)*(1-t) + B*(3*(1-t)*(1-t) - 6*(1 - t)*t) +  C*(6*(1 - t)*t - 3*t*t) + 3*D*t*t;
+return tan.Normal();
+
+// the four coefficients ..
+// A = x3 - 3 * x2 + 3 * x1 - x0
+// B = 3 * x2 - 6 * x1 + 3 * x0
+// C = 3 * x1 - 3 * x0
+// D = x0
+//
+// and then...
+// Vx = 3At2 + 2Bt + C
+
+// first calcuate what are usually know as the coeffients,
+// they are trivial based on the four control points:
+
+//double C1x = (D.X - (3.0 * C.X) + (3.0 * B.X) - A.X);
+//double C2x = ((3.0 * C.X) - (6.0 * B.X) + (3.0 * A.X));
+//double C3x = ((3.0 * B.X) - (3.0 * A.X));
+//double C4x = (A.X); // (not needed for this calculation)
+
+//double C1y = (D.Y - (3.0 * C.Y) + (3.0 * B.Y) - A.Y);
+//double C2y = ((3.0 * C.Y) - (6.0 * B.Y) + (3.0 * A.Y));
+//double C3y = ((3.0 * B.Y) - (3.0 * A.Y));
+//double C4y = (A.Y); // (not needed for this calculation)
+
+// finally it is easy to calculate the slope element, using those coefficients:
+
+//Vector2D vec(((3.0 * C1x * t * t) + (2.0 * C2x * t) + C3x), ((3.0 * C1y * t * t) + (2.0 * C2y * t) + C3y));
+
+//vec.Normalize();
+//return vec;
+// note that this routine works for both the x and y side;
+// simply run this routine twice, once for x once for y
+// note that there are sometimes said to be 8 (not 4) coefficients,
+// these are simply the four for x and four for y, calculated as above in each case.
+}
+bool                                    BezierRep::IsDiscontinuity(int t) const
+{
+if (t == 0 || t == MaxIndex()) {
+if (p.front() != p.back())
+return true;
+
+auto inTan  = (p[1]-p[0]).Normal();
+auto outTan = (p[p.size()-2]-p[p.size()-1]).Normal();
+if (_isnan(inTan.X) || _isnan(outTan.X) || inTan.Dot(outTan) > -0.998)
+return true;
+}
+
+return false;
+}
+Point2D                                 BezierRep::Reflect(const Point2D &srcPt) const
+{
+auto nrstT  = NearestT(srcPt);
+auto nrst   = (*this)[nrstT];
+if (nrstT < 1e-4)
+nrstT = 0;
+if ((MaxIndex()-nrstT) < 1e-4)
+nrstT = static_cast<double>(MaxIndex());
+if (nrstT == 0 || nrstT == MaxIndex() || (p.size()== 4 && p[0]==p[1] && p[2]==p[3])) {
+LnSeg seg(nrst, Tangent(nrstT));
+nrst = seg.LnClosestPoint(srcPt);
+}
+auto normal = Normal(nrstT);
+auto offset = (nrst - srcPt).Length();
+if (normal.Dot(srcPt-nrst) > 0)
+normal = -normal;
+return nrst + normal * offset;
+}
+BezierRep                               BezierRep::Reflect(const BezierRep &b)   const {
+std::vector<Point2D> reflected;
+for (auto &p : b.p) {
+reflected.push_back(Reflect(p));
+}
+return BezierRep(reflected);
+}
+
+//
+// ReflectAndClip - Clips one curve against another, then reflects the clipped segments.
+// This returns two lists of reflected segments corresponding to reflections of segments which were on the same side as the
+// initial point of the stroke (relative to the reflection axis) and those which which were on the opposite side.
+//
+std::vector<std::vector<std::tuple<BezierRep,BezierRep>>>  BezierRep::ReflectAndClip(const BezierRep &b)   const
+{
+BezierRep testRep = *this;
+if (MaxIndex() == 1 && p[0]==p[1] && p[2]==p[3]) {
+Vector2D dir = p[3]-p[0];
+std::vector<Point2D> pts;
+pts.push_back(p[0] - 10000 * dir);
+pts.push_back(p[0] - 10000 * dir);
+pts.push_back(p[3] + 10000 * dir);
+pts.push_back(p[3] + 10000 * dir);
+testRep = BezierRep(pts);
+}
+auto ints = Find_intersections(testRep, b);
+
+std::vector<std::vector<BezierRep>> flipSets;
+std::vector<BezierRep>              fragments[2];
+if (ints.size() == 0) {
+fragments[0].push_back(b);
+flipSets.push_back(fragments[0]);
+} else {
+double start = 0;
+int which = 0;
+for (auto &i: ints) {
+auto split = b.Split(start, std::get<1>(i));
+fragments[which++%2].push_back(split);
+start = std::get<1>(i);
+}
+fragments[which++%2].push_back(b.Split(start, static_cast<double>(b.MaxIndex())));
+
+}
+
+std::vector<std::vector<std::tuple<BezierRep,BezierRep>>> mirroredSides;
+for (auto &side: fragments) {
+std::vector<std::tuple<BezierRep,BezierRep>> mirrors;
+for (auto &f : side)
+mirrors.push_back(std::tuple<BezierRep,BezierRep>(f, Reflect(f)));
+mirroredSides.push_back(mirrors);
+}
+return mirroredSides;
+}
+#ifdef later
+if (!_isnan(influenceDistance) && influenceDistance < 0) {
+spliceL[2] = (spliceL[3] += v);
+if (moveLock) {
+moveLock->Side |= 1;
+moveLock->Cusp = true;
+}
+}
+else if (influenceDistance > 0 && l <= leftStep.T && spliceL[2] != spliceL[3]) {
+auto lTan = (spliceL[2]-spliceL[3]);
+spliceL[2]  = spliceL[3] + Mat::Rotate(ctrlPtRotate) * lTan * ctrlPtScale + v;
+spliceL[3] += v;
+
+LnSeg   tangent(spliceL[3], spliceL[2]), otherTangent(spliceL[0], spliceL[1]);
+auto    inter = otherTangent.LnIntersection(tangent);
+if (inter != Point2D::Null()) {
+auto    tandir = tangent.ClosestFraction(inter) <=0 ? -tangent.Direction() : tangent.Direction();
+auto    aspect = (spliceL[3]-inter).Length() / (spliceL[0]-spliceL[3]).Length() / .7071;
+auto    modinter = spliceL[3] + tandir * (spliceL[3]-inter).Length()*std::min(1.0,.5519/aspect);
+
+auto targetFrac = (modinter-spliceL[3]).Length();
+auto target     = spliceL[3] + targetFrac * tangent.Direction();
+spliceL[2] = (target * std::min(1.0, lextra/25) + spliceL[2] * std::max(0.0, 1-lextra/25));
+}
+} else {
+auto lTan = (spliceL[2]-spliceL[3]);
+if (lTan == Vector2D() && influenceDistance > 0) {
+if (moveLock)
+moveLock->Cusp = false;
+spliceL[2] = spliceL[3] + v - smoothParTangent.Normal()*lextra;
+} else
+spliceL[2]  = spliceL[3] + Mat::Rotate(ctrlPtRotate) * lTan * ctrlPtScale + v;
+spliceL[3] += v;
+}
+#endif
+#if 0
+if (!_isnan(influenceDistance) && influenceDistance < 0) {
+spliceR[1] = (spliceR[0] += v);
+if (moveLock) {
+moveLock->Side |= 2;
+moveLock->Cusp = true;
+}
+}
+else
+if (influenceDistance > 0 && r>=rightStep.T && spliceR[1] != spliceR[0]) {
+
+auto rTan = (spliceR[1]-spliceR[0]);
+spliceR[1]  = spliceR[0] + Mat::Rotate(ctrlPtRotate) * rTan* ctrlPtScale + v;
+spliceR[0] += v;
+
+LnSeg   tangent(spliceR[0], spliceR[1]), otherTangent(spliceR[3], spliceR[2]);
+auto    inter  = otherTangent.LnIntersection(tangent);
+if (inter != Point2D::Null()) {
+auto    tandir = tangent.ClosestFraction(inter) <=0 ? -tangent.Direction() : tangent.Direction();
+//auto    aspect = (spliceR[0]-inter).Length() / (spliceR[3]-inter).Length();
+auto    aspect = (spliceR[0]-inter).Length() / (spliceR[0]-spliceR[3]).Length() / .7071;
+auto    modinter = spliceR[0] + tandir * (spliceR[0]-inter).Length()*std::min(1.0,.5519/aspect);
+
+auto    targetFrac = (modinter-spliceR[0]).Length();
+auto    target     = spliceR[0] + targetFrac * tangent.Direction();
+spliceR[1]      = (target * std::min(1.0, rextra/25)  + spliceR[1] * std::max(0.0, 1-rextra/25));
+}
+} else {
+auto rTan = (spliceR[1]-spliceR[0]);
+if (rTan == Vector2D() && influenceDistance > 0) {
+if (moveLock)
+moveLock->Cusp = false;
+spliceR[1] = spliceR[0] + v + smoothParTangent.Normal()*rextra;
+} else
+spliceR[1]  = spliceR[0] + Mat::Rotate(ctrlPtRotate) * rTan* ctrlPtScale + v;
+spliceR[0] += v;
+}
+#endif
+
+*/
\ No newline at end of file
diff --git a/src/client/views/InkControlPtHandles.tsx b/src/client/views/InkControlPtHandles.tsx
index df803ba31..a91e74c44 100644
--- a/src/client/views/InkControlPtHandles.tsx
+++ b/src/client/views/InkControlPtHandles.tsx
@@ -102,7 +102,7 @@ export class InkControlPtHandles extends React.Component<InkControlProps> {
     @action
     onDelete = (e: KeyboardEvent) => {
         if (["-", "Backspace", "Delete"].includes(e.key)) {
-            InkStrokeProperties.Instance?.deletePoints(this.props.inkView);
+            InkStrokeProperties.Instance?.deletePoints(this.props.inkView, e.shiftKey);
             e.stopPropagation();
         }
     }
diff --git a/src/client/views/InkStrokeProperties.ts b/src/client/views/InkStrokeProperties.ts
index 8ad4864f9..4808bbc77 100644
--- a/src/client/views/InkStrokeProperties.ts
+++ b/src/client/views/InkStrokeProperties.ts
@@ -1,16 +1,18 @@
 import { Bezier } from "bezier-js";
 import { action, observable, reaction } from "mobx";
-import { Doc, Opt, DocListCast } from "../../fields/Doc";
+import { Doc, NumListCast, Opt } from "../../fields/Doc";
 import { InkData, InkField, InkTool, PointData } from "../../fields/InkField";
 import { List } from "../../fields/List";
 import { listSpec } from "../../fields/Schema";
 import { Cast, NumCast } from "../../fields/Types";
+import { Point } from "../../pen-gestures/ndollar";
 import { DocumentType } from "../documents/DocumentTypes";
+import { FitOneCurve } from "../util/bezierFit";
 import { CurrentUserUtils } from "../util/CurrentUserUtils";
+import { DocumentManager } from "../util/DocumentManager";
 import { undoBatch } from "../util/UndoManager";
 import { InkingStroke } from "./InkingStroke";
 import { DocumentView } from "./nodes/DocumentView";
-import { DocumentManager } from "../util/DocumentManager";
 
 export class InkStrokeProperties {
     static Instance: InkStrokeProperties | undefined;
@@ -139,18 +141,35 @@ export class InkStrokeProperties {
      */
     @undoBatch
     @action
-    deletePoints = (inkView: DocumentView) => this.applyFunction(inkView, (view: DocumentView, ink: InkData) => {
+    deletePoints = (inkView: DocumentView, preserve: boolean) => this.applyFunction(inkView, (view: DocumentView, ink: InkData) => {
         const doc = view.rootDoc;
-        const newPoints: { X: number, Y: number }[] = [];
-        const toRemove = Math.floor((this._currentPoint + 2) / 4);
-        const last = this._currentPoint === ink.length - 1;
-        for (let i = 0; i < ink.length; i++) {
-            if (Math.floor((i + 2) / 4) !== toRemove && (toRemove !== 0 || i > 3)) {
-                newPoints.push({ X: ink[i].X, Y: ink[i].Y });
+        const newPoints = ink.slice();
+        const brokenIndices = NumListCast(doc.brokenInkIndices);
+        if (preserve || this._currentPoint === 0 || this._currentPoint === ink.length - 1 || brokenIndices.includes(this._currentPoint)) {
+            newPoints.splice(this._currentPoint === 0 ? 0 : this._currentPoint === ink.length - 1 ? this._currentPoint - 3 : this._currentPoint - 2, 4);
+        } else {
+            const start = this._currentPoint === 0 ? 0 : this._currentPoint - 4;
+            const splicedPoints = ink.slice(start, start + (this._currentPoint === 0 || this._currentPoint === ink.length - 1 ? 4 : 8));
+            var samples: Point[] = [];
+            var startDir = { x: 0, y: 0 };
+            var endDir = { x: 0, y: 0 };
+            for (var i = 0; i < splicedPoints.length / 4; i++) {
+                var bez = new Bezier(splicedPoints.slice(i * 4, i * 4 + 4).map(p => ({ x: p.X, y: p.Y })));
+                if (i === 0) startDir = bez.derivative(0);
+                if (i === splicedPoints.length / 4 - 1) endDir = bez.derivative(1);
+                for (var t = 0; t < (i === splicedPoints.length / 4 - 1 ? 1 + 1e-7 : 1); t += 0.05) {
+                    var pt = bez.compute(t);
+                    samples.push(new Point(pt.x, pt.y));
+                }
+            }
+            const { finalCtrls, error } = FitOneCurve(samples, { X: startDir.x, Y: startDir.y }, { X: endDir.x, Y: endDir.y });
+            if (error < 100) {
+                newPoints.splice(this._currentPoint - 4, 8, ...finalCtrls);
+            } else {
+                newPoints.splice(this._currentPoint - 2, 4);
             }
         }
-        doc.brokenInkIndices = new List(Cast(doc.brokenInkIndices, listSpec("number"), []).map(control => control >= toRemove * 4 ? control - 4 : control));
-        if (last) newPoints.splice(newPoints.length - 3, 2);
+        doc.brokenInkIndices = new List(brokenIndices.map(control => control >= this._currentPoint ? control - 4 : control));
         this._currentPoint = -1;
         return newPoints.length < 4 ? undefined : newPoints;
     }, true)
diff --git a/src/client/views/InkingStroke.tsx b/src/client/views/InkingStroke.tsx
index ecc82a580..8ff080f81 100644
--- a/src/client/views/InkingStroke.tsx
+++ b/src/client/views/InkingStroke.tsx
@@ -132,7 +132,7 @@ export class InkingStroke extends ViewBoxBaseComponent<FieldViewProps, InkDocume
                     } else if (isEditing) {
                         this._nearestT && this._nearestSeg !== undefined && InkStrokeProperties.Instance?.addPoints(this.props.docViewPath().lastElement(), this._nearestT, this._nearestSeg, this.inkScaledData().inkData.slice());
                     }
-                }), isEditing, isEditing, () => wasSelected && InkStrokeProperties.Instance && (InkStrokeProperties.Instance._currentPoint = -1));
+                }), isEditing, isEditing, action(() => wasSelected && InkStrokeProperties.Instance && (InkStrokeProperties.Instance._currentPoint = -1)));
         }
     }