diff options
author | Paul Davis <paul@linuxaudiosystems.com> | 2014-03-05 11:37:13 -0500 |
---|---|---|
committer | Paul Davis <paul@linuxaudiosystems.com> | 2014-03-05 11:38:30 -0500 |
commit | 58a30da03d31f5a30faddd76efe4e76b242f6687 (patch) | |
tree | 9ff5470017eae8497e033bc59e825dbe21c91aff /libs/canvas/curve.cc | |
parent | f3300ec03ce82628aa4eb57565b1c05acb2915e9 (diff) |
use a centripetal catmull-rom curve to smooth ArdourCanvas::Curve
See http://en.wikipedia.org/wiki/Centripetal_Catmull-Rom to understand the benefits of this.
Diffstat (limited to 'libs/canvas/curve.cc')
-rw-r--r-- | libs/canvas/curve.cc | 370 |
1 files changed, 195 insertions, 175 deletions
diff --git a/libs/canvas/curve.cc b/libs/canvas/curve.cc index 54d36bbc9b..128aea5462 100644 --- a/libs/canvas/curve.cc +++ b/libs/canvas/curve.cc @@ -17,6 +17,7 @@ */ +#include <cmath> #include <exception> #include <algorithm> @@ -31,24 +32,9 @@ Curve::Curve (Group* parent) , PolyItem (parent) , Fill (parent) , n_samples (0) - , n_segments (512) + , points_per_segment (16) + , curve_type (CatmullRomCentripetal) { - set_n_samples (256); -} - -/** Set the number of points to compute when we smooth the data points into a - * curve. - */ -void -Curve::set_n_samples (Points::size_type n) -{ - /* this only changes our appearance rather than the bounding box, so we - just need to schedule a redraw rather than notify the parent of any - changes - */ - n_samples = n; - samples.assign (n_samples, Duple (0.0, 0.0)); - interpolate (); } /** When rendering the curve, we will always draw a fixed number of straight @@ -57,13 +43,14 @@ Curve::set_n_samples (Points::size_type n) * render. */ void -Curve::set_n_segments (uint32_t n) +Curve::set_points_per_segment (uint32_t n) { /* this only changes our appearance rather than the bounding box, so we just need to schedule a redraw rather than notify the parent of any changes */ - n_segments = n; + points_per_segment = n; + interpolate (); redraw (); } @@ -85,171 +72,204 @@ Curve::set (Points const& p) void Curve::interpolate () { - Points::size_type npoints = _points.size (); - - if (npoints < 3) { - return; - } - - Duple p; - double boundary; - - const double xfront = _points.front().x; - const double xextent = _points.back().x - xfront; - - /* initialize boundary curve points */ - - p = _points.front(); - boundary = round (((p.x - xfront)/xextent) * (n_samples - 1)); - - for (Points::size_type i = 0; i < boundary; ++i) { - samples[i] = Duple (i, p.y); - } - - p = _points.back(); - boundary = round (((p.x - xfront)/xextent) * (n_samples - 1)); - - for (Points::size_type i = boundary; i < n_samples; ++i) { - samples[i] = Duple (i, p.y); - } - - for (int i = 0; i < (int) npoints - 1; ++i) { - - Points::size_type p1, p2, p3, p4; - - p1 = max (i - 1, 0); - p2 = i; - p3 = i + 1; - p4 = min (i + 2, (int) npoints - 1); + samples.clear (); + interpolate (_points, points_per_segment, CatmullRomCentripetal, false, samples); + n_samples = samples.size(); +} - smooth (p1, p2, p3, p4, xfront, xextent); - } - - /* make sure that actual data points are used with their exact values */ +/* 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. + */ - for (Points::const_iterator p = _points.begin(); p != _points.end(); ++p) { - uint32_t idx = (((*p).x - xfront)/xextent) * (n_samples - 1); - samples[idx].y = (*p).y; - } +/** + * 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]); } -/* - * This function calculates the curve values between the control points - * p2 and p3, taking the potentially existing neighbors p1 and p4 into - * account. - * - * This function uses a cubic bezier curve for the individual segments and - * calculates the necessary intermediate control points depending on the - * neighbor curve control points. +/** + * This method will calculate the Catmull-Rom interpolation curve, returning + * it as a list of Coord coordinate objects. This method in particular + * adds the first and last control points which are not visible, but required + * for calculating the spline. * + * @param coordinates The list of original straight line points to calculate + * an interpolation from. + * @param points_per_segment The integer number of equally spaced points to + * return along each curve. The actual distance between each + * point will depend on the spacing between the control points. + * @return The list of interpolated coordinates. + * @param curve_type Chordal (stiff), Uniform(floppy), or Centripetal(medium) + * @throws gov.ca.water.shapelite.analysis.CatmullRomException if + * points_per_segment is less than 2. */ void -Curve::smooth (Points::size_type p1, Points::size_type p2, Points::size_type p3, Points::size_type p4, - double xfront, double xextent) +Curve::interpolate (const Points& coordinates, uint32_t points_per_segment, SplineType curve_type, bool closed, Points& results) { - int i; - double x0, x3; - double y0, y1, y2, y3; - double dx, dy; - double slope; - - /* the outer control points for the bezier curve. */ - - x0 = _points[p2].x; - y0 = _points[p2].y; - x3 = _points[p3].x; - y3 = _points[p3].y; - - /* - * the x values of the inner control points are fixed at - * x1 = 2/3*x0 + 1/3*x3 and x2 = 1/3*x0 + 2/3*x3 - * this ensures that the x values increase linearily with the - * parameter t and enables us to skip the calculation of the x - * values altogehter - just calculate y(t) evenly spaced. - */ - - dx = x3 - x0; - dy = y3 - y0; - - if (dx <= 0) { - /* error? */ - return; - } - - if (p1 == p2 && p3 == p4) { - /* No information about the neighbors, - * calculate y1 and y2 to get a straight line - */ - y1 = y0 + dy / 3.0; - y2 = y0 + dy * 2.0 / 3.0; - - } else if (p1 == p2 && p3 != p4) { - - /* only the right neighbor is available. Make the tangent at the - * right endpoint parallel to the line between the left endpoint - * and the right neighbor. Then point the tangent at the left towards - * the control handle of the right tangent, to ensure that the curve - * does not have an inflection point. - */ - slope = (_points[p4].y - y0) / (_points[p4].x - x0); - - y2 = y3 - slope * dx / 3.0; - y1 = y0 + (y2 - y0) / 2.0; - - } else if (p1 != p2 && p3 == p4) { - - /* see previous case */ - slope = (y3 - _points[p1].y) / (x3 - _points[p1].x); - - y1 = y0 + slope * dx / 3.0; - y2 = y3 + (y1 - y3) / 2.0; - - - } else /* (p1 != p2 && p3 != p4) */ { - - /* Both neighbors are available. Make the tangents at the endpoints - * parallel to the line between the opposite endpoint and the adjacent - * neighbor. - */ - - slope = (y3 - _points[p1].y) / (x3 - _points[p1].x); - - y1 = y0 + slope * dx / 3.0; - - slope = (_points[p4].y - y0) / (_points[p4].x - x0); - - y2 = y3 - slope * dx / 3.0; - } - - /* - * finally calculate the y(t) values for the given bezier values. We can - * use homogenously distributed values for t, since x(t) increases linearily. - */ - - dx = dx / xextent; - - int limit = round (dx * (n_samples - 1)); - const int idx_offset = round (((x0 - xfront)/xextent) * (n_samples - 1)); - - for (i = 0; i <= limit; i++) { - double y, t; - Points::size_type index; - - t = i / dx / (n_samples - 1); - - y = y0 * (1-t) * (1-t) * (1-t) + - 3 * y1 * (1-t) * (1-t) * t + - 3 * y2 * (1-t) * t * t + - y3 * t * t * t; - - index = i + idx_offset; - - if (index < n_samples) { - Duple d (i, max (y, 0.0)); - samples[index] = d; - } - } + if (points_per_segment < 2) { + return; + } + + // Cannot interpolate curves given only two points. Two points + // is best represented as a simple line segment. + if (coordinates.size() < 3) { + results = coordinates; + return; + } + + // Copy the incoming coordinates. We need to modify it during interpolation + Points vertices = coordinates; + + // Test whether the shape is open or closed by checking to see if + // the first point intersects with the last point. M and Z are ignored. + if (closed) { + // Use the second and second from last points as control points. + // get the second point. + Duple p2 = vertices[1]; + // get the point before the last point + Duple pn1 = vertices[vertices.size() - 2]; + + // insert the second from the last point as the first point in the list + // because when the shape is closed it keeps wrapping around to + // the second point. + vertices.insert(vertices.begin(), pn1); + // add the second point to the end. + vertices.push_back(p2); + } else { + // The shape is open, so use control points that simply extend + // the first and last segments + + // Get the change in x and y between the first and second coordinates. + double dx = vertices[1].x - vertices[0].x; + double dy = vertices[1].y - vertices[0].y; + + // Then using the change, extrapolate backwards to find a control point. + double x1 = vertices[0].x - dx; + double y1 = vertices[0].y - dy; + + // Actaully create the start point from the extrapolated values. + Duple start (x1, y1); + + // Repeat for the end control point. + int n = vertices.size() - 1; + dx = vertices[n].x - vertices[n - 1].x; + dy = vertices[n].y - vertices[n - 1].y; + double xn = vertices[n].x + dx; + double yn = vertices[n].y + dy; + Duple end (xn, yn); + + // insert the start control point at the start of the vertices list. + vertices.insert (vertices.begin(), start); + + // append the end control ponit to the end of the vertices list. + vertices.push_back (end); + } + + // When looping, remember that each cycle requires 4 points, starting + // with i and ending with i+3. So we don't loop through all the points. + + for (Points::size_type i = 0; i < vertices.size() - 3; i++) { + + // Actually calculate the Catmull-Rom curve for one segment. + Points r; + + _interpolate (vertices, i, points_per_segment, curve_type, r); + + // Since the middle points are added twice, once for each bordering + // segment, we only add the 0 index result point for the first + // segment. Otherwise we will have duplicate points. + + if (results.size() > 0) { + r.erase (r.begin()); + } + + // Add the coordinates for the segment to the result list. + + results.insert (results.end(), r.begin(), r.end()); + } } /** Given a fractional position within the x-axis range of the |