libs/canvas/curve.cc

   1 /*
   2     Copyright (C) 2013 Paul Davis
   3
   4     This program is free software; you can redistribute it and/or modify
   5     it under the terms of the GNU General Public License as published by
   6     the Free Software Foundation; either version 2 of the License, or
   7     (at your option) any later version.
   8
   9     This program is distributed in the hope that it will be useful,
  10     but WITHOUT ANY WARRANTY; without even the implied warranty of
  11     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  12     GNU General Public License for more details.
  13
  14     You should have received a copy of the GNU General Public License
  15     along with this program; if not, write to the Free Software
  16     Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
  17
  18 */
  19
  20 #include <cmath>
  21 #include <exception>
  22 #include <algorithm>
  23
  24 #include "canvas/curve.h"
  25
  26 using namespace ArdourCanvas;
  27 using std::min;
  28 using std::max;
  29
  30 Curve::Curve (Group* parent)
  31         : Item (parent)
  32         , PolyItem (parent)
  33         , Fill (parent)
  34         , n_samples (0)
  35         , points_per_segment (16)
  36         , curve_type (CatmullRomCentripetal)
  37         , curve_fill (None)
  38 {
  39 }
  40
  41 /** When rendering the curve, we will always draw a fixed number of straight
  42  * line segments to span the x-axis extent of the curve. More segments:
  43  * smoother visual rendering. Less rendering: closer to a visibily poly-line
  44  * render.
  45  */
  46 void
  47 Curve::set_points_per_segment (uint32_t n)
  48 {
  49         /* this only changes our appearance rather than the bounding box, so we
  50            just need to schedule a redraw rather than notify the parent of any
  51            changes
  52         */
  53         points_per_segment = n;
  54         interpolate ();
  55         redraw ();
  56 }
  57
  58 void
  59 Curve::compute_bounding_box () const
  60 {
  61         PolyItem::compute_bounding_box ();
  62
  63         /* possibly add extents of any point indicators here if we ever do that */
  64 }
  65
  66 void
  67 Curve::set (Points const& p)
  68 {
  69         PolyItem::set (p);
  70         interpolate ();
  71 }
  72
  73 void
  74 Curve::interpolate ()
  75 {
  76         samples.clear ();
  77         interpolate (_points, points_per_segment, CatmullRomCentripetal, false, samples);
  78         n_samples = samples.size();
  79 }
  80
  81 /* Cartmull-Rom code from http://stackoverflow.com/questions/9489736/catmull-rom-curve-with-no-cusps-and-no-self-intersections/19283471#19283471
  82  *
  83  * Thanks to Ted for his Java version, which I translated into Ardour-idiomatic
  84  * C++ here.
  85  */
  86
  87 /**
  88  * Calculate the same values but introduces the ability to "parameterize" the t
  89  * values used in the calculation. This is based on Figure 3 from
  90  * http://www.cemyuksel.com/research/catmullrom_param/catmullrom.pdf
  91  *
  92  * @param p An array of double values of length 4, where interpolation
  93  * occurs from p1 to p2.
  94  * @param time An array of time measures of length 4, corresponding to each
  95  * p value.
  96  * @param t the actual interpolation ratio from 0 to 1 representing the
  97  * position between p1 and p2 to interpolate the value.
  98  */
  99 static double
 100 __interpolate (double p[4], double time[4], double t)
 101 {
 102         const double L01 = p[0] * (time[1] - t) / (time[1] - time[0]) + p[1] * (t - time[0]) / (time[1] - time[0]);
 103         const double L12 = p[1] * (time[2] - t) / (time[2] - time[1]) + p[2] * (t - time[1]) / (time[2] - time[1]);
 104         const double L23 = p[2] * (time[3] - t) / (time[3] - time[2]) + p[3] * (t - time[2]) / (time[3] - time[2]);
 105         const double L012 = L01 * (time[2] - t) / (time[2] - time[0]) + L12 * (t - time[0]) / (time[2] - time[0]);
 106         const double L123 = L12 * (time[3] - t) / (time[3] - time[1]) + L23 * (t - time[1]) / (time[3] - time[1]);
 107         const double C12 = L012 * (time[2] - t) / (time[2] - time[1]) + L123 * (t - time[1]) / (time[2] - time[1]);
 108         return C12;
 109 }
 110
 111 /**
 112  * Given a list of control points, this will create a list of points_per_segment
 113  * points spaced uniformly along the resulting Catmull-Rom curve.
 114  *
 115  * @param points The list of control points, leading and ending with a
 116  * coordinate that is only used for controling the spline and is not visualized.
 117  * @param index The index of control point p0, where p0, p1, p2, and p3 are
 118  * used in order to create a curve between p1 and p2.
 119  * @param points_per_segment The total number of uniformly spaced interpolated
 120  * points to calculate for each segment. The larger this number, the
 121  * smoother the resulting curve.
 122  * @param curve_type Clarifies whether the curve should use uniform, chordal
 123  * or centripetal curve types. Uniform can produce loops, chordal can
 124  * produce large distortions from the original lines, and centripetal is an
 125  * optimal balance without spaces.
 126  * @return the list of coordinates that define the CatmullRom curve
 127  * between the points defined by index+1 and index+2.
 128  */
 129 static void
 130 _interpolate (const Points& points, Points::size_type index, int points_per_segment, Curve::SplineType curve_type, Points& results)
 131 {
 132         double x[4];
 133         double y[4];
 134         double time[4];
 135
 136         for (int i = 0; i < 4; i++) {
 137                 x[i] = points[index + i].x;
 138                 y[i] = points[index + i].y;
 139                 time[i] = i;
 140         }
 141
 142         double tstart = 1;
 143         double tend = 2;
 144
 145         if (curve_type != Curve::CatmullRomUniform) {
 146                 double total = 0;
 147                 for (int i = 1; i < 4; i++) {
 148                         double dx = x[i] - x[i - 1];
 149                         double dy = y[i] - y[i - 1];
 150                         if (curve_type == Curve::CatmullRomCentripetal) {
 151                                 total += pow (dx * dx + dy * dy, .25);
 152                         } else {
 153                                 total += pow (dx * dx + dy * dy, .5);
 154                         }
 155                         time[i] = total;
 156                 }
 157                 tstart = time[1];
 158                 tend = time[2];
 159         }
 160
 161         int segments = points_per_segment - 1;
 162         results.push_back (points[index + 1]);
 163
 164         for (int i = 1; i < segments; i++) {
 165                 double xi = __interpolate (x, time, tstart + (i * (tend - tstart)) / segments);
 166                 double yi = __interpolate (y, time, tstart + (i * (tend - tstart)) / segments);
 167                 results.push_back (Duple (xi, yi));
 168         }
 169
 170         results.push_back (points[index + 2]);
 171 }
 172
 173 /**
 174  * This method will calculate the Catmull-Rom interpolation curve, returning
 175  * it as a list of Coord coordinate objects.  This method in particular
 176  * adds the first and last control points which are not visible, but required
 177  * for calculating the spline.
 178  *
 179  * @param coordinates The list of original straight line points to calculate
 180  * an interpolation from.
 181  * @param points_per_segment The integer number of equally spaced points to
 182  * return along each curve.  The actual distance between each
 183  * point will depend on the spacing between the control points.
 184  * @return The list of interpolated coordinates.
 185  * @param curve_type Chordal (stiff), Uniform(floppy), or Centripetal(medium)
 186  * @throws gov.ca.water.shapelite.analysis.CatmullRomException if
 187  * points_per_segment is less than 2.
 188  */
 189
 190 void
 191 Curve::interpolate (const Points& coordinates, uint32_t points_per_segment, SplineType curve_type, bool closed, Points& results)
 192 {
 193         if (points_per_segment < 2) {
 194                 return;
 195         }
 196
 197         // Cannot interpolate curves given only two points.  Two points
 198         // is best represented as a simple line segment.
 199         if (coordinates.size() < 3) {
 200                 results = coordinates;
 201                 return;
 202         }
 203
 204         // Copy the incoming coordinates. We need to modify it during interpolation
 205         Points vertices = coordinates;
 206
 207         // Test whether the shape is open or closed by checking to see if
 208         // the first point intersects with the last point.  M and Z are ignored.
 209         if (closed) {
 210                 // Use the second and second from last points as control points.
 211                 // get the second point.
 212                 Duple p2 = vertices[1];
 213                 // get the point before the last point
 214                 Duple pn1 = vertices[vertices.size() - 2];
 215
 216                 // insert the second from the last point as the first point in the list
 217                 // because when the shape is closed it keeps wrapping around to
 218                 // the second point.
 219                 vertices.insert(vertices.begin(), pn1);
 220                 // add the second point to the end.
 221                 vertices.push_back(p2);
 222         } else {
 223                 // The shape is open, so use control points that simply extend
 224                 // the first and last segments
 225
 226                 // Get the change in x and y between the first and second coordinates.
 227                 double dx = vertices[1].x - vertices[0].x;
 228                 double dy = vertices[1].y - vertices[0].y;
 229
 230                 // Then using the change, extrapolate backwards to find a control point.
 231                 double x1 = vertices[0].x - dx;
 232                 double y1 = vertices[0].y - dy;
 233
 234                 // Actaully create the start point from the extrapolated values.
 235                 Duple start (x1, y1);
 236
 237                 // Repeat for the end control point.
 238                 int n = vertices.size() - 1;
 239                 dx = vertices[n].x - vertices[n - 1].x;
 240                 dy = vertices[n].y - vertices[n - 1].y;
 241                 double xn = vertices[n].x + dx;
 242                 double yn = vertices[n].y + dy;
 243                 Duple end (xn, yn);
 244
 245                 // insert the start control point at the start of the vertices list.
 246                 vertices.insert (vertices.begin(), start);
 247
 248                 // append the end control ponit to the end of the vertices list.
 249                 vertices.push_back (end);
 250         }
 251
 252         // When looping, remember that each cycle requires 4 points, starting
 253         // with i and ending with i+3.  So we don't loop through all the points.
 254
 255         for (Points::size_type i = 0; i < vertices.size() - 3; i++) {
 256
 257                 // Actually calculate the Catmull-Rom curve for one segment.
 258                 Points r;
 259
 260                 _interpolate (vertices, i, points_per_segment, curve_type, r);
 261
 262                 // Since the middle points are added twice, once for each bordering
 263                 // segment, we only add the 0 index result point for the first
 264                 // segment.  Otherwise we will have duplicate points.
 265
 266                 if (results.size() > 0) {
 267                         r.erase (r.begin());
 268                 }
 269
 270                 // Add the coordinates for the segment to the result list.
 271
 272                 results.insert (results.end(), r.begin(), r.end());
 273         }
 274 }
 275
 276 void
 277 Curve::render (Rect const & area, Cairo::RefPtr<Cairo::Context> context) const
 278 {
 279         if (!_outline || _points.size() < 2 || !_bounding_box) {
 280                 return;
 281         }
 282
 283         Rect self = item_to_window (_bounding_box.get());
 284         boost::optional<Rect> d = self.intersection (area);
 285         assert (d);
 286         Rect draw = d.get ();
 287
 288         /* Our approach is to always draw n_segments across our total size.
 289          *
 290          * This is very inefficient if we are asked to only draw a small
 291          * section of the curve. For now we rely on cairo clipping to help
 292          * with this.
 293          */
 294
 295
 296         setup_outline_context (context);
 297
 298         if (_points.size() == 2) {
 299
 300                 /* straight line */
 301
 302                 Duple window_space;
 303
 304                 window_space = item_to_window (_points.front());
 305                 context->move_to (window_space.x, window_space.y);
 306                 window_space = item_to_window (_points.back());
 307                 context->line_to (window_space.x, window_space.y);
 308
 309
 310                 switch (curve_fill) {
 311                         case None:
 312                                 context->stroke();
 313                                 break;
 314                         case Inside:
 315                                 context->stroke_preserve ();
 316                                 window_space = item_to_window (Duple(_points.back().x, draw.height()));
 317                                 context->line_to (window_space.x, window_space.y);
 318                                 window_space = item_to_window (Duple(_points.front().x, draw.height()));
 319                                 context->line_to (window_space.x, window_space.y);
 320                                 context->close_path();
 321                                 setup_fill_context(context);
 322                                 context->fill ();
 323                                 break;
 324                         case Outside:
 325                                 context->stroke_preserve ();
 326                                 window_space = item_to_window (Duple(_points.back().x, 0.0));
 327                                 context->line_to (window_space.x, window_space.y);
 328                                 window_space = item_to_window (Duple(_points.front().x, 0.0));
 329                                 context->line_to (window_space.x, window_space.y);
 330                                 context->close_path();
 331                                 setup_fill_context(context);
 332                                 context->fill ();
 333                                 break;
 334                 }
 335
 336         } else {
 337
 338                 /* curve of at least 3 points */
 339
 340                 /* x-axis limits of the curve, in window space coordinates */
 341
 342                 Duple w1 = item_to_window (Duple (_points.front().x, 0.0));
 343                 Duple w2 = item_to_window (Duple (_points.back().x, 0.0));
 344
 345                 /* clamp actual draw to area bound by points, rather than our bounding box which is slightly different */
 346
 347                 context->save ();
 348                 context->rectangle (draw.x0, draw.y0, draw.width(), draw.height());
 349                 context->clip ();
 350
 351                 /* expand drawing area by several pixels on each side to avoid cairo stroking effects at the boundary.
 352                    they will still occur, but cairo's clipping will hide them.
 353                  */
 354
 355                 draw = draw.expand (4.0);
 356
 357                 /* now clip it to the actual points in the curve */
 358
 359                 if (draw.x0 < w1.x) {
 360                         draw.x0 = w1.x;
 361                 }
 362
 363                 if (draw.x1 >= w2.x) {
 364                         draw.x1 = w2.x;
 365                 }
 366
 367                 /* find left and right-most sample */
 368                 Duple window_space;
 369                 Points::size_type left = 0;
 370                 Points::size_type right = n_samples;
 371
 372                 for (Points::size_type idx = 0; idx < n_samples - 1; ++idx) {
 373                         left = idx;
 374                         window_space = item_to_window (Duple (samples[idx].x, 0.0));
 375                         if (window_space.x >= draw.x0) break;
 376                 }
 377                 for (Points::size_type idx = n_samples; idx > left + 1; --idx) {
 378                         window_space = item_to_window (Duple (samples[idx].x, 0.0));
 379                         if (window_space.x <= draw.x1) break;
 380                         right = idx;
 381                 }
 382
 383                 /* draw line between samples */
 384                 window_space = item_to_window (Duple (samples[left].x, samples[left].y));
 385                 context->move_to (window_space.x, window_space.y);
 386                 Coord last_x = round(window_space.x);
 387                 for (uint32_t idx = left + 1; idx < right; ++idx) {
 388                         window_space = item_to_window (Duple (samples[idx].x, samples[idx].y));
 389                         if (last_x == round(window_space.x)) continue;
 390                         last_x = round(window_space.x);
 391                         context->line_to (last_x - .5 , window_space.y);
 392                 }
 393
 394                 switch (curve_fill) {
 395                         case None:
 396                                 context->stroke();
 397                                 break;
 398                         case Inside:
 399                                 context->stroke_preserve ();
 400                                 window_space = item_to_window (Duple (samples[right-1].x, draw.height()));
 401                                 context->line_to (window_space.x, window_space.y);
 402                                 window_space = item_to_window (Duple (samples[left].x, draw.height()));
 403                                 context->line_to (window_space.x, window_space.y);
 404                                 context->close_path();
 405                                 setup_fill_context(context);
 406                                 context->fill ();
 407                                 break;
 408                         case Outside:
 409                                 context->stroke_preserve ();
 410                                 window_space = item_to_window (Duple (samples[right-1].x, 0.0));
 411                                 context->line_to (window_space.x, window_space.y);
 412                                 window_space = item_to_window (Duple (samples[left].x, 0.0));
 413                                 context->line_to (window_space.x, window_space.y);
 414                                 context->close_path();
 415                                 setup_fill_context(context);
 416                                 context->fill ();
 417                                 break;
 418                 }
 419                 context->restore ();
 420         }
 421
 422 #if 0
 423         /* add points */
 424         setup_outline_context (context);
 425         for (Points::const_iterator p = _points.begin(); p != _points.end(); ++p) {
 426                 Duple window_space (item_to_window (*p));
 427                 context->arc (window_space.x, window_space.y, 5.0, 0.0, 2 * M_PI);
 428                 context->stroke ();
 429         }
 430 #endif
 431 }
 432
 433 bool
 434 Curve::covers (Duple const & pc) const
 435 {
 436         Duple point = canvas_to_item (pc);
 437
 438         /* O(N) N = number of points, and not accurate */
 439
 440         for (Points::const_iterator p = _points.begin(); p != _points.end(); ++p) {
 441
 442                 const Coord dx = point.x - (*p).x;
 443                 const Coord dy = point.y - (*p).y;
 444                 const Coord dx2 = dx * dx;
 445                 const Coord dy2 = dy * dy;
 446
 447                 if ((dx2 < 2.0 && dy2 < 2.0) || (dx2 + dy2 < 4.0)) {
 448                         return true;
 449                 }
 450         }
 451
 452         return false;
 453 }