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