*/
+#include <cmath>
#include <exception>
#include <algorithm>
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
{
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<Cairo::Context> 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<Cairo::Context> 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<double> 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<Cairo::Context> 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<Rect> 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<double> const & rhs)
-{
- std::vector<double>::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<double>::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<double>::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
{
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) {
projects / ardour.git / blobdiff