X-Git-Url: https://main.carlh.net/gitweb/?a=blobdiff_plain;f=libs%2Fcanvas%2Fcurve.cc;h=128aea5462c6f90132d41e9d03f4a334405e743e;hb=621887cfaac4bef4b6849037c0d72f7e7b66fa03;hp=de988ee4b255bba524f85a28951ca585b07d83c0;hpb=2a7ed69c28c5c4606ff13b3605b9bc9c3eba607d;p=ardour.git diff --git a/libs/canvas/curve.cc b/libs/canvas/curve.cc index de988ee4b2..128aea5462 100644 --- a/libs/canvas/curve.cc +++ b/libs/canvas/curve.cc @@ -17,6 +17,7 @@ */ +#include #include #include @@ -29,8 +30,28 @@ using std::max; Curve::Curve (Group* parent) : Item (parent) , PolyItem (parent) + , Fill (parent) + , n_samples (0) + , points_per_segment (16) + , curve_type (CatmullRomCentripetal) { +} +/** When rendering the curve, we will always draw a fixed number of straight + * line segments to span the x-axis extent of the curve. More segments: + * smoother visual rendering. Less rendering: closer to a visibily poly-line + * render. + */ +void +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 + */ + points_per_segment = n; + interpolate (); + redraw (); } void @@ -38,176 +59,359 @@ Curve::compute_bounding_box () const { PolyItem::compute_bounding_box (); - if (_bounding_box) { - - bool have1 = false; - bool have2 = false; - - Rect bbox1; - Rect bbox2; - - for (Points::const_iterator i = first_control_points.begin(); i != first_control_points.end(); ++i) { - if (have1) { - bbox1.x0 = min (bbox1.x0, i->x); - bbox1.y0 = min (bbox1.y0, i->y); - bbox1.x1 = max (bbox1.x1, i->x); - bbox1.y1 = max (bbox1.y1, i->y); - } else { - bbox1.x0 = bbox1.x1 = i->x; - bbox1.y0 = bbox1.y1 = i->y; - have1 = true; - } - } - - for (Points::const_iterator i = second_control_points.begin(); i != second_control_points.end(); ++i) { - if (have2) { - bbox2.x0 = min (bbox2.x0, i->x); - bbox2.y0 = min (bbox2.y0, i->y); - bbox2.x1 = max (bbox2.x1, i->x); - bbox2.y1 = max (bbox2.y1, i->y); - } else { - bbox2.x0 = bbox2.x1 = i->x; - bbox2.y0 = bbox2.y1 = i->y; - have2 = true; - } - } - - Rect u = bbox1.extend (bbox2); - _bounding_box = u.extend (_bounding_box.get()); - } - - _bounding_box_dirty = false; + /* possibly add extents of any point indicators here if we ever do that */ } void Curve::set (Points const& p) { PolyItem::set (p); - - first_control_points.clear (); - second_control_points.clear (); - - compute_control_points (_points, first_control_points, second_control_points); + interpolate (); } void -Curve::render (Rect const & area, Cairo::RefPtr context) const +Curve::interpolate () { - if (_outline) { - setup_outline_context (context); - render_path (area, context); - context->stroke (); - } + 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]); } -void -Curve::render_path (Rect const & area, Cairo::RefPtr context) const +/** + * 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::interpolate (const Points& coordinates, uint32_t points_per_segment, SplineType curve_type, bool closed, Points& results) { - PolyItem::render_curve (area, context, first_control_points, second_control_points); + 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()); + } } -void -Curve::compute_control_points (Points const& knots, - Points& firstControlPoints, - Points& secondControlPoints) +/** Given a fractional position within the x-axis range of the + * curve, return the corresponding y-axis value + */ + +double +Curve::map_value (double x) const { - Points::size_type n = knots.size() - 1; - - if (n < 1) { - return; - } - - if (n == 1) { - /* Special case: Bezier curve should be a straight line. */ - - Duple d; + if (x > 0.0 && x < 1.0) { + + double f; + Points::size_type index; - d.x = (2.0 * knots[0].x + knots[1].x) / 3; - d.y = (2.0 * knots[0].y + knots[1].y) / 3; - firstControlPoints.push_back (d); + /* linearly interpolate between two of our smoothed "samples" + */ - d.x = 2.0 * firstControlPoints[0].x - knots[0].x; - d.y = 2.0 * firstControlPoints[0].y - knots[0].y; - secondControlPoints.push_back (d); + x = x * (n_samples - 1); + index = (Points::size_type) x; // XXX: should we explicitly use floor()? + f = x - index; + + return (1.0 - f) * samples[index].y + f * samples[index+1].y; - return; - } - - // Calculate first Bezier control points - // Right hand side vector - - std::vector rhs; - - rhs.assign (n, 0); - - // Set right hand side X values - - for (Points::size_type i = 1; i < n - 1; ++i) { - rhs[i] = 4 * knots[i].x + 2 * knots[i + 1].x; + } else if (x >= 1.0) { + return samples.back().y; + } else { + return samples.front().y; } - rhs[0] = knots[0].x + 2 * knots[1].x; - rhs[n - 1] = (8 * knots[n - 1].x + knots[n].x) / 2.0; - - // Get first control points X-values - double* x = solve (rhs); +} - // Set right hand side Y values - for (Points::size_type i = 1; i < n - 1; ++i) { - rhs[i] = 4 * knots[i].y + 2 * knots[i + 1].y; +void +Curve::render (Rect const & area, Cairo::RefPtr context) const +{ + if (!_outline || _points.size() < 2 || !_bounding_box) { + return; } - rhs[0] = knots[0].y + 2 * knots[1].y; - rhs[n - 1] = (8 * knots[n - 1].y + knots[n].y) / 2.0; - - // Get first control points Y-values - double* y = solve (rhs); + + Rect self = item_to_window (_bounding_box.get()); + boost::optional d = self.intersection (area); + assert (d); + Rect draw = d.get (); + + /* Our approach is to always draw n_segments across our total size. + * + * This is very inefficient if we are asked to only draw a small + * section of the curve. For now we rely on cairo clipping to help + * with this. + */ - for (Points::size_type i = 0; i < n; ++i) { - - firstControlPoints.push_back (Duple (x[i], y[i])); + + setup_outline_context (context); + + if (_points.size() == 2) { + + /* straight line */ + + Duple window_space; + + window_space = item_to_window (_points.front()); + context->move_to (window_space.x, window_space.y); + window_space = item_to_window (_points.back()); + context->line_to (window_space.x, window_space.y); + + context->stroke (); + + } else { + + /* curve of at least 3 points */ + + /* x-axis limits of the curve, in window space coordinates */ + + Duple w1 = item_to_window (Duple (_points.front().x, 0.0)); + Duple w2 = item_to_window (Duple (_points.back().x, 0.0)); + + /* clamp actual draw to area bound by points, rather than our bounding box which is slightly different */ + + context->save (); + context->rectangle (draw.x0, draw.y0, draw.width(), draw.height()); + context->clip (); + + /* expand drawing area by several pixels on each side to avoid cairo stroking effects at the boundary. + they will still occur, but cairo's clipping will hide them. + */ + + draw = draw.expand (4.0); + + /* now clip it to the actual points in the curve */ - if (i < n - 1) { - secondControlPoints.push_back (Duple (2 * knots [i + 1].x - x[i + 1], - 2 * knots[i + 1].y - y[i + 1])); - } else { - secondControlPoints.push_back (Duple ((knots [n].x + x[n - 1]) / 2, - (knots[n].y + y[n - 1]) / 2)); + if (draw.x0 < w1.x) { + draw.x0 = w1.x; } - } - - delete [] x; - delete [] y; -} -/** Solves a tridiagonal system for one of coordinates (x or y) - * of first Bezier control points. - */ + if (draw.x1 >= w2.x) { + draw.x1 = w2.x; + } -double* -Curve::solve (std::vector const & rhs) -{ - std::vector::size_type n = rhs.size(); - double* x = new double[n]; // Solution vector. - double* tmp = new double[n]; // Temp workspace. - - double b = 2.0; + /* full width of the curve */ + const double xextent = _points.back().x - _points.front().x; + /* Determine where the first drawn point will be */ + Duple item_space = window_to_item (Duple (draw.x0, 0)); /* y value is irrelevant */ + /* determine the fractional offset of this location into the overall extent of the curve */ + const double xfract_offset = (item_space.x - _points.front().x)/xextent; + const uint32_t pixels = draw.width (); + Duple window_space; - x[0] = rhs[0] / b; + /* draw the first point */ - for (std::vector::size_type i = 1; i < n; i++) { - // Decomposition and forward substitution. - tmp[i] = 1 / b; - b = (i < n - 1 ? 4.0 : 3.5) - tmp[i]; - x[i] = (rhs[i] - x[i - 1]) / b; - } - - for (std::vector::size_type i = 1; i < n; i++) { - // Backsubstitution - x[n - i - 1] -= tmp[n - i] * x[n - i]; + for (uint32_t pixel = 0; pixel < pixels; ++pixel) { + + /* fractional distance into the total horizontal extent of the curve */ + double xfract = xfract_offset + (pixel / xextent); + /* compute vertical coordinate (item-space) at that location */ + double y = map_value (xfract); + + /* convert to window space for drawing */ + window_space = item_to_window (Duple (0.0, y)); /* x-value is irrelevant */ + + /* we are moving across the draw area pixel-by-pixel */ + window_space.x = draw.x0 + pixel; + + /* plot this point */ + if (pixel == 0) { + context->move_to (window_space.x, window_space.y); + } else { + context->line_to (window_space.x, window_space.y); + } + } + + context->stroke (); + context->restore (); } - delete [] tmp; +#if 0 + /* add points */ - return x; + setup_fill_context (context); + for (Points::const_iterator p = _points.begin(); p != _points.end(); ++p) { + Duple window_space (item_to_window (*p)); + context->arc (window_space.x, window_space.y, 5.0, 0.0, 2 * M_PI); + context->stroke (); + } +#endif } bool @@ -215,7 +419,7 @@ Curve::covers (Duple const & pc) const { Duple point = canvas_to_item (pc); - /* XXX Hellaciously expensive ... */ + /* O(N) N = number of points, and not accurate */ for (Points::const_iterator p = _points.begin(); p != _points.end(); ++p) {