HDK
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
UT_FitCubicImpl.h
Go to the documentation of this file.
1 /*
2  * PROPRIETARY INFORMATION. This software is proprietary to
3  * Side Effects Software Inc., and is not to be reproduced,
4  * transmitted, or disclosed in any way without written permission.
5  *
6  * NAME: UT_FitCubicImpl.h (C++)
7  *
8  * COMMENTS: Given a group within a detail, fit curves to
9  * each polygon. Original fitting code from "Graphic Gems":
10  * Extended to cubic curves.
11  *
12  * An Algorithm for Automatically Fitting Digitized Curves
13  * by Philip J. Schneider
14  * mostly from "Graphics Gems", Academic Press, 1990
15  *
16  */
17 
18 #ifndef __UT_FitCubicImpl__
19 #define __UT_FitCubicImpl__
20 
21 #ifndef __UT_FitCubic_H__
22 #error Do not include this file directly. Include UT_FitCubic.h instead.
23 #endif
24 
25 
26 #include "UT_Vector2.h"
27 #include "UT_FitCubic.h"
28 #include "UT_Interrupt.h"
29 #include "UT_Spline.h"
30 #include <SYS/SYS_Math.h>
31 
32 #define ACUTE_ANGLE 1.6
33 #define MAXITERATIONS 4
34 
35 
36 
37 /*
38  * allocate a new node in the linked list
39  * containing the non-redundant points of the curve
40  */
41 
42 template <typename T>
43 typename UT_FitCubicT<T>::Span *
44 UT_FitCubicT<T>::appendCurve(CubicCurve fcurve, int islinear)
45 {
46  Span *node = new Span;
47 
48  if (node)
49  {
50  node->point[0] = fcurve[0]; // last point always redundant
51  node->point[1] = fcurve[1];
52  node->point[2] = fcurve[2];
53  node->isLinear = islinear;
54  node->next = 0;
55  }
56 
57  return (node);
58 }
59 
60 
61 /*
62  * FitCurve :
63  * Fit a curve to a set of digitized points
64  *
65  */
66 
67 template <typename T>
69 {
70  myHead = 0;
71  myContainsCurve = 0;
72  myError2 = 1.0;
73  myNSpans = 0;
74  myFitType = UT_FIT_BEZIER;
75  myCurveType = UT_FIT_BOTH;
76 }
77 
78 template <typename T>
80 {
81  destroySolution();
82 }
83 
84 template <typename T>
85 void
87 {
88  Span *temp;
89 
90  while(myHead)
91  {
92  temp = myHead;
93  myHead = myHead->next;
94  delete temp;
95  };
96 
97  myContainsCurve = 0;
98 }
99 
100 
101 template <typename T>
102 int
103 UT_FitCubicT<T>::fitCurve(Vector2 *d, int npts, int closed, fpreal error2,
104  bool preserveExtrema, bool use_kupan_slopes)
105 {
106  Vector2 tHat1, tHat2; /* Unit tangent vectors at endpoints */
107  Span *tail=0, *headS, *tailS;
108  int start, corner;
109  int last = npts - 1;
110  int nsegs;
111  UT_Interrupt *boss;
112 
113 
114  // remove points from corner to corner and fit curve to each segment
115 
116  start = 0;
117 
118  myError2 = error2;
119  myNumPoints = npts;
120  myClosed = closed;
121 
122  myNSpans = 0;
123 
124  destroySolution();
125 
126  boss = UTgetInterrupt();
127 
128  if(boss->opStart("Fitting Curve"))
129  {
130  do
131  {
132  corner = last;
133 
134  if (myCurveType == UT_FIT_POLYS || myCurveType == UT_FIT_BOTH)
135  findCorner(d, start, last, &corner);
136  else
137  corner = last;
138 
139  nsegs = 0;
140 
141  if (myCurveType == UT_FIT_POLYS)
142  nsegs = fitLine(d, start, corner, &headS, &tailS, 1);
143  else if (myCurveType == UT_FIT_BOTH)
144  nsegs = fitLine(d, start, corner, &headS, &tailS, 0);
145 
146  if (!nsegs &&
147  ((myCurveType == UT_FIT_BOTH) || (myCurveType == UT_FIT_CURVES)))
148  {
149  if (use_kupan_slopes)
150  tHat1 = computeMonotonicTangent(d, start);
151  else
152  tHat1 = computeCenterTangent(d, start);
153  tHat1.negate();
154  if (use_kupan_slopes)
155  tHat2 = computeMonotonicTangent(d, corner);
156  else
157  tHat2 = computeCenterTangent(d, corner);
158 
159  // use combination of both if at seam
160  if (closed && (corner==last) && !isSeamCorner(d, last))
161  {
162  tHat2 -= tHat1;
163  tHat2.normalize();
164  tHat1 = tHat2*(-1.0);
165  }
166  nsegs = fitCubic(d, start, corner, tHat1, tHat2,
167  &headS, &tailS, preserveExtrema,
168  use_kupan_slopes, boss);
169  }
170 
171  myNSpans += nsegs;
172 
173 
174  // join with existing curve
175  if(headS && tailS)
176  {
177  if (myHead)
178  {
179  tail->next = headS;
180  tail = tailS;
181  }
182  else
183  {
184  myHead = headS;
185  tail = tailS;
186  }
187  }
188 
189  start = corner;
190 
191  if(boss->opInterrupt())
192  break;
193 
194  } while(start < last);
195  }
196  boss->opEnd();
197 
198  return myNSpans;
199 
200 }
201 
202 
203 /*
204  * Search through the list of points
205  * If there appears to be an acute angle set splitPoint to center
206  *
207  */
208 
209 template <typename T>
210 void
211 UT_FitCubicT<T>::findCorner(Vector2 *d, int first, int last, int *splitPoint)
212 {
213  fpreal theta;
214  Vector2 v1, v2;
215  int i;
216 
217  // prev and next must be 2 steps away from center
218 
219  for (i=first+2; i<last-2; ++i)
220  {
221  // if angle(prev, i, next) < 95 then corner
222 
223  v1 = d[i+2] - d[i];
224  v1.normalize();
225 
226  v2 = d[i-2] - d[i];
227  v2.normalize();
228 
229  theta = SYSacos(v1.dot(v2));
230 
231  if (theta < ACUTE_ANGLE)
232  {
233  *splitPoint = i;
234  return;
235  }
236 
237  }
238 }
239 
240 
241 // check if the closed shape contains a corner at point 0
242 
243 template <typename T>
244 int
246 {
247  fpreal theta;
248  Vector2 v1, v2;
249 
250  // prev and next must be 2 steps away from center
251  // if angle(prev, i, next) < 95 then corner
252 
253  if (last<2) return 1; // not enought points
254 
255  v1 = d[2] - d[0];
256  v1.normalize();
257 
258  v2 = d[last-1] - d[0];
259  v2.normalize();
260 
261  theta = SYSacos(v1.dot(v2));
262 
263  if (theta < ACUTE_ANGLE)
264  return 1;
265  else
266  return 0;
267 }
268 
269 
270 
271 
272 /*
273  * fitLine :
274  * Fit a line segment to a set of digitized points
275  * Return the number of segments required (0 or 1)
276  */
277 
278 template <typename T>
279 int
280 UT_FitCubicT<T>::fitLine(Vector2 *d, int first, int last,
281  Span **head, Span **tail,
282  int recurse)
283 {
284  CubicCurve fcurve; /* Control pts of fitted curve*/
285  T *u = 0; /* Parameter values for point */
286  fpreal maxError2; /* Maximum fitting error */
287  int npts; /* Number of points in subset */
288  Vector2 temp1, temp2; /* tangents in using line approx */
289  fpreal dist;
290  int maxi;
291  Span *headL, *tailL; /* Preceeding subsection */
292  Span *headR, *tailR; /* Proceeding subsection */
293  int nsegs = 0;
294 
295  npts = last - first + 1;
296 
297  /* Use heuristic to to form line */
298 
299  dist = 1.0/3.0;
300  temp1 = (d[last]-d[first]);
301  temp2 = (d[first]-d[last]);
302  fcurve[0] = d[first];
303  fcurve[3] = d[last];
304  fcurve[1] = fcurve[0] + temp1*dist;
305  fcurve[2] = fcurve[3] + temp2*dist;
306 
307  /* if fewer than 2 points use this line */
308 
309  if (npts <= 2)
310  {
311  *head = *tail = appendCurve(fcurve, 1);
312  return 1;
313  }
314 
315  /* else check if line fits well enough */
316 
317  u = chordLengthParameterize(d, first, last);
318  maxError2 = computeLineError(d, first, last, fcurve, u, &maxi);
319  delete [] u;
320 
321  if (maxError2 < myError2)
322  {
323  *tail = *head = appendCurve(fcurve, 1);
324  return 1;
325  }
326  else if (recurse)
327  {
328  nsegs = fitLine(d, first, maxi, &headL, &tailL, recurse);
329  nsegs += fitLine(d, maxi, last, &headR, &tailR, recurse);
330 
331  if(headL && tailL && headR && tailR)
332  {
333  *head = headL;
334  tailL->next = headR;
335  *tail = tailR;
336  }
337  else if(headL && tailL)
338  {
339  *head = headL;
340  tailL->next = 0;
341  *tail = tailL;
342  }
343  else
344  {
345  *head = 0;
346  *tail = 0;
347  }
348  }
349 
350  return nsegs;
351 }
352 
353 
354 
355 
356 /*
357  * FitCubic :
358  * Fit a curve to a (sub)set of digitized points
359  * Return the number of segments required
360  */
361 
362 #include <stdio.h>
363 
364 template <typename T>
365 int
366 UT_FitCubicT<T>::fitCubic(Vector2 *d, int first, int last, Vector2 tHat1,
367  Vector2 tHat2, Span **head, Span **tail,
368  bool preserveExtrema, bool use_kupan_slopes, UT_Interrupt *boss)
369 {
370  CubicCurve rcurve; /* Control pts of fitted curve*/
371  T *u = 0; /* Parameter values for point */
372  fpreal maxError2; /* Maximum fitting error */
373  int splitPoint; /* Point to split point set at */
374  int npts; /* Number of points in subset */
375  fpreal iterationError;/* Error below which try iterating */
376  Vector2 tHatCenter; /* Unit tangent vector at splitPoint */
377  int i;
378  Span *headL, *tailL; /* Preceeding subsection */
379  Span *headR, *tailR; /* Proceeding subsection */
380  int nsegs; /* Number of curves required */
381 
382  /* if fewer than 2 points use heuristic*/
383 
384  npts = last - first + 1;
385  if (npts <= 2)
386  {
387  Vector2 tempv(d[first] - d[last]);
388  fpreal dist = tempv.length() / 3.0;
389  rcurve[0] = d[first];
390  rcurve[3] = d[last];
391  rcurve[1] = rcurve[0] + tHat1*dist;
392  rcurve[2] = rcurve[3] + tHat2*dist;
393  *head = *tail = appendCurve(rcurve, 1);
394  return 1;
395  }
396 
397  /* attempt to fit curve */
398 
399  if (myFitType == UT_FIT_BEZIER)
400  {
401  u = chordLengthParameterize(d, first, last);
402  generateBezier(d, first, last, u, tHat1, tHat2, rcurve);
403  maxError2 = computeBezierError(d, first, last, rcurve, u, &splitPoint);
404  }
405  else
406  {
407  simpleCurve(d, first, last, tHat1, tHat2, rcurve);
408  maxError2 = computeCubicError(d, first, last, rcurve, &splitPoint,
409  preserveExtrema);
410  }
411 
412 
413  if (maxError2 < myError2)
414  {
415  if (u) delete [] u; // RB: stop memory leak
416  *tail = *head = appendCurve(rcurve, 0);
417  myContainsCurve = 1;
418  return 1;
419  }
420 
421 
422  /* If error not too large, try some reparameterization */
423  /* and iteration */
424 
425  if (myFitType == UT_FIT_BEZIER)
426  {
427  T *uPrime = 0; /* Improved parameter values */
428 
429  iterationError = myError2 * myError2;
430  if (maxError2 < iterationError)
431  {
432  for (i = 0; i < MAXITERATIONS; i++)
433  {
434  uPrime = reparameterize(d, first, last, u, rcurve);
435  generateBezier(d, first, last, uPrime, tHat1, tHat2, rcurve);
436  maxError2 = computeBezierError(d, first, last, rcurve,
437  uPrime, &splitPoint);
438  if (maxError2 < myError2)
439  {
440  delete [] u; // RB: stop memory leaks
441  delete [] uPrime;
442  *tail = *head = appendCurve(rcurve, 0);
443  myContainsCurve = 1;
444  return 1;
445  }
446 
447  delete [] u;
448  u = uPrime;
449  }
450  }
451 
452  if (uPrime) delete [] uPrime;
453  }
454  else
455  {
456  delete [] u;
457  }
458 
459 
460  /* Fitting failed -- split at max error point and fit recursively */
461  if(!boss->opInterrupt())
462  {
463  if (use_kupan_slopes)
464  tHatCenter = computeMonotonicTangent(d, splitPoint);
465  else
466  tHatCenter = computeCenterTangent(d, splitPoint);
467 
468  nsegs = fitCubic(d, first, splitPoint, tHat1, tHatCenter, &headL,
469  &tailL, preserveExtrema, use_kupan_slopes, boss);
470 
471  tHatCenter *= -1.0;
472 
473  nsegs += fitCubic(d, splitPoint, last, tHatCenter, tHat2, &headR,
474  &tailR, preserveExtrema, use_kupan_slopes, boss);
475 
476  /* Concatenate two pieces together */
477 
478  if(headL && tailL && headR && tailR)
479  {
480  *head = headL;
481  tailL->next = headR;
482  *tail = tailR;
483  }
484  else if(headL && tailL)
485  {
486  *head = headL;
487  tailL->next = 0;
488  *tail = tailL;
489  }
490  else
491  {
492  *head = 0;
493  *tail = 0;
494  }
495  }
496  else
497  {
498  nsegs = 0;
499  tail =0;
500  head =0;
501  }
502 
503  return nsegs;
504 
505 }
506 
507 
508 /*
509  * generateBezier :
510  * Use least-squares method to find Bezier control points for region.
511  *
512  */
513 
514 
515 /*
516  * use the Wu/Barsky heuristic
517  */
518 
519 template <typename T>
520 void
521 UT_FitCubicT<T>::simpleCurve(Vector2 *d, int first, int last,
522  Vector2 tHat1, Vector2 tHat2,
523  CubicCurve rcurve)
524 {
525  Vector2 tempv = Vector2(d[first] - d[last]);
526  fpreal dist = tempv.length() / 3.0;
527 
528  rcurve[0] = d[first];
529  rcurve[3] = d[last];
530  rcurve[1] = rcurve[0] + tHat1*dist;
531  rcurve[2] = rcurve[3] + tHat2*dist;
532  return;
533 }
534 
535 template <typename T>
536 void
537 UT_FitCubicT<T>::generateBezier(Vector2 *d, int first, int last,
538  T *uPrime, Vector2 tHat1,
539  Vector2 tHat2, CubicCurve rcurve)
540 {
541  int i;
542  Vector2 *A0, *A1; /* Precomputed rhs for eqn */
543  int npts; /* Number of pts in sub-curve */
544  T C[2][2]; /* Matrix C */
545  T X[2]; /* Matrix X */
546  fpreal det_C0_C1, /* Determinants of matrices */
547  det_C0_X,
548  det_X_C1;
549  fpreal alpha_l, alpha_r; /* Alpha values, left and right */
550  Vector2 tmp; /* Utility variable */
551 
552  npts = last - first + 1;
553 
554  /* Allocate 2-columned A array */
555 
556  A0 = new Vector2[npts];
557  if (!A0) return;
558 
559  A1 = new Vector2[npts];
560  if (!A1)
561  {
562  delete [] A0;
563  return;
564  }
565 
566  /* Compute the A's */
567  for (i=0; i<npts; i++) {
568  A0[i] = tHat1*BEZ_MULT_1(uPrime[i]);
569  A1[i] = tHat2*BEZ_MULT_2(uPrime[i]);
570  }
571 
572  /* Create the C and X matrices */
573  C[0][0] = 0.0;
574  C[0][1] = 0.0;
575  C[1][0] = 0.0;
576  C[1][1] = 0.0;
577  X[0] = 0.0;
578  X[1] = 0.0;
579 
580  for (i = 0; i < npts; i++)
581  {
582  C[0][0] += A0[i].dot(A0[i]);
583  C[0][1] += A0[i].dot(A1[i]);
584  C[1][0] = C[0][1];
585  C[1][1] += A1[i].dot(A1[i]);
586 
587  tmp = d[last]*BEZ_MULT_3(uPrime[i]);
588  tmp += d[last]*BEZ_MULT_2(uPrime[i]);
589  tmp += d[first]*BEZ_MULT_1(uPrime[i]);
590  tmp += d[first]*BEZ_MULT_0(uPrime[i]);
591  tmp = d[first+i] - tmp;
592 
593  X[0] += A0[i].dot(tmp);
594  X[1] += A1[i].dot(tmp);
595  }
596 
597  /* Finished with A array */
598  delete [] A0;
599  delete [] A1;
600 
601  /* Compute the determinants of C and X */
602  det_C0_C1 = C[0][0] * C[1][1] - C[1][0] * C[0][1];
603  det_C0_X = C[0][0] * X[1] - C[0][1] * X[0];
604  det_X_C1 = X[0] * C[1][1] - X[1] * C[0][1];
605 
606  /* Finally, derive alpha values */
607  if (SYSequalZero(det_C0_C1))
608  {
609  alpha_l = alpha_r = 0;
610  }
611  else
612  {
613  alpha_l = det_X_C1 / det_C0_C1;
614  alpha_r = det_C0_X / det_C0_C1;
615  }
616 
617  if (alpha_l < 0.0 || alpha_r < 0.0)
618  {
619  simpleCurve(d, first, last, tHat1, tHat2, rcurve);
620  }
621  else
622  {
623  /* First and last control points of the curve are */
624  /* positioned exactly at the first and last data points */
625  /* Control points 1 and 2 are positioned an alpha distance out */
626  /* on the tangent vectors, left and right, respectively */
627 
628  rcurve[0] = d[first];
629  rcurve[3] = d[last];
630  rcurve[1] = rcurve[0] + tHat1*alpha_l;
631  rcurve[2] = rcurve[3] + tHat2*alpha_r;
632  }
633 }
634 
635 
636 /*
637  * reparameterize:
638  * Given set of points and their parameterization, try to find
639  * a better parameterization.
640  *
641  */
642 
643 template <typename T>
644 T *
645 UT_FitCubicT<T>::reparameterize(Vector2 *d, int first, int last, T *u,
646  CubicCurve fcurve)
647 {
648  int npts = last-first+1;
649  int i;
650  T *uPrime; /* New parameter values */
651 
652  uPrime = new T[npts];
653 
654  for (i = first; i <= last; i++)
655  uPrime[i-first] = newtonRaphsonRootFind(fcurve, d[i], u[i- first]);
656 
657  return (uPrime);
658 }
659 
660 
661 
662 /*
663  * newtonRaphsonRootFind :
664  * Use Newton-Raphson iteration to find better root.
665  */
666 
667 template <typename T>
668 fpreal
670 {
671  fpreal numerator, denominator;
672  Vector2 Q1[3], Q2[2]; /* Q' and Q'' */
673  Vector2 Q_u, Q1_u, Q2_u; /*u evaluated at Q, Q', & Q'' */
674  fpreal uPrime; /* Improved u */
675  int i;
676 
677  /* Compute Q(u) */
678  if (myFitType == UT_FIT_BEZIER)
679  Q_u = calcBezier3(Q, u);
680 
681  /* Generate control vertices for Q' */
682  for (i = 0; i <= 2; i++)
683  Q1[i] = (Q[i+1] - Q[i]) * 3.0;
684 
685  /* Generate control vertices for Q'' */
686  for (i = 0; i <= 1; i++)
687  Q2[i] = (Q1[i+1] - Q1[i]) * 2.0;
688 
689  /* Compute Q'(u) and Q''(u) */
690  if (myFitType == UT_FIT_BEZIER)
691  Q1_u = calcBezier2(Q1, u);
692 
693  Q2_u = Q2[0]*(1.0-u);
694  Q2_u += Q2[1]*u;
695 
696  /* Compute f(u)/f'(u) */
697  numerator = Q1_u.dot(Q_u - P);
698  denominator = Q1_u.length2() + Q2_u.dot(Q_u - P);
699 
700 
701  /* u = u - f(u)/f'(u) */
702 
703  if (!SYSequalZero(denominator))
704  uPrime = u - (numerator/denominator);
705  else
706  uPrime = u;
707 
708  return (uPrime);
709 }
710 
711 
712 
713 /*
714  * Bezier :
715  * Evaluate a Bezier curve at a particular parameter value
716  *
717  */
718 
719 template <typename T>
720 typename UT_FitCubicT<T>::Vector2
722 {
723  fpreal u = 1.0 - t;
724  Vector2 Vtemp[4]; /* Local copy of control points */
725 
726  /* Triangle computation */
727 
728  Vtemp[0] = V[0]*u;
729  Vtemp[0] += V[1]*t;
730 
731  Vtemp[1] = V[1]*u;
732  Vtemp[1] += V[2]*t;
733 
734  Vtemp[2] = V[2]*u;
735  Vtemp[2] += V[3]*t;
736 
737  Vtemp[0] = Vtemp[0]*u;
738  Vtemp[0] += Vtemp[1]*t;
739 
740  Vtemp[1] = Vtemp[1]*u;
741  Vtemp[1] += Vtemp[2]*t;
742 
743  Vtemp[0] = Vtemp[0]*u;
744  Vtemp[0] += Vtemp[1]*t;
745 
746  return Vtemp[0];
747 }
748 
749 
750 template <typename T>
751 typename UT_FitCubicT<T>::Vector2
753 {
754  fpreal u = 1.0 - t;
755  Vector2 Vtemp[2];
756 
757  /* Triangle computation */
758 
759  Vtemp[0] = V[0]*u;
760  Vtemp[0] += V[1]*t;
761 
762  Vtemp[1] = V[1]*u;
763  Vtemp[1] += V[2]*t;
764 
765  Vtemp[0] = Vtemp[0]*u;
766  Vtemp[0] += Vtemp[1]*t;
767 
768  return Vtemp[0];
769 }
770 
771 template <typename T>
772 typename UT_FitCubicT<T>::Vector2
774 {
775  fpreal a, b, c, d;
776  fpreal y0, y1, m0, m1;
777  fpreal x, ml, mh;
778  fpreal xl, xh, yh, yl;
779  Vector2 vec;
780 
781  xl = V[0].x();
782  xh = V[3].x();
783 
784  yl = V[0].y();
785  yh = V[3].y();
786 
787  ml = (V[1].y() - V[0].y()) / (V[1].x() - V[0].x());
788  mh = (V[3].y() - V[2].y()) / (V[3].x() - V[2].x());
789 
790  // m0 = ml * (xh-xl)/(yh-yl);
791  // m1 = mh * (xh-xl)/(yh-yl);
792 
793  m0 = ml * (xh-xl);
794  m1 = mh * (xh-xl);
795 
796  y0 = yl;
797  y1 = yh;
798 
799  a = m0 + m1 + 2*(y0 - y1);
800  b = 3*(y1-y0) - m1 - 2*m0;
801  c = m0;
802  d = y1;
803 
804  x = (t - xl) / (xh - xl);
805 
806  vec.assign(t, a*x*x*x + b*x*x + c*x + d - (yh-yl));
807 
808  return vec;
809 }
810 
811 
812 /*
813  * computeCenterTangent :
814  * Approximate unit tangents at "center" of digitized curve
815  */
816 
817 
818 template <typename T>
819 typename UT_FitCubicT<T>::Vector2
821 {
822  Vector2 tHatCenter;
823  int i1, i2;
824 
825  i1 = center - 1;
826  i2 = center + 1;
827 
828  if (i1 < 0)
829  i1 = myClosed ? (myNumPoints - 1) : 0;
830 
831  if (i2 >= myNumPoints)
832  i2 = myClosed ? 0 : (myNumPoints - 1);
833 
834  tHatCenter = d[i1] - d[i2];
835  tHatCenter.normalize();
836 
837  return tHatCenter;
838 }
839 
840 template <typename T>
841 typename UT_FitCubicT<T>::Vector2
843 {
844  Vector2 tangent;
845  int i1 = center - 1;
846  int i2 = center + 1;
847 
848  if (i1 < 0)
849  {
850  if (i2 >= myNumPoints)
851  tangent = Vector2(1, 0);
852  else
853  tangent = (d[center] - d[i2]);
854  }
855  else if (i2 >= myNumPoints)
856  {
857  tangent = (d[i1] - d[center]);
858  }
859  else
860  {
861  fpreal64 v0, v1, v2;
862  fpreal64 t0, t1, t2;
863 
864  v0 = d[i1].y();
865  v1 = d[center].y();
866  v2 = d[i2].y();
867  t0 = d[i1].x();
868  t1 = d[center].x();
869  t2 = d[i2].x();
870 
871  fpreal64 slope = UT_Spline::getMonotoneSlopePA(v1, t1, v0, t0, v2, t2);
872  tangent = Vector2(1, slope);
873  }
874 
875  tangent.normalize();
876  return tangent;
877 }
878 
879 
880 /*
881  * chordLengthParameterize :
882  * Assign parameter values to digitized points
883  * using relative distances between points.
884  */
885 
886 template <typename T>
887 T *
888 UT_FitCubicT<T>::chordLengthParameterize(Vector2 *d, int first, int last)
889 {
890  int i;
891  T *u; /* Parameterization */
892  Vector2 tempv;
893 
894  u = new T[last-first+1];
895 
896  u[0] = 0.0;
897 
898  for (i = first+1; i <= last; i++)
899  {
900  tempv = d[i] - d[i-1];
901  u[i-first] = u[i-first-1] + tempv.length();
902  }
903 
904  for (i = first+1; i <= last; i++)
905  u[i-first] = u[i-first] / u[last-first];
906 
907  return(u);
908 }
909 
910 
911 
912 /*
913  * computeError :
914  * Find the maximum squared distance of digitized points
915  * to fitted curve.
916  */
917 
918 template <typename T>
919 fpreal
920 UT_FitCubicT<T>::computeBezierError(Vector2 *d, int first, int last,
921  CubicCurve rcurve, T *u, int *splitPoint)
922 {
923  int i;
924  fpreal maxDist; /* Maximum error */
925  fpreal dist; /* Current error */
926  Vector2 P; /* Point on curve */
927  Vector2 v; /* Vector from point to curve */
928 
929  maxDist = 0.0;
930  *splitPoint = int((first+last) * 0.5);
931 
932  for (i = first + 1; i < last; i++)
933  {
934  P = calcBezier3(rcurve, u[i-first]);
935  v = P - d[i];
936  dist = v.length2();
937  if (dist >= maxDist)
938  {
939  maxDist = dist;
940  *splitPoint = i;
941  }
942  }
943 
944  return maxDist;
945 }
946 
947 template <typename T>
948 fpreal
949 UT_FitCubicT<T>::computeCubicError(Vector2 *d, int first, int last,
950  CubicCurve rcurve, int *splitPoint,
951  bool preserveExtrema)
952 {
953  int i;
954  fpreal maxDist; /* Maximum error */
955  fpreal dist; /* Current error */
956  Vector2 P; /* Point on curve */
957 
958  // Local min/max variables
959  int iLocal = first;
960  fpreal localVal = d[first].y();
961  int typeLocal = SYSisGreater(d[first+1].y(), localVal) ? +1 :
962  (SYSisLess(d[first+1].y(), localVal) ? -1 : 0);
963 
964  maxDist = 0.0;
965  *splitPoint = int((first+last) * 0.5);
966 
967  for (i = first + 1; i < last; i++)
968  {
969  P = calcCubic(rcurve, d[i].x());
970  dist = SYSpow(P.y() - d[i].y(), T(2.0));
971  if (dist >= maxDist)
972  {
973  maxDist = dist;
974  *splitPoint = i;
975  }
976 
977  if (preserveExtrema)
978  {
979  if ((SYSisLess(d[i].y(), localVal) && typeLocal == -1)
980  || (SYSisGreater(d[i].y(), localVal) && typeLocal == 1)
981  || (SYSisEqual(localVal, d[i].y()) && typeLocal == 0))
982  {
983  localVal = d[i].y();
984  iLocal = i;
985  }
986  else
987  {
988  *splitPoint = iLocal;
989  return myError2;
990  }
991  }
992  }
993 
994  return maxDist;
995 }
996 
997 
998 /*
999  * computeLineError :
1000  * Find the maximum squared distance of digitized points
1001  * to fitted line.
1002  */
1003 
1004 template <typename T>
1005 fpreal
1006 UT_FitCubicT<T>::computeLineError(Vector2 *d, int first, int last,
1007  CubicCurve rcurve, T *u, int *maxi)
1008 {
1009  int i;
1010  fpreal maxDist; /* Maximum error */
1011  fpreal dist; /* Current error */
1012  Vector2 P; /* Point on line */
1013  Vector2 v; /* Vector from point to line */
1014 
1015  maxDist = 0.0;
1016  if (maxi)
1017  *maxi = first;
1018 
1019  for (i = first + 1; i < last; i++)
1020  {
1021  P = rcurve[0]*(1-u[i-first]);
1022  P += rcurve[3]*u[i-first];
1023  v = P - d[i];
1024  dist = v.length2();
1025  if (dist >= maxDist)
1026  {
1027  maxDist = dist;
1028  if (maxi)
1029  *maxi = i;
1030  }
1031  }
1032  return (maxDist);
1033 }
1034 
1035 #endif // __UT_FitCubicImpl__
GLint first
Definition: glcorearb.h:405
typedef int(APIENTRYP RE_PFNGLXSWAPINTERVALSGIPROC)(int)
GA_API const UT_StringHolder dist
#define MAXITERATIONS
const GLdouble * v
Definition: glcorearb.h:837
UT_StringArray JOINTS head
GLuint start
Definition: glcorearb.h:475
GLint GLint i2
Definition: glad.h:2724
GLboolean GLboolean GLboolean GLboolean a
Definition: glcorearb.h:1222
X
Definition: ImathEuler.h:183
GLint y
Definition: glcorearb.h:103
GLfloat GLfloat GLfloat v2
Definition: glcorearb.h:818
2D Vector class.
Definition: UT_Vector2.h:159
double fpreal64
Definition: SYS_Types.h:201
constexpr SYS_FORCE_INLINE T & x() noexcept
Definition: UT_Vector2.h:423
int opInterrupt(int percent=-1)
GLdouble y1
Definition: glad.h:2349
static Vector2 calcCubic(Vector2 *V, fpreal t)
Definition: Types.h:285
GLint i1
Definition: glad.h:2724
GLboolean GLboolean GLboolean b
Definition: glcorearb.h:1222
GLint GLenum GLint x
Definition: glcorearb.h:409
static fpreal64 getMonotoneSlopePA(fpreal64 v_cur, fpreal64 t_cur, fpreal64 v_prev, fpreal64 t_prev, fpreal64 v_next, fpreal64 t_next)
GLdouble t
Definition: glad.h:2397
GLfloat v0
Definition: glcorearb.h:816
bool SYSequalZero(const UT_Vector3T< T > &v)
Definition: UT_Vector3.h:1069
__hostdev__ uint64_t last(uint32_t i) const
Definition: NanoVDB.h:5976
constexpr SYS_FORCE_INLINE void negate() noexcept
Definition: UT_Vector2.h:284
fpreal64 fpreal
Definition: SYS_Types.h:277
UT_API UT_Interrupt * UTgetInterrupt()
Obtain global UT_Interrupt singleton.
void opEnd(int opid=-1)
GLfloat GLfloat v1
Definition: glcorearb.h:817
S dot(const V &rhs) const
Return the dot product of two vectors.
Definition: Types.h:229
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.
Definition: UT_Vector2.h:446
int fitCurve(Vector2 *d, int nPts, int closed, fpreal error2, bool preserveExtrema, bool use_kupan_slopes=false)
#define ACUTE_ANGLE
bool SYSisEqual(const UT_Vector2T< T > &a, const UT_Vector2T< T > &b, S tol=SYS_FTOLERANCE)
Componentwise equality.
Definition: UT_Vector2.h:674
SYS_FORCE_INLINE UT_StorageMathFloat_t< T > normalize() noexcept
Definition: UT_Vector2.h:309
constexpr SYS_FORCE_INLINE T & y() noexcept
Definition: UT_Vector2.h:425