18 #ifndef __UT_FitCubicImpl__
19 #define __UT_FitCubicImpl__
21 #ifndef __UT_FitCubic_H__
22 #error Do not include this file directly. Include UT_FitCubic.h instead.
31 #define ACUTE_ANGLE 1.6
32 #define MAXITERATIONS 4
45 Span *node =
new Span;
49 node->point[0] = fcurve[0];
50 node->point[1] = fcurve[1];
51 node->point[2] = fcurve[2];
52 node->isLinear = islinear;
92 myHead = myHead->next;
100 template <
typename T>
103 bool preserveExtrema)
106 Span *tail=0, *headS, *tailS;
127 if(boss->
opStart(
"Fitting Curve"))
134 findCorner(d, start, last, &corner);
141 nsegs = fitLine(d, start, corner, &headS, &tailS, 1);
143 nsegs = fitLine(d, start, corner, &headS, &tailS, 0);
148 tHat1 = computeCenterTangent(d, start);
150 tHat2 = computeCenterTangent(d, corner);
153 if (closed && (corner==last) && !isSeamCorner(d, last))
157 tHat1 = tHat2*(-1.0);
159 nsegs = fitCubic(d, start, corner, tHat1, tHat2,
160 &headS, &tailS, preserveExtrema, boss);
186 }
while(start < last);
201 template <
typename T>
211 for (i=first+2; i<last-2; ++i)
221 theta = SYSacos(v1.
dot(v2));
235 template <
typename T>
245 if (last<2)
return 1;
250 v2 = d[last-1] - d[0];
253 theta = SYSacos(v1.
dot(v2));
270 template <
typename T>
273 Span **
head, Span **tail,
287 npts = last - first + 1;
292 temp1 = (d[last]-d[
first]);
293 temp2 = (d[
first]-d[last]);
294 fcurve[0] = d[
first];
296 fcurve[1] = fcurve[0] + temp1*
dist;
297 fcurve[2] = fcurve[3] + temp2*
dist;
303 *head = *tail = appendCurve(fcurve, 1);
309 u = chordLengthParameterize(d, first, last);
310 maxError2 = computeLineError(d, first, last, fcurve, u, &maxi);
313 if (maxError2 < myError2)
315 *tail = *head = appendCurve(fcurve, 1);
320 nsegs = fitLine(d, first, maxi, &headL, &tailL, recurse);
321 nsegs += fitLine(d, maxi, last, &headR, &tailR, recurse);
323 if(headL && tailL && headR && tailR)
329 else if(headL && tailL)
356 template <
typename T>
359 Vector2 tHat2, Span **head, Span **tail,
376 npts = last - first + 1;
379 Vector2 tempv(d[first] - d[last]);
380 fpreal dist = tempv.length() / 3.0;
381 rcurve[0] = d[
first];
383 rcurve[1] = rcurve[0] + tHat1*
dist;
384 rcurve[2] = rcurve[3] + tHat2*
dist;
385 *head = *tail = appendCurve(rcurve, 1);
393 u = chordLengthParameterize(d, first, last);
394 generateBezier(d, first, last, u, tHat1, tHat2, rcurve);
395 maxError2 = computeBezierError(d, first, last, rcurve, u, &splitPoint);
399 simpleCurve(d, first, last, tHat1, tHat2, rcurve);
400 maxError2 = computeCubicError(d, first, last, rcurve, &splitPoint,
405 if (maxError2 < myError2)
408 *tail = *head = appendCurve(rcurve, 0);
421 iterationError = myError2 * myError2;
422 if (maxError2 < iterationError)
426 uPrime = reparameterize(d, first, last, u, rcurve);
427 generateBezier(d, first, last, uPrime, tHat1, tHat2, rcurve);
428 maxError2 = computeBezierError(d, first, last, rcurve,
429 uPrime, &splitPoint);
430 if (maxError2 < myError2)
434 *tail = *head = appendCurve(rcurve, 0);
444 if (uPrime)
delete [] uPrime;
455 tHatCenter = computeCenterTangent(d, splitPoint);
457 nsegs = fitCubic(d, first, splitPoint, tHat1, tHatCenter, &headL,
458 &tailL, preserveExtrema, boss);
462 nsegs += fitCubic(d, splitPoint, last, tHatCenter, tHat2, &headR,
463 &tailR, preserveExtrema, boss);
467 if(headL && tailL && headR && tailR)
473 else if(headL && tailL)
508 template <
typename T>
515 fpreal dist = tempv.length() / 3.0;
517 rcurve[0] = d[
first];
519 rcurve[1] = rcurve[0] + tHat1*
dist;
520 rcurve[2] = rcurve[3] + tHat2*
dist;
524 template <
typename T>
528 Vector2 tHat2, CubicCurve rcurve)
541 npts = last - first + 1;
556 for (i=0; i<npts; i++) {
557 A0[i] = tHat1*BEZ_MULT_1(uPrime[i]);
558 A1[i] = tHat2*BEZ_MULT_2(uPrime[i]);
569 for (i = 0; i < npts; i++)
571 C[0][0] += A0[i].
dot(A0[i]);
572 C[0][1] += A0[i].
dot(A1[i]);
574 C[1][1] += A1[i].
dot(A1[i]);
576 tmp = d[last]*BEZ_MULT_3(uPrime[i]);
577 tmp += d[last]*BEZ_MULT_2(uPrime[i]);
578 tmp += d[
first]*BEZ_MULT_1(uPrime[i]);
579 tmp += d[
first]*BEZ_MULT_0(uPrime[i]);
580 tmp = d[first+i] - tmp;
582 X[0] += A0[i].
dot(tmp);
583 X[1] += A1[i].
dot(tmp);
591 det_C0_C1 = C[0][0] * C[1][1] - C[1][0] * C[0][1];
592 det_C0_X = C[0][0] * X[1] - C[0][1] * X[0];
593 det_X_C1 = X[0] * C[1][1] - X[1] * C[0][1];
598 alpha_l = alpha_r = 0;
602 alpha_l = det_X_C1 / det_C0_C1;
603 alpha_r = det_C0_X / det_C0_C1;
606 if (alpha_l < 0.0 || alpha_r < 0.0)
608 simpleCurve(d, first, last, tHat1, tHat2, rcurve);
617 rcurve[0] = d[
first];
619 rcurve[1] = rcurve[0] + tHat1*alpha_l;
620 rcurve[2] = rcurve[3] + tHat2*alpha_r;
632 template <
typename T>
637 int npts = last-first+1;
641 uPrime =
new T[npts];
643 for (i = first; i <= last; i++)
644 uPrime[i-first] = newtonRaphsonRootFind(fcurve, d[i], u[i- first]);
656 template <
typename T>
668 Q_u = calcBezier3(Q, u);
671 for (i = 0; i <= 2; i++)
672 Q1[i] = (Q[i+1] - Q[i]) * 3.0;
675 for (i = 0; i <= 1; i++)
676 Q2[i] = (Q1[i+1] - Q1[i]) * 2.0;
680 Q1_u = calcBezier2(Q1, u);
682 Q2_u = Q2[0]*(1.0-u);
686 numerator = Q1_u.
dot(Q_u - P);
687 denominator = Q1_u.length2() + Q2_u.
dot(Q_u - P);
708 template <
typename T>
726 Vtemp[0] = Vtemp[0]*u;
727 Vtemp[0] += Vtemp[1]*
t;
729 Vtemp[1] = Vtemp[1]*u;
730 Vtemp[1] += Vtemp[2]*
t;
732 Vtemp[0] = Vtemp[0]*u;
733 Vtemp[0] += Vtemp[1]*
t;
739 template <
typename T>
754 Vtemp[0] = Vtemp[0]*u;
755 Vtemp[0] += Vtemp[1]*
t;
760 template <
typename T>
776 ml = (V[1].
y() - V[0].
y()) / (V[1].
x() - V[0].
x());
777 mh = (V[3].
y() - V[2].
y()) / (V[3].
x() - V[2].
x());
788 a = m0 + m1 + 2*(y0 -
y1);
789 b = 3*(y1-
y0) - m1 - 2*m0;
793 x = (t - xl) / (xh - xl);
795 vec.
assign(t, a*x*x*x + b*x*x + c*x + d - (yh-yl));
807 template <
typename T>
818 i1 = myClosed ? (myNumPoints - 1) : 0;
820 if (i2 >= myNumPoints)
821 i2 = myClosed ? 0 : (myNumPoints - 1);
823 tHatCenter = d[i1] - d[i2];
824 tHatCenter.normalize();
836 template <
typename T>
844 u =
new T[last-first+1];
848 for (i = first+1; i <= last; i++)
850 tempv = d[i] - d[i-1];
851 u[i-
first] = u[i-first-1] + tempv.length();
854 for (i = first+1; i <= last; i++)
855 u[i-first] = u[i-first] / u[last-first];
868 template <
typename T>
871 CubicCurve rcurve, T *u,
int *splitPoint)
880 *splitPoint =
int((first+last) * 0.5);
882 for (i = first + 1; i < last; i++)
884 P = calcBezier3(rcurve, u[i-first]);
897 template <
typename T>
900 CubicCurve rcurve,
int *splitPoint,
901 bool preserveExtrema)
911 int typeLocal = d[first+1].y() > localVal ? +1 :
912 (d[first+1].y() < localVal ? -1 : 0);
915 *splitPoint =
int((first+last) * 0.5);
917 for (i = first + 1; i < last; i++)
919 P = calcCubic(rcurve, d[i].
x());
920 dist = SYSpow(P.y() - d[i].y(),
T(2.0));
929 if ((d[i].
y() < localVal && typeLocal == -1)
930 || (d[i].y() > localVal && typeLocal == 1)
931 || (d[i].
y() == localVal && typeLocal == 0))
938 *splitPoint = iLocal;
954 template <
typename T>
957 CubicCurve rcurve, T *u,
int *maxi)
969 for (i = first + 1; i < last; i++)
971 P = rcurve[0]*(1-u[i-
first]);
972 P += rcurve[3]*u[i-
first];
985 #endif // __UT_FitCubicImpl__
GA_API const UT_StringHolder dist
GLboolean GLboolean GLboolean b
UT_StringArray JOINTS head
int opInterrupt(int percent=-1)
GLfloat GLfloat GLfloat v2
INT32 INT32 * denominator
GLboolean GLboolean GLboolean GLboolean a
GLuint GLfloat GLfloat y0
static Vector2 calcCubic(Vector2 *V, fpreal t)
typedef int(WINAPI *PFNWGLRELEASEPBUFFERDCARBPROC)(HPBUFFERARB hPbuffer
GLuint GLfloat GLfloat GLfloat GLfloat y1
bool SYSequalZero(const UT_Vector3T< T > &v)
constexpr SYS_FORCE_INLINE void negate() noexcept
UT_API UT_Interrupt * UTgetInterrupt()
Obtain global UT_Interrupt singleton.
S dot(const V &rhs) const
Return the dot product of two vectors.
int opStart(const char *opname=0, int npoll=0, int immediate_dialog=0, int *opid=0)
void assign(T xx=0.0f, T yy=0.0f)
Set the values of the vector components.
int fitCurve(Vector2 *d, int nPts, int closed, fpreal error2, bool preserveExtrema)
SYS_FORCE_INLINE UT_StorageMathFloat_t< T > normalize() noexcept