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 "pbd/xml++.h"
  24
  25 #include "canvas/curve.h"
  26
  27 using namespace ArdourCanvas;
  28 using std::min;
  29 using std::max;
  30
  31 Curve::Curve (Group* parent)
  32         : Item (parent)
  33         , PolyItem (parent)
  34 {
  35
  36 }
  37
  38 void
  39 Curve::compute_bounding_box () const
  40 {
  41         PolyItem::compute_bounding_box ();
  42
  43         if (_bounding_box) {
  44
  45                 bool have1 = false;
  46                 bool have2 = false;
  47
  48                 Rect bbox1;
  49                 Rect bbox2;
  50
  51                 for (Points::const_iterator i = first_control_points.begin(); i != first_control_points.end(); ++i) {
  52                         if (have1) {
  53                                 bbox1.x0 = min (bbox1.x0, i->x);
  54                                 bbox1.y0 = min (bbox1.y0, i->y);
  55                                 bbox1.x1 = max (bbox1.x1, i->x);
  56                                 bbox1.y1 = max (bbox1.y1, i->y);
  57                         } else {
  58                                 bbox1.x0 = bbox1.x1 = i->x;
  59                                 bbox1.y0 = bbox1.y1 = i->y;
  60                                 have1 = true;
  61                         }
  62                 }
  63
  64                 for (Points::const_iterator i = second_control_points.begin(); i != second_control_points.end(); ++i) {
  65                         if (have2) {
  66                                 bbox2.x0 = min (bbox2.x0, i->x);
  67                                 bbox2.y0 = min (bbox2.y0, i->y);
  68                                 bbox2.x1 = max (bbox2.x1, i->x);
  69                                 bbox2.y1 = max (bbox2.y1, i->y);
  70                         } else {
  71                                 bbox2.x0 = bbox2.x1 = i->x;
  72                                 bbox2.y0 = bbox2.y1 = i->y;
  73                                 have2 = true;
  74                         }
  75                 }
  76
  77                 Rect u = bbox1.extend (bbox2);
  78                 _bounding_box = u.extend (_bounding_box.get());
  79         }
  80
  81         _bounding_box_dirty = false;
  82 }
  83
  84 void
  85 Curve::set (Points const& p)
  86 {
  87         PolyItem::set (p);
  88
  89         first_control_points.clear ();
  90         second_control_points.clear ();
  91
  92         compute_control_points (_points, first_control_points, second_control_points);
  93 }
  94
  95 void
  96 Curve::render (Rect const & area, Cairo::RefPtr<Cairo::Context> context) const
  97 {
  98         if (_outline) {
  99                 setup_outline_context (context);
 100                 render_path (area, context);
 101                 context->stroke ();
 102         }
 103 }
 104
 105 XMLNode *
 106 Curve::get_state () const
 107 {
 108         XMLNode* node = new XMLNode ("PolyLine");
 109         add_poly_item_state (node);
 110         add_outline_state (node);
 111         return node;
 112 }
 113
 114 void
 115 Curve::set_state (XMLNode const * node)
 116 {
 117         set_poly_item_state (node);
 118         set_outline_state (node);
 119 }
 120
 121 void
 122 Curve::render_path (Rect const & area, Cairo::RefPtr<Cairo::Context> context) const
 123 {
 124         PolyItem::render_curve (area, context, first_control_points, second_control_points);
 125 }
 126
 127 void
 128 Curve::compute_control_points (Points const& knots,
 129                                Points& firstControlPoints,
 130                                Points& secondControlPoints)
 131 {
 132         Points::size_type n = knots.size() - 1;
 133
 134         if (n < 1) {
 135                 return;
 136         }
 137
 138         if (n == 1) {
 139                 /* Special case: Bezier curve should be a straight line. */
 140
 141                 Duple d;
 142
 143                 d.x = (2.0 * knots[0].x + knots[1].x) / 3;
 144                 d.y = (2.0 * knots[0].y + knots[1].y) / 3;
 145                 firstControlPoints.push_back (d);
 146
 147                 d.x = 2.0 * firstControlPoints[0].x - knots[0].x;
 148                 d.y = 2.0 * firstControlPoints[0].y - knots[0].y;
 149                 secondControlPoints.push_back (d);
 150
 151                 return;
 152         }
 153
 154         // Calculate first Bezier control points
 155         // Right hand side vector
 156
 157         std::vector<double> rhs;
 158
 159         rhs.assign (n, 0);
 160
 161         // Set right hand side X values
 162
 163         for (Points::size_type i = 1; i < n - 1; ++i) {
 164                 rhs[i] = 4 * knots[i].x + 2 * knots[i + 1].x;
 165         }
 166         rhs[0] = knots[0].x + 2 * knots[1].x;
 167         rhs[n - 1] = (8 * knots[n - 1].x + knots[n].x) / 2.0;
 168
 169         // Get first control points X-values
 170         double* x = solve (rhs);
 171
 172         // Set right hand side Y values
 173         for (Points::size_type i = 1; i < n - 1; ++i) {
 174                 rhs[i] = 4 * knots[i].y + 2 * knots[i + 1].y;
 175         }
 176         rhs[0] = knots[0].y + 2 * knots[1].y;
 177         rhs[n - 1] = (8 * knots[n - 1].y + knots[n].y) / 2.0;
 178
 179         // Get first control points Y-values
 180         double* y = solve (rhs);
 181
 182         for (Points::size_type i = 0; i < n; ++i) {
 183
 184                 firstControlPoints.push_back (Duple (x[i], y[i]));
 185
 186                 if (i < n - 1) {
 187                         secondControlPoints.push_back (Duple (2 * knots [i + 1].x - x[i + 1],
 188                                                               2 * knots[i + 1].y - y[i + 1]));
 189                 } else {
 190                         secondControlPoints.push_back (Duple ((knots [n].x + x[n - 1]) / 2,
 191                                                               (knots[n].y + y[n - 1]) / 2));
 192                 }
 193         }
 194
 195         delete [] x;
 196         delete [] y;
 197 }
 198
 199 /** Solves a tridiagonal system for one of coordinates (x or y)
 200  *  of first Bezier control points.
 201  */
 202
 203 double*
 204 Curve::solve (std::vector<double> const & rhs)
 205 {
 206         std::vector<double>::size_type n = rhs.size();
 207         double* x = new double[n]; // Solution vector.
 208         double* tmp = new double[n]; // Temp workspace.
 209
 210         double b = 2.0;
 211
 212         x[0] = rhs[0] / b;
 213
 214         for (std::vector<double>::size_type i = 1; i < n; i++) {
 215                 // Decomposition and forward substitution.
 216                 tmp[i] = 1 / b;
 217                 b = (i < n - 1 ? 4.0 : 3.5) - tmp[i];
 218                 x[i] = (rhs[i] - x[i - 1]) / b;
 219         }
 220
 221         for (std::vector<double>::size_type i = 1; i < n; i++) {
 222                 // Backsubstitution
 223                 x[n - i - 1] -= tmp[n - i] * x[n - i];
 224         }
 225
 226         delete [] tmp;
 227
 228         return x;
 229 }