/* eslint-disable no-use-before-define */
/* eslint-disable prefer-destructuring */
/* eslint-disable no-param-reassign */
/* eslint-disable camelcase */
import { Point } from '../../pen-gestures/ndollar';

export enum SVGType {
    Rect = 'rect',
    Path = 'path',
    Circle = 'circle',
    Ellipse = 'ellipse',
    Line = 'line',
    Polygon = 'polygon',
}

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) {
        let ret = 0;
        const minx = Math.max(other.TopLeft.X * axis.X + other.TopLeft.Y * axis.Y, this.TopLeft.X * axis.X + this.TopLeft.Y * axis.Y);
        const 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[]) {
        const 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) [r.minx, r.maxx] = [r.maxx, r.minx];
            if (r.miny > r.maxy) [r.miny, r.maxy] = [r.maxy, r.miny];

            p.forEach(pt => {
                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;
    }
}

export function Distance(p: Point) {
    return Math.sqrt(p.X * p.X + p.Y * p.Y);
}
export function Normalize(p: Point) {
    const len = Distance(p);
    return new Point(p.X / len, p.Y / len);
}

function ReparameterizeBezier(d: Point[], first: number, last: number, u: number[], bezCurve: Point[]) {
    const uPrime = new Array<number>(last - first + 1); //  New parameter values

    for (let 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[]) {
    let maxError = 0; //  Maximum error
    let splitPoint2D = (last - first + 1) / 2;
    for (let i = first + 1; i < last; i++) {
        const P = [0, 0]; //  point on curve
        EvalBezierFast(bezCurve, u[i - first], P);
        const dx = P[0] - d[i].X; //  offset from point to curve
        const dy = P[1] - d[i].Y;
        const 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) {
    const u = new Array<number>(last - first + 1); //  Parameterization

    let prev = 0.0;
    u[0] = prev;
    for (let i = first + 1; i <= last; i++) {
        const lastd = d[i - 1];
        const curd = d[i];
        const dx = lastd.X - curd.X;
        const dy = lastd.Y - curd.Y;
        prev = u[i - first] = prev + Math.sqrt(dx * dx + dy * dy);
    }

    const ulastfirst = u[last - first];
    for (let i = first + 1; i <= last; i++) {
        u[i - first] /= ulastfirst;
    }

    return u;
}
/*
 *  B0, B1, B2, B3 :
 *  Bezier multipliers
 */
function B0(u: number) {
    const tmp = 1.0 - u;
    return tmp * tmp * tmp;
}
function B1(u: number) {
    const tmp = 1.0 - u;
    return 3 * u * tmp * tmp;
}
function B2(u: number) {
    const tmp = 1.0 - u;
    return 3 * u * u * tmp;
}
function B3(u: number) {
    return u * u * u;
}
function bounds(p: Point[]) {
    const 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] = [r.maxx, r.minx];
    if (r.miny > r.maxy) [r.miny, r.maxy] = [r.maxy, r.miny]; // swap min & max

    for (let 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[]) {
    const sz = 4;
    const Vtemp = new Array<Array<Point>>(4);
    for (let i = 0; i < 4; i++) Vtemp[i] = new Array<Point>(4);

    /* Copy control points */
    // std::copy(p.begin(), p.end(), Vtemp[0]);
    for (let i = 0; i < sz; i++) {
        Vtemp[0][i].X = p[i].X;
        Vtemp[0][i].Y = p[i].Y;
    }

    /* Triangle computation */
    for (let i = 1; i < sz; i++) {
        for (let j = 0; j < sz - i; j++) {
            const a = Vtemp[i - 1][j];
            const 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 (let j = 0; j < sz; j++) {
            left[j].X = Vtemp[j][0].X;
            left[j].Y = Vtemp[j][0].Y;
        }
    }
    if (right) {
        for (let 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.
 */
// eslint-disable-next-line @typescript-eslint/no-unused-vars
function recursively_intersect(a: Point[], t0: number, t1: number, deptha: number, b: Point[], u0: number, u1: number, depthb: number, parameters: number[][]) {
    if (deptha > 0) {
        const a1 = new Array<Point>(4);
        const a2 = new Array<Point>(4);
        splitCubic(a, 0.5, a1, a2);
        const tmid = (t0 + t1) * 0.5;
        deptha--;
        if (depthb > 0) {
            const b1 = new Array<Point>(4);
            const b2 = new Array<Point>(4);
            splitCubic(b, 0.5, b1, b2);
            const 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) {
        const b1 = new Array<Point>(4);
        const b2 = new Array<Point>(4);
        splitCubic(b, 0.5, b1, b2);
        const 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);
        }
    } // Both segments are fully subdivided; now do line segments
    else {
        const xlk = a[3].X - a[0].X;
        const ylk = a[3].Y - a[0].Y;
        const xnm = b[3].X - b[0].X;
        const ynm = b[3].Y - b[0].Y;
        const xmk = b[0].X - a[0].X;
        const ymk = b[0].Y - a[0].Y;
        const det = xnm * ylk - ynm * xlk;
        if (1.0 + det === 1.0) {
            return;
        }
        const detinv = 1.0 / det;
        const s = (xnm * ymk - ynm * xmk) * detinv;
        const t = (xlk * ymk - ylk * xmk) * detinv;
        if (s < 0.0 || s > 1.0 || t < 0.0 || t > 1.0 || isNaN(s) || 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;
    }

    const 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 (let i = 0; i <= degree; i++) {
        Vtemp[i].X = V[i].X;
        Vtemp[i].Y = V[i].Y;
    }

    /* Triangle computation */
    for (let i = 1; i <= degree; i++) {
        for (let 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[]) {
    const n = 3;
    const u = 1.0 - t;
    let bc = 1;
    let tn = 1;
    let tmpX = p[0].X * u;
    let tmpY = p[0].Y * u;
    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 *= 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) {
    const use = 1;
    const 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) {
    const use = 1;
    const 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);
    }
    const V1 = ComputeLeftTangent(d, center); // d[center] - d[center-1];
    const V2 = ComputeRightTangent(d, center); // d[center] - d[center + 1];
    let 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 */) {
    const nPts = last - first + 1; // Number of pts in sub-curve
    const Ax = new Array<number>(nPts * 2); // Precomputed rhs for eqn  //std::vector<Vector2D> A(nPts * 2);
    const Ay = new Array<number>(nPts * 2); // Precomputed rhs for eqn  //std::vector<Vector2D> A(nPts * 2);

    /* Compute the A's  */
    for (let i = 0; i < nPts; i++) {
        const uprime = uPrime[i];
        const b1 = B1(uprime);
        const 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  */
    const C = [
        [0, 0],
        [0, 0],
    ];
    const df = d[first];
    const dl = d[last];

    const X = [0, 0]; // Matrix X
    for (let 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]);
        const uprime = uPrime[i];
        const b0plb1 = B0(uprime) + B1(uprime);
        const b2plb3 = B2(uprime) + B3(uprime);
        const df1 = d[first + i];
        const tmpX = df1.X - (df.X * b0plb1 + dl.X * b2plb3);
        const 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  */
    const det_C0_C1 = C[0][0] * C[1][1] - C[1][0] * C[0][1] || C[0][0] * C[1][1] * 10e-12;
    const det_C0_X = C[0][0] * X[1] - C[0][1] * X[0];
    const det_X_C1 = X[0] * C[1][1] - X[1] * C[0][1];

    /* Finally, derive alpha values */
    let alpha_l = det_C0_C1 === 0 ? 0.0 : det_X_C1 / det_C0_C1;
    let 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. */
    const segLength = Math.sqrt((df.X - dl.X) * (df.X - dl.X) + (df.Y - dl.Y) * (df.Y - dl.Y));
    const 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) {
    const Q1 = [new Point(0, 0), new Point(0, 0), new Point(0, 0)];
    const Q2 = [new Point(0, 0), new Point(0, 0)]; //  Q' and Q''
    const Q_u = [0, 0];
    const Q1_u = [0, 0];
    const Q2_u = [0, 0]; // u evaluated at Q, Q', & Q''

    /* Compute Q(u) */
    let uPrime: number; //  Improved u
    EvalBezierFast(Q, u, Q_u);

    /* Generate control vertices for Q' */
    for (let 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 (let 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) */
    const numerator = (Q_u[0] - P.X) * Q1_u[0] + (Q_u[1] - P.Y) * Q1_u[1];
    const 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[]) {
    const bezCurve = new Array<Point>(4); //  Control points of fitted Bezier curve
    const maxIterations = 4; //  Max times to try iterating

    const iterationError = error * error; //  Error below which you try iterating
    const nPts = last - first + 1; //  Number of points in subset

    /*  Use heuristic if region only has two points in it */
    if (nPts === 2) {
        const 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 */
    let u = ChordLengthParameterize(d, first, last);
    GenerateBezier(d, first, last, u, tHat1, tHat2, bezCurve);

    /*  Find max deviation of points to fitted curve */
    const { 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 (let i = 0; i < maxIterations; i++) {
            const uPrime = ReparameterizeBezier(d, first, last, u, bezCurve); //  Improved parameter values
            GenerateBezier(d, first, last, uPrime, tHat1, tHat2, bezCurve);
            const { maxError: maximumError } = ComputeMaxError(d, first, last, bezCurve, uPrime);
            if (maximumError < 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 */
    const tHatCenter = splitPoint2D >= last - 1 ? ComputeRightTangent(d, splitPoint2D) : ComputeCenterTangent(d, splitPoint2D);
    FitCubic(d, first, splitPoint2D, tHat1, tHatCenter, error, result);
    const negThatCenter = new Point(-tHatCenter.X, -tHatCenter.Y);
    FitCubic(d, splitPoint2D, last, negThatCenter, tHat2, error, result);
}
/**
 * Convert polyline coordinates to a (multi) segment bezier curve
 * @param d - polyline coordinates
 * @param error - how much error to allow in fitting (measured in pixels)
 * @returns
 */
export function FitCurve(d: Point[], error: number) {
    const tHat1 = ComputeLeftTangent(d, 0); //  Unit tangent vectors at endpoints
    const tHat2 = ComputeRightTangent(d, d.length - 1);
    const result = [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);
    let u = ChordLengthParameterize(d, 0, d.length - 1);
    const 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 */
    let finalCtrls = bezCurveCtrls.slice();
    let { maxError: error } = ComputeMaxError(d, 0, d.length - 1, bezCurveCtrls, u);
    for (let i = 0; i < 10; i++) {
        const uPrime = ReparameterizeBezier(d, 0, d.length - 1, u, bezCurveCtrls); //  Improved parameter values
        GenerateBezier(d, 0, d.length - 1, uPrime, tHat1, tHat2, bezCurveCtrls);
        const { maxError } = ComputeMaxError(d, 0, d.length - 1, bezCurveCtrls, uPrime);
        if (maxError < error) {
            error = maxError;
            finalCtrls = bezCurveCtrls.slice();
        }
        u = uPrime;
    }
    return { finalCtrls, error };
}

// alpha determines how far away the tangents are, or the "tightness" of the bezier
export function GenerateControlPoints(coordinates: Point[], alpha = 0.1) {
    const firstEnd = coordinates.length ? [coordinates[0], coordinates[0]] : [];
    const lastEnd = coordinates.length ? [coordinates.lastElement(), coordinates.lastElement()] : [];
    const points: Point[] = coordinates.slice(1, coordinates.length - 1).flatMap((pt, index, inkData) => {
        const prevPt: Point = index === 0 ? firstEnd[0] : inkData[index - 1];
        const nextPt: Point = index === inkData.length - 1 ? lastEnd[0] : inkData[index + 1];
        if (prevPt.X === nextPt.X) {
            const verticalDist = nextPt.Y - prevPt.Y;
            return [{ X: pt.X, Y: pt.Y - alpha * verticalDist }, pt, pt, { X: pt.X, Y: pt.Y + alpha * verticalDist }];
        } else if (prevPt.Y === nextPt.Y) {
            const horizDist = nextPt.X - prevPt.X;
            return [{ X: pt.X - alpha * horizDist, Y: pt.Y }, pt, pt, { X: pt.X + alpha * horizDist, Y: pt.Y }];
        }
        // tangent vectors between the adjacent points
        const tanX = nextPt.X - prevPt.X;
        const tanY = nextPt.Y - prevPt.Y;
        const ctrlPt1: Point = { X: pt.X - alpha * tanX, Y: pt.Y - alpha * tanY };
        const ctrlPt2: Point = { X: pt.X + alpha * tanX, Y: pt.Y + alpha * tanY };
        return [ctrlPt1, pt, pt, ctrlPt2];
    });
    return [...firstEnd, ...points, ...lastEnd];
}

export function SVGToBezier(name: SVGType, attributes: any): Point[] {
    switch (name) {
        case 'line': {
            const x1 = parseInt(attributes.x1);
            const x2 = parseInt(attributes.x2);
            const y1 = parseInt(attributes.y1);
            const y2 = parseInt(attributes.y2);
            return [
                { X: x1, Y: y1 },
                { X: x1, Y: y1 },
                { X: x2, Y: y2 },
                { X: x2, Y: y2 },
            ];
        }
        case 'circle':
        case 'ellipse': {
            const c = 0.551915024494;
            const centerX = parseInt(attributes.cx);
            const centerY = parseInt(attributes.cy);
            const radiusX = parseInt(attributes.rx) || parseInt(attributes.r);
            const radiusY = parseInt(attributes.ry) || parseInt(attributes.r);
            return [
                { X: centerX, Y: centerY + radiusY },
                { X: centerX + c * radiusX, Y: centerY + radiusY },
                { X: centerX + radiusX, Y: centerY + c * radiusY },
                { X: centerX + radiusX, Y: centerY },
                { X: centerX + radiusX, Y: centerY },
                { X: centerX + radiusX, Y: centerY - c * radiusY },
                { X: centerX + c * radiusX, Y: centerY - radiusY },
                { X: centerX, Y: centerY - radiusY },
                { X: centerX, Y: centerY - radiusY },
                { X: centerX - c * radiusX, Y: centerY - radiusY },
                { X: centerX - radiusX, Y: centerY - c * radiusY },
                { X: centerX - radiusX, Y: centerY },
                { X: centerX - radiusX, Y: centerY },
                { X: centerX - radiusX, Y: centerY + c * radiusY },
                { X: centerX - c * radiusX, Y: centerY + radiusY },
                { X: centerX, Y: centerY + radiusY },
            ];
        }
        case 'rect': {
            const x = parseInt(attributes.x);
            const y = parseInt(attributes.y);
            const width = parseInt(attributes.width);
            const height = parseInt(attributes.height);
            return [
                { X: x, Y: y },
                { X: x, Y: y },
                { X: x + width, Y: y },
                { X: x + width, Y: y },
                { X: x + width, Y: y },
                { X: x + width, Y: y },
                { X: x + width, Y: y + height },
                { X: x + width, Y: y + height },
                { X: x + width, Y: y + height },
                { X: x + width, Y: y + height },
                { X: x, Y: y + height },
                { X: x, Y: y + height },
                { X: x, Y: y + height },
                { X: x, Y: y + height },
                { X: x, Y: y },
                { X: x, Y: y },
            ];
        }
        case 'path': {
            const coordList: Point[] = [];
            const startPt = attributes.d.match(/M(-?\d+\.?\d*),(-?\d+\.?\d*)/);
            coordList.push({ X: parseInt(startPt[1]), Y: parseInt(startPt[2]) });
            const matches: RegExpMatchArray[] = Array.from(
                attributes.d.matchAll(/Q(-?\d+\.?\d*),(-?\d+\.?\d*) (-?\d+\.?\d*),(-?\d+\.?\d*)|C(-?\d+\.?\d*),(-?\d+\.?\d*) (-?\d+\.?\d*),(-?\d+\.?\d*) (-?\d+\.?\d*),(-?\d+\.?\d*)|L(-?\d+\.?\d*),(-?\d+\.?\d*)/g)
            );
            let lastPt: Point = startPt;
            matches.forEach(match => {
                if (match[0].startsWith('Q')) {
                    coordList.push({ X: parseInt(match[1]), Y: parseInt(match[2]) });
                    coordList.push({ X: parseInt(match[1]), Y: parseInt(match[2]) });
                    coordList.push({ X: parseInt(match[3]), Y: parseInt(match[4]) });
                    coordList.push({ X: parseInt(match[3]), Y: parseInt(match[4]) });
                    lastPt = { X: parseInt(match[3]), Y: parseInt(match[4]) };
                } else if (match[0].startsWith('C')) {
                    coordList.push({ X: parseInt(match[5]), Y: parseInt(match[6]) });
                    coordList.push({ X: parseInt(match[7]), Y: parseInt(match[8]) });
                    coordList.push({ X: parseInt(match[9]), Y: parseInt(match[10]) });
                    coordList.push({ X: parseInt(match[9]), Y: parseInt(match[10]) });
                    lastPt = { X: parseInt(match[9]), Y: parseInt(match[10]) };
                } else {
                    coordList.push(lastPt);
                    coordList.push({ X: parseInt(match[11]), Y: parseInt(match[12]) });
                    coordList.push({ X: parseInt(match[11]), Y: parseInt(match[12]) });
                    coordList.push({ X: parseInt(match[11]), Y: parseInt(match[12]) });
                    lastPt = { X: parseInt(match[11]), Y: parseInt(match[12]) };
                }
            });
            const hasZ = attributes.d.match(/Z/);
            if (hasZ) {
                coordList.push(lastPt);
                coordList.push({ X: parseInt(startPt[1]), Y: parseInt(startPt[2]) });
                coordList.push({ X: parseInt(startPt[1]), Y: parseInt(startPt[2]) });
            } else {
                coordList.pop();
            }
            return coordList;
        }
        case 'polygon': {
            const coords: RegExpMatchArray[] = Array.from(attributes.points.matchAll(/(-?\d+\.?\d*),(-?\d+\.?\d*)/g));
            let list: Point[] = [];
            coords.forEach(coord => {
                list.push({ X: parseInt(coord[1]), Y: parseInt(coord[2]) });
                list.push({ X: parseInt(coord[1]), Y: parseInt(coord[2]) });
                list.push({ X: parseInt(coord[1]), Y: parseInt(coord[2]) });
                list.push({ X: parseInt(coord[1]), Y: parseInt(coord[2]) });
            });
            const firstPts = list.splice(0, 2);
            list = list.concat(firstPts);
            return list;
        }
        default:
            return [];
    }
}

/*
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

*/