+ samples.clear ();
+ interpolate (_points, points_per_segment, CatmullRomCentripetal, false, samples);
+ n_samples = samples.size();
+}
+
+/* Cartmull-Rom code from http://stackoverflow.com/questions/9489736/catmull-rom-curve-with-no-cusps-and-no-self-intersections/19283471#19283471
+ *
+ * Thanks to Ted for his Java version, which I translated into Ardour-idiomatic
+ * C++ here.
+ */
+
+/**
+ * Calculate the same values but introduces the ability to "parameterize" the t
+ * values used in the calculation. This is based on Figure 3 from
+ * http://www.cemyuksel.com/research/catmullrom_param/catmullrom.pdf
+ *
+ * @param p An array of double values of length 4, where interpolation
+ * occurs from p1 to p2.
+ * @param time An array of time measures of length 4, corresponding to each
+ * p value.
+ * @param t the actual interpolation ratio from 0 to 1 representing the
+ * position between p1 and p2 to interpolate the value.
+ */
+static double
+__interpolate (double p[4], double time[4], double t)
+{
+ const double L01 = p[0] * (time[1] - t) / (time[1] - time[0]) + p[1] * (t - time[0]) / (time[1] - time[0]);
+ const double L12 = p[1] * (time[2] - t) / (time[2] - time[1]) + p[2] * (t - time[1]) / (time[2] - time[1]);
+ const double L23 = p[2] * (time[3] - t) / (time[3] - time[2]) + p[3] * (t - time[2]) / (time[3] - time[2]);
+ const double L012 = L01 * (time[2] - t) / (time[2] - time[0]) + L12 * (t - time[0]) / (time[2] - time[0]);
+ const double L123 = L12 * (time[3] - t) / (time[3] - time[1]) + L23 * (t - time[1]) / (time[3] - time[1]);
+ const double C12 = L012 * (time[2] - t) / (time[2] - time[1]) + L123 * (t - time[1]) / (time[2] - time[1]);
+ return C12;
+}
+
+/**
+ * Given a list of control points, this will create a list of points_per_segment
+ * points spaced uniformly along the resulting Catmull-Rom curve.
+ *
+ * @param points The list of control points, leading and ending with a
+ * coordinate that is only used for controling the spline and is not visualized.
+ * @param index The index of control point p0, where p0, p1, p2, and p3 are
+ * used in order to create a curve between p1 and p2.
+ * @param points_per_segment The total number of uniformly spaced interpolated
+ * points to calculate for each segment. The larger this number, the
+ * smoother the resulting curve.
+ * @param curve_type Clarifies whether the curve should use uniform, chordal
+ * or centripetal curve types. Uniform can produce loops, chordal can
+ * produce large distortions from the original lines, and centripetal is an
+ * optimal balance without spaces.
+ * @return the list of coordinates that define the CatmullRom curve
+ * between the points defined by index+1 and index+2.
+ */
+static void
+_interpolate (const Points& points, Points::size_type index, int points_per_segment, Curve::SplineType curve_type, Points& results)
+{
+ double x[4];
+ double y[4];
+ double time[4];
+
+ for (int i = 0; i < 4; i++) {
+ x[i] = points[index + i].x;
+ y[i] = points[index + i].y;
+ time[i] = i;
+ }
+
+ double tstart = 1;
+ double tend = 2;
+
+ if (curve_type != Curve::CatmullRomUniform) {
+ double total = 0;
+ for (int i = 1; i < 4; i++) {
+ double dx = x[i] - x[i - 1];
+ double dy = y[i] - y[i - 1];
+ if (curve_type == Curve::CatmullRomCentripetal) {
+ total += pow (dx * dx + dy * dy, .25);
+ } else {
+ total += pow (dx * dx + dy * dy, .5);
+ }
+ time[i] = total;
+ }
+ tstart = time[1];
+ tend = time[2];
+ }
+
+ int segments = points_per_segment - 1;
+ results.push_back (points[index + 1]);
+
+ for (int i = 1; i < segments; i++) {
+ double xi = __interpolate (x, time, tstart + (i * (tend - tstart)) / segments);
+ double yi = __interpolate (y, time, tstart + (i * (tend - tstart)) / segments);
+ results.push_back (Duple (xi, yi));
+ }
+
+ results.push_back (points[index + 2]);