Paul Davis
58a30da03d
See http://en.wikipedia.org/wiki/Centripetal_Catmull-Rom to understand the benefits of this.
438 lines
14 KiB
C++
438 lines
14 KiB
C++
/*
|
|
Copyright (C) 2013 Paul Davis
|
|
|
|
This program is free software; you can redistribute it and/or modify
|
|
it under the terms of the GNU General Public License as published by
|
|
the Free Software Foundation; either version 2 of the License, or
|
|
(at your option) any later version.
|
|
|
|
This program is distributed in the hope that it will be useful,
|
|
but WITHOUT ANY WARRANTY; without even the implied warranty of
|
|
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
|
GNU General Public License for more details.
|
|
|
|
You should have received a copy of the GNU General Public License
|
|
along with this program; if not, write to the Free Software
|
|
Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
|
|
|
|
*/
|
|
|
|
#include <cmath>
|
|
#include <exception>
|
|
#include <algorithm>
|
|
|
|
#include "canvas/curve.h"
|
|
|
|
using namespace ArdourCanvas;
|
|
using std::min;
|
|
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
|
|
Curve::compute_bounding_box () const
|
|
{
|
|
PolyItem::compute_bounding_box ();
|
|
|
|
/* possibly add extents of any point indicators here if we ever do that */
|
|
}
|
|
|
|
void
|
|
Curve::set (Points const& p)
|
|
{
|
|
PolyItem::set (p);
|
|
interpolate ();
|
|
}
|
|
|
|
void
|
|
Curve::interpolate ()
|
|
{
|
|
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]);
|
|
}
|
|
|
|
/**
|
|
* 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)
|
|
{
|
|
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
|
|
* curve, return the corresponding y-axis value
|
|
*/
|
|
|
|
double
|
|
Curve::map_value (double x) const
|
|
{
|
|
if (x > 0.0 && x < 1.0) {
|
|
|
|
double f;
|
|
Points::size_type index;
|
|
|
|
/* linearly interpolate between two of our smoothed "samples"
|
|
*/
|
|
|
|
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;
|
|
|
|
} else if (x >= 1.0) {
|
|
return samples.back().y;
|
|
} else {
|
|
return samples.front().y;
|
|
}
|
|
}
|
|
|
|
void
|
|
Curve::render (Rect const & area, Cairo::RefPtr<Cairo::Context> context) const
|
|
{
|
|
if (!_outline || _points.size() < 2 || !_bounding_box) {
|
|
return;
|
|
}
|
|
|
|
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.
|
|
*/
|
|
|
|
|
|
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 (draw.x0 < w1.x) {
|
|
draw.x0 = w1.x;
|
|
}
|
|
|
|
if (draw.x1 >= w2.x) {
|
|
draw.x1 = w2.x;
|
|
}
|
|
|
|
/* 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;
|
|
|
|
/* draw the first point */
|
|
|
|
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 ();
|
|
}
|
|
|
|
#if 0
|
|
/* add points */
|
|
|
|
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
|
|
Curve::covers (Duple const & pc) const
|
|
{
|
|
Duple point = canvas_to_item (pc);
|
|
|
|
/* O(N) N = number of points, and not accurate */
|
|
|
|
for (Points::const_iterator p = _points.begin(); p != _points.end(); ++p) {
|
|
|
|
const Coord dx = point.x - (*p).x;
|
|
const Coord dy = point.y - (*p).y;
|
|
const Coord dx2 = dx * dx;
|
|
const Coord dy2 = dy * dy;
|
|
|
|
if ((dx2 < 2.0 && dy2 < 2.0) || (dx2 + dy2 < 4.0)) {
|
|
return true;
|
|
}
|
|
}
|
|
|
|
return false;
|
|
}
|