fork(37) download
  1. #include <iostream>
  2. #include <cmath>
  3.  
  4. using namespace std;
  5.  
  6. struct CubicPoly
  7. {
  8. float c0, c1, c2, c3;
  9.  
  10. float eval(float t)
  11. {
  12. float t2 = t*t;
  13. float t3 = t2 * t;
  14. return c0 + c1*t + c2*t2 + c3*t3;
  15. }
  16. };
  17.  
  18. /*
  19.  * Compute coefficients for a cubic polynomial
  20.  * p(s) = c0 + c1*s + c2*s^2 + c3*s^3
  21.  * such that
  22.  * p(0) = x0, p(1) = x1
  23.  * and
  24.  * p'(0) = t0, p'(1) = t1.
  25.  */
  26. void InitCubicPoly(float x0, float x1, float t0, float t1, CubicPoly &p)
  27. {
  28. p.c0 = x0;
  29. p.c1 = t0;
  30. p.c2 = -3*x0 + 3*x1 - 2*t0 - t1;
  31. p.c3 = 2*x0 - 2*x1 + t0 + t1;
  32. }
  33.  
  34. // standard Catmull-Rom spline: interpolate between x1 and x2 with previous/following points x1/x4
  35. // (we don't need this here, but it's for illustration)
  36. void InitCatmullRom(float x0, float x1, float x2, float x3, CubicPoly &p)
  37. {
  38. // Catmull-Rom with tension 0.5
  39. InitCubicPoly(x1, x2, 0.5f*(x2-x0), 0.5f*(x3-x1), p);
  40. }
  41.  
  42. // compute coefficients for a nonuniform Catmull-Rom spline
  43. void InitNonuniformCatmullRom(float x0, float x1, float x2, float x3, float dt0, float dt1, float dt2, CubicPoly &p)
  44. {
  45. // compute tangents when parameterized in [t1,t2]
  46. float t1 = (x1 - x0) / dt0 - (x2 - x0) / (dt0 + dt1) + (x2 - x1) / dt1;
  47. float t2 = (x2 - x1) / dt1 - (x3 - x1) / (dt1 + dt2) + (x3 - x2) / dt2;
  48.  
  49. // rescale tangents for parametrization in [0,1]
  50. t1 *= dt1;
  51. t2 *= dt1;
  52.  
  53. InitCubicPoly(x1, x2, t1, t2, p);
  54. }
  55.  
  56. struct Vec2D
  57. {
  58. Vec2D(float _x, float _y) : x(_x), y(_y) {}
  59. float x, y;
  60. };
  61.  
  62. float VecDistSquared(const Vec2D& p, const Vec2D& q)
  63. {
  64. float dx = q.x - p.x;
  65. float dy = q.y - p.y;
  66. return dx*dx + dy*dy;
  67. }
  68.  
  69. void InitCentripetalCR(const Vec2D& p0, const Vec2D& p1, const Vec2D& p2, const Vec2D& p3,
  70. CubicPoly &px, CubicPoly &py)
  71. {
  72. float dt0 = powf(VecDistSquared(p0, p1), 0.25f);
  73. float dt1 = powf(VecDistSquared(p1, p2), 0.25f);
  74. float dt2 = powf(VecDistSquared(p2, p3), 0.25f);
  75.  
  76. // safety check for repeated points
  77. if (dt1 < 1e-4f) dt1 = 1.0f;
  78. if (dt0 < 1e-4f) dt0 = dt1;
  79. if (dt2 < 1e-4f) dt2 = dt1;
  80.  
  81. InitNonuniformCatmullRom(p0.x, p1.x, p2.x, p3.x, dt0, dt1, dt2, px);
  82. InitNonuniformCatmullRom(p0.y, p1.y, p2.y, p3.y, dt0, dt1, dt2, py);
  83. }
  84.  
  85.  
  86. int main()
  87. {
  88. Vec2D p0(0,0), p1(1,1), p2(1.1,1), p3(2,0);
  89. CubicPoly px, py;
  90. InitCentripetalCR(p0, p1, p2, p3, px, py);
  91. for (int i = 0; i <= 10; ++i)
  92. cout << px.eval(0.1f*i) << " " << py.eval(0.1f*i) << endl;
  93. }
Success #stdin #stdout 0s 3296KB
stdin
Standard input is empty
stdout
1 1
1.01254 1.00505
1.02346 1.00902
1.03316 1.01189
1.04203 1.01365
1.05046 1.01428
1.05886 1.01377
1.06762 1.0121
1.07713 1.00926
1.08779 1.00523
1.1 1