aboutsummaryrefslogtreecommitdiff
path: root/src/client/util
diff options
context:
space:
mode:
Diffstat (limited to 'src/client/util')
-rw-r--r--src/client/util/bezierFit.ts1446
1 files changed, 1446 insertions, 0 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
generated by cgit v1.2.3-70-g09d2 (git 2.46.0) at 2025-10-07 14:13:48 +0000