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 <exception>
  21 #include <algorithm>
  22
  23 #include "canvas/curve.h"
  24
  25 using namespace ArdourCanvas;
  26 using std::min;
  27 using std::max;
  28
  29 Curve::Curve (Group* parent)
  30         : Item (parent)
  31         , PolyItem (parent)
  32         , Fill (parent)
  33         , n_samples (0)
  34         , n_segments (512)
  35 {
  36         set_n_samples (256);
  37 }
  38
  39 /** Set the number of points to compute when we smooth the data points into a
  40  * curve.
  41  */
  42 void
  43 Curve::set_n_samples (Points::size_type n)
  44 {
  45         /* this only changes our appearance rather than the bounding box, so we
  46            just need to schedule a redraw rather than notify the parent of any
  47            changes
  48         */
  49         n_samples = n;
  50         samples.assign (n_samples, Duple (0.0, 0.0));
  51         interpolate ();
  52 }
  53
  54 /** When rendering the curve, we will always draw a fixed number of straight
  55  * line segments to span the x-axis extent of the curve. More segments:
  56  * smoother visual rendering. Less rendering: closer to a visibily poly-line
  57  * render.
  58  */
  59 void
  60 Curve::set_n_segments (uint32_t n)
  61 {
  62         /* this only changes our appearance rather than the bounding box, so we
  63            just need to schedule a redraw rather than notify the parent of any
  64            changes
  65         */
  66         n_segments = n;
  67         redraw ();
  68 }
  69
  70 void
  71 Curve::compute_bounding_box () const
  72 {
  73         PolyItem::compute_bounding_box ();
  74
  75         /* possibly add extents of any point indicators here if we ever do that */
  76 }
  77
  78 void
  79 Curve::set (Points const& p)
  80 {
  81         PolyItem::set (p);
  82         interpolate ();
  83 }
  84
  85 void
  86 Curve::interpolate ()
  87 {
  88         Points::size_type npoints = _points.size ();
  89
  90         if (npoints < 3) {
  91                 return;
  92         }
  93
  94         Duple p;
  95         double boundary;
  96
  97         const double xfront = _points.front().x;
  98         const double xextent = _points.back().x - xfront;
  99
 100         /* initialize boundary curve points */
 101
 102         p = _points.front();
 103         boundary = round (((p.x - xfront)/xextent) * (n_samples - 1));
 104
 105         for (Points::size_type i = 0; i < boundary; ++i) {
 106                 samples[i] = Duple (i, p.y);
 107         }
 108
 109         p = _points.back();
 110         boundary = round (((p.x - xfront)/xextent) * (n_samples - 1));
 111
 112         for (Points::size_type i = boundary; i < n_samples; ++i) {
 113                 samples[i] = Duple (i, p.y);
 114         }
 115
 116         for (int i = 0; i < (int) npoints - 1; ++i) {
 117
 118                 Points::size_type p1, p2, p3, p4;
 119
 120                 p1 = max (i - 1, 0);
 121                 p2 = i;
 122                 p3 = i + 1;
 123                 p4 = min (i + 2, (int) npoints - 1);
 124
 125                 smooth (p1, p2, p3, p4, xfront, xextent);
 126         }
 127
 128         /* make sure that actual data points are used with their exact values */
 129
 130         for (Points::const_iterator p = _points.begin(); p != _points.end(); ++p) {
 131                 uint32_t idx = (((*p).x - xfront)/xextent) * (n_samples - 1);
 132                 samples[idx].y = (*p).y;
 133         }
 134 }
 135
 136 /*
 137  * This function calculates the curve values between the control points
 138  * p2 and p3, taking the potentially existing neighbors p1 and p4 into
 139  * account.
 140  *
 141  * This function uses a cubic bezier curve for the individual segments and
 142  * calculates the necessary intermediate control points depending on the
 143  * neighbor curve control points.
 144  *
 145  */
 146
 147 void
 148 Curve::smooth (Points::size_type p1, Points::size_type p2, Points::size_type p3, Points::size_type p4,
 149                double xfront, double xextent)
 150 {
 151         int i;
 152         double x0, x3;
 153         double y0, y1, y2, y3;
 154         double dx, dy;
 155         double slope;
 156
 157         /* the outer control points for the bezier curve. */
 158
 159         x0 = _points[p2].x;
 160         y0 = _points[p2].y;
 161         x3 = _points[p3].x;
 162         y3 = _points[p3].y;
 163
 164         /*
 165          * the x values of the inner control points are fixed at
 166          * x1 = 2/3*x0 + 1/3*x3   and  x2 = 1/3*x0 + 2/3*x3
 167          * this ensures that the x values increase linearily with the
 168          * parameter t and enables us to skip the calculation of the x
 169          * values altogehter - just calculate y(t) evenly spaced.
 170          */
 171
 172         dx = x3 - x0;
 173         dy = y3 - y0;
 174
 175         if (dx <= 0) {
 176                 /* error? */
 177                 return;
 178         }
 179
 180         if (p1 == p2 && p3 == p4) {
 181                 /* No information about the neighbors,
 182                  * calculate y1 and y2 to get a straight line
 183                  */
 184                 y1 = y0 + dy / 3.0;
 185                 y2 = y0 + dy * 2.0 / 3.0;
 186
 187         } else if (p1 == p2 && p3 != p4) {
 188
 189                 /* only the right neighbor is available. Make the tangent at the
 190                  * right endpoint parallel to the line between the left endpoint
 191                  * and the right neighbor. Then point the tangent at the left towards
 192                  * the control handle of the right tangent, to ensure that the curve
 193                  * does not have an inflection point.
 194                  */
 195                 slope = (_points[p4].y - y0) / (_points[p4].x - x0);
 196
 197                 y2 = y3 - slope * dx / 3.0;
 198                 y1 = y0 + (y2 - y0) / 2.0;
 199
 200         } else if (p1 != p2 && p3 == p4) {
 201
 202                 /* see previous case */
 203                 slope = (y3 - _points[p1].y) / (x3 - _points[p1].x);
 204
 205                 y1 = y0 + slope * dx / 3.0;
 206                 y2 = y3 + (y1 - y3) / 2.0;
 207
 208
 209         } else /* (p1 != p2 && p3 != p4) */ {
 210
 211                 /* Both neighbors are available. Make the tangents at the endpoints
 212                  * parallel to the line between the opposite endpoint and the adjacent
 213                  * neighbor.
 214                  */
 215
 216                 slope = (y3 - _points[p1].y) / (x3 - _points[p1].x);
 217
 218                 y1 = y0 + slope * dx / 3.0;
 219
 220                 slope = (_points[p4].y - y0) / (_points[p4].x - x0);
 221
 222                 y2 = y3 - slope * dx / 3.0;
 223         }
 224
 225         /*
 226          * finally calculate the y(t) values for the given bezier values. We can
 227          * use homogenously distributed values for t, since x(t) increases linearily.
 228          */
 229
 230         dx = dx / xextent;
 231
 232         int limit = round (dx * (n_samples - 1));
 233         const int idx_offset = round (((x0 - xfront)/xextent) * (n_samples - 1));
 234
 235         for (i = 0; i <= limit; i++) {
 236                 double y, t;
 237                 Points::size_type index;
 238
 239                 t = i / dx / (n_samples - 1);
 240
 241                 y =     y0 * (1-t) * (1-t) * (1-t) +
 242                         3 * y1 * (1-t) * (1-t) * t     +
 243                         3 * y2 * (1-t) * t     * t     +
 244                         y3 * t     * t     * t;
 245
 246                 index = i + idx_offset;
 247
 248                 if (index < n_samples) {
 249                         Duple d (i, max (y, 0.0));
 250                         samples[index] = d;
 251                 }
 252         }
 253 }
 254
 255 /** Given a fractional position within the x-axis range of the
 256  * curve, return the corresponding y-axis value
 257  */
 258
 259 double
 260 Curve::map_value (double x) const
 261 {
 262         if (x > 0.0 && x < 1.0) {
 263
 264                 double f;
 265                 Points::size_type index;
 266
 267                 /* linearly interpolate between two of our smoothed "samples"
 268                  */
 269
 270                 x = x * (n_samples - 1);
 271                 index = (Points::size_type) x; // XXX: should we explicitly use floor()?
 272                 f = x - index;
 273
 274                 return (1.0 - f) * samples[index].y + f * samples[index+1].y;
 275
 276         } else if (x >= 1.0) {
 277                 return samples.back().y;
 278         } else {
 279                 return samples.front().y;
 280         }
 281 }
 282
 283 void
 284 Curve::render (Rect const & area, Cairo::RefPtr<Cairo::Context> context) const
 285 {
 286         if (!_outline || _points.size() < 2 || !_bounding_box) {
 287                 return;
 288         }
 289
 290         Rect self = item_to_window (_bounding_box.get());
 291         boost::optional<Rect> d = self.intersection (area);
 292         assert (d);
 293         Rect draw = d.get ();
 294
 295         /* Our approach is to always draw n_segments across our total size.
 296          *
 297          * This is very inefficient if we are asked to only draw a small
 298          * section of the curve. For now we rely on cairo clipping to help
 299          * with this.
 300          */
 301
 302
 303         setup_outline_context (context);
 304
 305         if (_points.size() == 2) {
 306
 307                 /* straight line */
 308
 309                 Duple window_space;
 310
 311                 window_space = item_to_window (_points.front());
 312                 context->move_to (window_space.x, window_space.y);
 313                 window_space = item_to_window (_points.back());
 314                 context->line_to (window_space.x, window_space.y);
 315
 316                 context->stroke ();
 317
 318         } else {
 319
 320                 /* curve of at least 3 points */
 321
 322                 /* x-axis limits of the curve, in window space coordinates */
 323
 324                 Duple w1 = item_to_window (Duple (_points.front().x, 0.0));
 325                 Duple w2 = item_to_window (Duple (_points.back().x, 0.0));
 326
 327                 /* clamp actual draw to area bound by points, rather than our bounding box which is slightly different */
 328
 329                 context->save ();
 330                 context->rectangle (draw.x0, draw.y0, draw.width(), draw.height());
 331                 context->clip ();
 332
 333                 /* expand drawing area by several pixels on each side to avoid cairo stroking effects at the boundary.
 334                    they will still occur, but cairo's clipping will hide them.
 335                  */
 336
 337                 draw = draw.expand (4.0);
 338
 339                 /* now clip it to the actual points in the curve */
 340
 341                 if (draw.x0 < w1.x) {
 342                         draw.x0 = w1.x;
 343                 }
 344
 345                 if (draw.x1 >= w2.x) {
 346                         draw.x1 = w2.x;
 347                 }
 348
 349                 /* full width of the curve */
 350                 const double xextent = _points.back().x - _points.front().x;
 351                 /* Determine where the first drawn point will be */
 352                 Duple item_space = window_to_item (Duple (draw.x0, 0)); /* y value is irrelevant */
 353                 /* determine the fractional offset of this location into the overall extent of the curve */
 354                 const double xfract_offset = (item_space.x - _points.front().x)/xextent;
 355                 const uint32_t pixels = draw.width ();
 356                 Duple window_space;
 357
 358                 /* draw the first point */
 359
 360                 for (uint32_t pixel = 0; pixel < pixels; ++pixel) {
 361
 362                         /* fractional distance into the total horizontal extent of the curve */
 363                         double xfract = xfract_offset + (pixel / xextent);
 364                         /* compute vertical coordinate (item-space) at that location */
 365                         double y = map_value (xfract);
 366
 367                         /* convert to window space for drawing */
 368                         window_space = item_to_window (Duple (0.0, y)); /* x-value is irrelevant */
 369
 370                         /* we are moving across the draw area pixel-by-pixel */
 371                         window_space.x = draw.x0 + pixel;
 372
 373                         /* plot this point */
 374                         if (pixel == 0) {
 375                                 context->move_to (window_space.x, window_space.y);
 376                         } else {
 377                                 context->line_to (window_space.x, window_space.y);
 378                         }
 379                 }
 380
 381                 context->stroke ();
 382                 context->restore ();
 383         }
 384
 385 #if 0
 386         /* add points */
 387
 388         setup_fill_context (context);
 389         for (Points::const_iterator p = _points.begin(); p != _points.end(); ++p) {
 390                 Duple window_space (item_to_window (*p));
 391                 context->arc (window_space.x, window_space.y, 5.0, 0.0, 2 * M_PI);
 392                 context->stroke ();
 393         }
 394 #endif
 395 }
 396
 397 bool
 398 Curve::covers (Duple const & pc) const
 399 {
 400         Duple point = canvas_to_item (pc);
 401
 402         /* O(N) N = number of points, and not accurate */
 403
 404         for (Points::const_iterator p = _points.begin(); p != _points.end(); ++p) {
 405
 406                 const Coord dx = point.x - (*p).x;
 407                 const Coord dy = point.y - (*p).y;
 408                 const Coord dx2 = dx * dx;
 409                 const Coord dy2 = dy * dy;
 410
 411                 if ((dx2 < 2.0 && dy2 < 2.0) || (dx2 + dy2 < 4.0)) {
 412                         return true;
 413                 }
 414         }
 415
 416         return false;
 417 }