API tools faq
paste
Advertisement
Guest User

Untitled

a guest
Jul 24th, 2013
346
0
Never
Add comment
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
text 14.61 KB | None | 0 0
raw download clone embed print report
  1. /*
  2. An Algorithm for Automatically Fitting Digitized Curves
  3. by Philip J. Schneider
  4. from "Graphics Gems", Academic Press, 1990
  5. */
  6. public static class FitCurves
  7. {
  8. /* Fit the Bezier curves */
  9.  
  10. private const int MAXPOINTS = 10000;
  11. public static List<Point> FitCurve(Point[] d, double error)
  12. {
  13. Vector tHat1, tHat2; /* Unit tangent vectors at endpoints */
  14.  
  15. tHat1 = ComputeLeftTangent(d, 0);
  16. tHat2 = ComputeRightTangent(d, d.Length - 1);
  17. List<Point> result = new List<Point>();
  18. FitCubic(d, 0, d.Length - 1, tHat1, tHat2, error,result);
  19. return result;
  20. }
  21.  
  22.  
  23. private static void FitCubic(Point[] d, int first, int last, Vector tHat1, Vector tHat2, double error,List<Point> result)
  24. {
  25. Point[] bezCurve; /*Control points of fitted Bezier curve*/
  26. double[] u; /* Parameter values for point */
  27. double[] uPrime; /* Improved parameter values */
  28. double maxError; /* Maximum fitting error */
  29. int splitPoint; /* Point to split point set at */
  30. int nPts; /* Number of points in subset */
  31. double iterationError; /*Error below which you try iterating */
  32. int maxIterations = 4; /* Max times to try iterating */
  33. Vector tHatCenter; /* Unit tangent vector at splitPoint */
  34. int i;
  35.  
  36. iterationError = error * error;
  37. nPts = last - first + 1;
  38.  
  39. /* Use heuristic if region only has two points in it */
  40. if(nPts == 2)
  41. {
  42. double dist = (d[first]-d[last]).Length / 3.0;
  43.  
  44. bezCurve = new Point[4];
  45. bezCurve[0] = d[first];
  46. bezCurve[3] = d[last];
  47. bezCurve[1] = (tHat1 * dist) + bezCurve[0];
  48. bezCurve[2] = (tHat2 * dist) + bezCurve[3];
  49.  
  50. result.Add(bezCurve[1]);
  51. result.Add(bezCurve[2]);
  52. result.Add(bezCurve[3]);
  53. return;
  54. }
  55.  
  56. /* Parameterize points, and attempt to fit curve */
  57. u = ChordLengthParameterize(d, first, last);
  58. bezCurve = GenerateBezier(d, first, last, u, tHat1, tHat2);
  59.  
  60. /* Find max deviation of points to fitted curve */
  61. maxError = ComputeMaxError(d, first, last, bezCurve, u,out splitPoint);
  62. if(maxError < error)
  63. {
  64. result.Add(bezCurve[1]);
  65. result.Add(bezCurve[2]);
  66. result.Add(bezCurve[3]);
  67. return;
  68. }
  69.  
  70.  
  71. /* If error not too large, try some reparameterization */
  72. /* and iteration */
  73. if(maxError < iterationError)
  74. {
  75. for(i = 0; i < maxIterations; i++)
  76. {
  77. uPrime = Reparameterize(d, first, last, u, bezCurve);
  78. bezCurve = GenerateBezier(d, first, last, uPrime, tHat1, tHat2);
  79. maxError = ComputeMaxError(d, first, last,
  80. bezCurve, uPrime,out splitPoint);
  81. if(maxError < error)
  82. {
  83. result.Add(bezCurve[1]);
  84. result.Add(bezCurve[2]);
  85. result.Add(bezCurve[3]);
  86. return;
  87. }
  88. u = uPrime;
  89. }
  90. }
  91.  
  92. /* Fitting failed -- split at max error point and fit recursively */
  93. tHatCenter = ComputeCenterTangent(d, splitPoint);
  94. FitCubic(d, first, splitPoint, tHat1, tHatCenter, error,result);
  95. tHatCenter.Negate();
  96. FitCubic(d, splitPoint, last, tHatCenter, tHat2, error,result);
  97. }
  98.  
  99. static Point[] GenerateBezier(Point[] d, int first, int last, double[] uPrime, Vector tHat1, Vector tHat2)
  100. {
  101. int i;
  102. Vector[,] A = new Vector[MAXPOINTS,2];/* Precomputed rhs for eqn */
  103.  
  104. int nPts; /* Number of pts in sub-curve */
  105. double[,] C = new double[2,2]; /* Matrix C */
  106. double[] X = new double[2]; /* Matrix X */
  107. double det_C0_C1, /* Determinants of matrices */
  108. det_C0_X,
  109. det_X_C1;
  110. double alpha_l, /* Alpha values, left and right */
  111. alpha_r;
  112. Vector tmp; /* Utility variable */
  113. Point[] bezCurve = new Point[4]; /* RETURN bezier curve ctl pts */
  114. nPts = last - first + 1;
  115.  
  116. /* Compute the A's */
  117. for (i = 0; i < nPts; i++) {
  118. Vector v1, v2;
  119. v1 = tHat1;
  120. v2 = tHat2;
  121. v1 *= B1(uPrime[i]);
  122. v2 *= B2(uPrime[i]);
  123. A[i,0] = v1;
  124. A[i,1] = v2;
  125. }
  126.  
  127. /* Create the C and X matrices */
  128. C[0,0] = 0.0;
  129. C[0,1] = 0.0;
  130. C[1,0] = 0.0;
  131. C[1,1] = 0.0;
  132. X[0] = 0.0;
  133. X[1] = 0.0;
  134.  
  135. for (i = 0; i < nPts; i++) {
  136. C[0,0] += V2Dot(A[i,0], A[i,0]);
  137. C[0,1] += V2Dot(A[i,0], A[i,1]);
  138. /* C[1][0] += V2Dot(&A[i][0], &A[i][9]);*/
  139. C[1,0] = C[0,1];
  140. C[1,1] += V2Dot(A[i,1], A[i,1]);
  141.  
  142. tmp = ((Vector)d[first + i] -
  143. (
  144. ((Vector)d[first] * B0(uPrime[i])) +
  145. (
  146. ((Vector)d[first] * B1(uPrime[i])) +
  147. (
  148. ((Vector)d[last] * B2(uPrime[i])) +
  149. ((Vector)d[last] * B3(uPrime[i]))))));
  150.  
  151.  
  152. X[0] += V2Dot(A[i,0], tmp);
  153. X[1] += V2Dot(A[i,1], tmp);
  154. }
  155.  
  156. /* Compute the determinants of C and X */
  157. det_C0_C1 = C[0,0] * C[1,1] - C[1,0] * C[0,1];
  158. det_C0_X = C[0,0] * X[1] - C[1,0] * X[0];
  159. det_X_C1 = X[0] * C[1,1] - X[1] * C[0,1];
  160.  
  161. /* Finally, derive alpha values */
  162. alpha_l = (det_C0_C1 == 0) ? 0.0 : det_X_C1 / det_C0_C1;
  163. alpha_r = (det_C0_C1 == 0) ? 0.0 : det_C0_X / det_C0_C1;
  164.  
  165. /* If alpha negative, use the Wu/Barsky heuristic (see text) */
  166. /* (if alpha is 0, you get coincident control points that lead to
  167. * divide by zero in any subsequent NewtonRaphsonRootFind() call. */
  168. double segLength = (d[first] - d[last]).Length;
  169. double epsilon = 1.0e-6 * segLength;
  170. if (alpha_l < epsilon || alpha_r < epsilon)
  171. {
  172. /* fall back on standard (probably inaccurate) formula, and subdivide further if needed. */
  173. double dist = segLength / 3.0;
  174. bezCurve[0] = d[first];
  175. bezCurve[3] = d[last];
  176. bezCurve[1] = (tHat1 * dist) + bezCurve[0];
  177. bezCurve[2] = (tHat2 * dist) + bezCurve[3];
  178. return (bezCurve);
  179. }
  180.  
  181. /* First and last control points of the Bezier curve are */
  182. /* positioned exactly at the first and last data points */
  183. /* Control points 1 and 2 are positioned an alpha distance out */
  184. /* on the tangent vectors, left and right, respectively */
  185. bezCurve[0] = d[first];
  186. bezCurve[3] = d[last];
  187. bezCurve[1] = (tHat1 * alpha_l) + bezCurve[0];
  188. bezCurve[2] = (tHat2 * alpha_r) + bezCurve[3];
  189. return (bezCurve);
  190. }
  191.  
  192. /*
  193. * Reparameterize:
  194. * Given set of points and their parameterization, try to find
  195. * a better parameterization.
  196. *
  197. */
  198. static double[] Reparameterize(Point[] d,int first,int last,double[] u,Point[] bezCurve)
  199. {
  200. int nPts = last-first+1;
  201. int i;
  202. double[] uPrime = new double[nPts]; /* New parameter values */
  203.  
  204. for (i = first; i <= last; i++) {
  205. uPrime[i-first] = NewtonRaphsonRootFind(bezCurve, d[i], u[i-first]);
  206. }
  207. return uPrime;
  208. }
  209.  
  210.  
  211.  
  212. /*
  213. * NewtonRaphsonRootFind :
  214. * Use Newton-Raphson iteration to find better root.
  215. */
  216. static double NewtonRaphsonRootFind(Point[] Q,Point P,double u)
  217. {
  218. double numerator, denominator;
  219. Point[] Q1 = new Point[3], Q2 = new Point[2]; /* Q' and Q'' */
  220. Point Q_u, Q1_u, Q2_u; /*u evaluated at Q, Q', & Q'' */
  221. double uPrime; /* Improved u */
  222. int i;
  223.  
  224. /* Compute Q(u) */
  225. Q_u = BezierII(3, Q, u);
  226.  
  227. /* Generate control vertices for Q' */
  228. for (i = 0; i <= 2; i++) {
  229. Q1[i].X = (Q[i+1].X - Q[i].X) * 3.0;
  230. Q1[i].Y = (Q[i+1].Y - Q[i].Y) * 3.0;
  231. }
  232.  
  233. /* Generate control vertices for Q'' */
  234. for (i = 0; i <= 1; i++) {
  235. Q2[i].X = (Q1[i+1].X - Q1[i].X) * 2.0;
  236. Q2[i].Y = (Q1[i+1].Y - Q1[i].Y) * 2.0;
  237. }
  238.  
  239. /* Compute Q'(u) and Q''(u) */
  240. Q1_u = BezierII(2, Q1, u);
  241. Q2_u = BezierII(1, Q2, u);
  242.  
  243. /* Compute f(u)/f'(u) */
  244. numerator = (Q_u.X - P.X) * (Q1_u.X) + (Q_u.Y - P.Y) * (Q1_u.Y);
  245. denominator = (Q1_u.X) * (Q1_u.X) + (Q1_u.Y) * (Q1_u.Y) +
  246. (Q_u.X - P.X) * (Q2_u.X) + (Q_u.Y - P.Y) * (Q2_u.Y);
  247. if (denominator == 0.0f) return u;
  248.  
  249. /* u = u - f(u)/f'(u) */
  250. uPrime = u - (numerator/denominator);
  251. return (uPrime);
  252. }
  253.  
  254.  
  255.  
  256. /*
  257. * Bezier :
  258. * Evaluate a Bezier curve at a particular parameter value
  259. *
  260. */
  261. static Point BezierII(int degree,Point[] V,double t)
  262. {
  263. int i, j;
  264. Point Q; /* Point on curve at parameter t */
  265. Point[] Vtemp; /* Local copy of control points */
  266.  
  267. /* Copy array */
  268. Vtemp = new Point[degree+1];
  269. for (i = 0; i <= degree; i++) {
  270. Vtemp[i] = V[i];
  271. }
  272.  
  273. /* Triangle computation */
  274. for (i = 1; i <= degree; i++) {
  275. for (j = 0; j <= degree-i; j++) {
  276. Vtemp[j].X = (1.0 - t) * Vtemp[j].X + t * Vtemp[j+1].X;
  277. Vtemp[j].Y = (1.0 - t) * Vtemp[j].Y + t * Vtemp[j+1].Y;
  278. }
  279. }
  280.  
  281. Q = Vtemp[0];
  282. return Q;
  283. }
  284.  
  285.  
  286. /*
  287. * B0, B1, B2, B3 :
  288. * Bezier multipliers
  289. */
  290. static double B0(double u)
  291. {
  292. double tmp = 1.0 - u;
  293. return (tmp * tmp * tmp);
  294. }
  295.  
  296.  
  297. static double B1(double u)
  298. {
  299. double tmp = 1.0 - u;
  300. return (3 * u * (tmp * tmp));
  301. }
  302.  
  303. static double B2(double u)
  304. {
  305. double tmp = 1.0 - u;
  306. return (3 * u * u * tmp);
  307. }
  308.  
  309. static double B3(double u)
  310. {
  311. return (u * u * u);
  312. }
  313.  
  314. /*
  315. * ComputeLeftTangent, ComputeRightTangent, ComputeCenterTangent :
  316. *Approximate unit tangents at endpoints and "center" of digitized curve
  317. */
  318. static Vector ComputeLeftTangent(Point[] d,int end)
  319. {
  320. Vector tHat1;
  321. tHat1 = d[end+1]- d[end];
  322. tHat1.Normalize();
  323. return tHat1;
  324. }
  325.  
  326. static Vector ComputeRightTangent(Point[] d,int end)
  327. {
  328. Vector tHat2;
  329. tHat2 = d[end-1] - d[end];
  330. tHat2.Normalize();
  331. return tHat2;
  332. }
  333.  
  334. static Vector ComputeCenterTangent(Point[] d,int center)
  335. {
  336. Vector V1, V2, tHatCenter = new Vector();
  337.  
  338. V1 = d[center-1] - d[center];
  339. V2 = d[center] - d[center+1];
  340. tHatCenter.X = (V1.X + V2.X)/2.0;
  341. tHatCenter.Y = (V1.Y + V2.Y)/2.0;
  342. tHatCenter.Normalize();
  343. return tHatCenter;
  344. }
  345.  
  346.  
  347. /*
  348. * ChordLengthParameterize :
  349. * Assign parameter values to digitized points
  350. * using relative distances between points.
  351. */
  352. static double[] ChordLengthParameterize(Point[] d,int first,int last)
  353. {
  354. int i;
  355. double[] u = new double[last-first+1]; /* Parameterization */
  356.  
  357. u[0] = 0.0;
  358. for (i = first+1; i <= last; i++) {
  359. u[i-first] = u[i-first-1] + (d[i-1] - d[i]).Length;
  360. }
  361.  
  362. for (i = first + 1; i <= last; i++) {
  363. u[i-first] = u[i-first] / u[last-first];
  364. }
  365.  
  366. return u;
  367. }
  368.  
  369.  
  370.  
  371.  
  372. /*
  373. * ComputeMaxError :
  374. * Find the maximum squared distance of digitized points
  375. * to fitted curve.
  376. */
  377. static double ComputeMaxError(Point[] d,int first,int last,Point[] bezCurve,double[] u,out int splitPoint)
  378. {
  379. int i;
  380. double maxDist; /* Maximum error */
  381. double dist; /* Current error */
  382. Point P; /* Point on curve */
  383. Vector v; /* Vector from point to curve */
  384.  
  385. splitPoint = (last - first + 1)/2;
  386. maxDist = 0.0;
  387. for (i = first + 1; i < last; i++) {
  388. P = BezierII(3, bezCurve, u[i-first]);
  389. v = P - d[i];
  390. dist = v.LengthSquared;
  391. if (dist >= maxDist) {
  392. maxDist = dist;
  393. splitPoint = i;
  394. }
  395. }
  396. return maxDist;
  397. }
  398.  
  399. private static double V2Dot(Vector a,Vector b)
  400. {
  401. return((a.X*b.X)+(a.Y*b.Y));
  402. }
  403.  
  404. }
  405.  
  406. public class Vector
  407. {
  408. public double X { get; set; }
  409. public double Y { get; set; }
  410. public Vector(double x=0, double y=0)
  411. {
  412. X = x;
  413. Y = y;
  414. }
  415.  
  416. public static implicit operator Vector(Point b)
  417. {
  418. return new Vector(b.X, b.Y);
  419. }
  420.  
  421. public static Point operator *(Vector left, double right)
  422. {
  423. return new Point(left.X * right, left.Y * right);
  424. }
  425. public static Vector operator -(Vector left, Point right)
  426. {
  427. return new Vector(left.X - right.X, left.Y - right.Y);
  428. }
  429.  
  430. internal void Negate()
  431. {
  432. X = -X;
  433. Y = -Y;
  434. }
  435.  
  436. internal void Normalize()
  437. {
  438. double factor = 1.0 / Math.Sqrt(LengthSquared);
  439. X *= factor;
  440. Y *= factor;
  441. }
  442.  
  443. public double LengthSquared { get { return X * X + Y * Y; } }
  444. }
  445.  
  446. public static double Length(Point a, Point b)
  447. {
  448. double x = a.X-b.X;
  449. double y = a.Y-b.Y;
  450. return Math.Sqrt(x*x+y*y);
  451. }
  452.  
  453. public static Point Add(Point a, Point b)
  454. {
  455. return new Point(a.X + b.X, a.Y + b.Y);
  456. }
  457.  
  458. public static Point Subtract(Point a, Point b)
  459. {
  460. return new Point(a.X - b.X, a.Y - b.Y);
  461. }
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement
Public Pastes
Advertisement
We use cookies for various purposes including analytics. By continuing to use Pastebin, you agree to our use of cookies as described in the Cookies Policy.  OK, I Understand
Not a member of Pastebin yet?
Sign Up, it unlocks many cool features!