2 files changed, 91 insertions, 2 deletions

diff --git a/src/live_effects/lpe-powerstroke-interpolators.h b/src/live_effects/lpe-powerstroke-interpolators.h
index 6f5b75af8..ee4951882 100644
--- a/src/live_effects/lpe-powerstroke-interpolators.h
+++ b/src/live_effects/lpe-powerstroke-interpolators.h
@@ -27,7 +27,8 @@ enum InterpolatorType {
   INTERP_LINEAR,
   INTERP_CUBICBEZIER,
   INTERP_CUBICBEZIER_JOHAN,
-  INTERP_SPIRO
+  INTERP_SPIRO,
+  INTERP_CENTRIPETAL_CATMULLROM
 };
 
 class Interpolator {
@@ -167,6 +168,91 @@ private:
     SpiroInterpolator& operator=(const SpiroInterpolator&);
 };
 
+class CentripetalCatmullRomInterpolator : public Interpolator {
+// the code in this class can certainly be optimized and made more terse, feel free    
+public:
+    CentripetalCatmullRomInterpolator() {};
+    virtual ~CentripetalCatmullRomInterpolator() {};
+
+    virtual Path interpolateToPath(std::vector<Point> const &points) const {
+        unsigned int n_points = points.size();
+
+        Geom::Path fit(points.front());
+
+        if (n_points < 2) return fit;
+        if (n_points < 3) {
+            // if only 2 points, return linear segment
+            fit.appendNew<Geom::LineSegment>(points[1]);
+            return fit;
+        }
+
+        // return n_points-1 cubic segments
+
+        // duplicate first point
+        fit.append(calc_bezier(points[0],points[0],points[1],points[2]));
+
+        for (std::size_t i = 0; i < n_points-2; ++i) {
+            Point p0 = points[i];
+            Point p1 = points[i+1];
+            Point p2 = points[i+2];
+            Point p3 = (i < n_points-3) ? points[i+3] : points[i+2];
+
+            fit.append(calc_bezier(p0, p1, p2, p3));
+            // this is quite wasteful: the distances between the same points are repeatedly calculated by calc_bezier
+        }
+
+        return fit;
+    };
+
+private:
+    CubicBezier calc_bezier(Point p0, Point p1, Point p2, Point p3) const {
+        // create interpolating bezier between p1 and p2
+
+        // Part of the code comes from StackOverflow user eriatarka84
+        // http://stackoverflow.com/a/23980479/2929337
+
+        // calculate time coords (deltas) of points
+        // the factor 0.25 can be generalized for other Catmull-Rom interpolation types 
+        // see alpha in Yuksel et al. "On the Parameterization of Catmull-Rom Curves", 
+        // --> http://www.cemyuksel.com/research/catmullrom_param/catmullrom.pdf
+        double dt0 = powf(distanceSq(p0, p1), 0.25);
+        double dt1 = powf(distanceSq(p1, p2), 0.25);
+        double dt2 = powf(distanceSq(p2, p3), 0.25);
+
+
+        // safety check for repeated points
+        double eps = Geom::EPSILON;
+        if (dt1 < eps)
+            dt1 = 1.0;
+        if (dt0 < eps)
+            dt0 = dt1;
+        if (dt2 < eps)
+            dt2 = dt1;
+
+        // compute tangents when parameterized in [t1,t2]
+        Point tan1 = (p1 - p0) / dt0 - (p2 - p0) / (dt0 + dt1) + (p2 - p1) / dt1;
+        Point tan2 = (p2 - p1) / dt1 - (p3 - p1) / (dt1 + dt2) + (p3 - p2) / dt2;
+        // rescale tangents for parametrization in [0,1]
+        tan1 *= dt1;
+        tan2 *= dt1;
+
+        // create bezier from tangents (this is already in 2geom somewhere, or should be moved to it)
+        // the tangent of a bezier curve is: B'(t) = 3(1-t)^2 (b1 - b0) + 6(1-t)t(b2-b1) + 3t^2(b3-b2)
+        // So we have to make sure that B'(0) = tan1  and  B'(1) = tan2, and we already know that b0=p1 and b3=p2
+        // tan1 = B'(0) = 3 (b1 - p1)  --> p1 + (tan1)/3 = b1
+        // tan2 = B'(1) = 3 (p2 - b2)  --> p2 - (tan2)/3 = b2
+
+        Point b0 = p1;
+        Point b1 = p1 + tan1 / 3;
+        Point b2 = p2 - tan2 / 3;
+        Point b3 = p2;
+
+        return CubicBezier(b0, b1, b2, b3);
+    }
+
+    CentripetalCatmullRomInterpolator(const CentripetalCatmullRomInterpolator&);
+    CentripetalCatmullRomInterpolator& operator=(const CentripetalCatmullRomInterpolator&);
+};
 
 Interpolator*
 Interpolator::create(InterpolatorType type) {
@@ -179,6 +265,8 @@ Interpolator::create(InterpolatorType type) {
         return new Geom::Interpolate::CubicBezierJohan();
       case INTERP_SPIRO:
         return new Geom::Interpolate::SpiroInterpolator();
+      case INTERP_CENTRIPETAL_CATMULLROM:
+        return new Geom::Interpolate::CentripetalCatmullRomInterpolator();
       default:
         return new Geom::Interpolate::Linear();
     }
diff --git a/src/live_effects/lpe-powerstroke.cpp b/src/live_effects/lpe-powerstroke.cpp
index 90b01aaa4..eed564b20 100644
--- a/src/live_effects/lpe-powerstroke.cpp
+++ b/src/live_effects/lpe-powerstroke.cpp
@@ -188,7 +188,8 @@ static const Util::EnumData<unsigned> InterpolatorTypeData[] = {
     {Geom::Interpolate::INTERP_LINEAR          , N_("Linear"), "Linear"},
     {Geom::Interpolate::INTERP_CUBICBEZIER          , N_("CubicBezierFit"), "CubicBezierFit"},
     {Geom::Interpolate::INTERP_CUBICBEZIER_JOHAN     , N_("CubicBezierJohan"), "CubicBezierJohan"},
-    {Geom::Interpolate::INTERP_SPIRO  , N_("SpiroInterpolator"), "SpiroInterpolator"}
+    {Geom::Interpolate::INTERP_SPIRO  , N_("SpiroInterpolator"), "SpiroInterpolator"},
+    {Geom::Interpolate::INTERP_CENTRIPETAL_CATMULLROM  , N_("Centripetal Catmull-Rom"), "CentripetalCatmullRom"}
 };
 static const Util::EnumDataConverter<unsigned> InterpolatorTypeConverter(InterpolatorTypeData, sizeof(InterpolatorTypeData)/sizeof(*InterpolatorTypeData));