libs/canvas/canvas/interpolated_curve.h

   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 #ifndef __CANVAS_INTERPOLATED_CURVE_H__
  20 #define __CANVAS_INTERPOLATED_CURVE_H__
  21
  22 #include "canvas/visibility.h"
  23 #include "canvas/types.h"
  24
  25 namespace ArdourCanvas {
  26
  27 class LIBCANVAS_API InterpolatedCurve
  28 {
  29 public:
  30     enum SplineType {
  31             CatmullRomUniform,
  32             CatmullRomCentripetal,
  33     };
  34
  35 protected:
  36
  37         /**
  38          * This method will calculate the Catmull-Rom interpolation curve, returning
  39          * it as a list of Coord coordinate objects.  This method in particular
  40          * adds the first and last control points which are not visible, but required
  41          * for calculating the spline.
  42          *
  43          * @param coordinates The list of original straight line points to calculate
  44          * an interpolation from.
  45          * @param points_per_segment The integer number of equally spaced points to
  46          * return along each curve.  The actual distance between each
  47          * point will depend on the spacing between the control points.
  48          * @return The list of interpolated coordinates.
  49          * @param curve_type Chordal (stiff), Uniform(floppy), or Centripetal(medium)
  50          * @throws gov.ca.water.shapelite.analysis.CatmullRomException if
  51          * points_per_segment is less than 2.
  52          */
  53         static void
  54                 interpolate (const Points& coordinates, uint32_t points_per_segment, SplineType curve_type, bool closed, Points& results)
  55         {
  56                 if (points_per_segment < 2) {
  57                         return;
  58                 }
  59
  60                 // Cannot interpolate curves given only two points.  Two points
  61                 // is best represented as a simple line segment.
  62                 if (coordinates.size() < 3) {
  63                         results = coordinates;
  64                         return;
  65                 }
  66
  67                 // Copy the incoming coordinates. We need to modify it during interpolation
  68                 Points vertices = coordinates;
  69
  70                 // Test whether the shape is open or closed by checking to see if
  71                 // the first point intersects with the last point.  M and Z are ignored.
  72                 if (closed) {
  73                         // Use the second and second from last points as control points.
  74                         // get the second point.
  75                         Duple p2 = vertices[1];
  76                         // get the point before the last point
  77                         Duple pn1 = vertices[vertices.size() - 2];
  78
  79                         // insert the second from the last point as the first point in the list
  80                         // because when the shape is closed it keeps wrapping around to
  81                         // the second point.
  82                         vertices.insert(vertices.begin(), pn1);
  83                         // add the second point to the end.
  84                         vertices.push_back(p2);
  85                 } else {
  86                         // The shape is open, so use control points that simply extend
  87                         // the first and last segments
  88
  89                         // Get the change in x and y between the first and second coordinates.
  90                         double dx = vertices[1].x - vertices[0].x;
  91                         double dy = vertices[1].y - vertices[0].y;
  92
  93                         // Then using the change, extrapolate backwards to find a control point.
  94                         double x1 = vertices[0].x - dx;
  95                         double y1 = vertices[0].y - dy;
  96
  97                         // Actaully create the start point from the extrapolated values.
  98                         Duple start (x1, y1);
  99
 100                         // Repeat for the end control point.
 101                         int n = vertices.size() - 1;
 102                         dx = vertices[n].x - vertices[n - 1].x;
 103                         dy = vertices[n].y - vertices[n - 1].y;
 104                         double xn = vertices[n].x + dx;
 105                         double yn = vertices[n].y + dy;
 106                         Duple end (xn, yn);
 107
 108                         // insert the start control point at the start of the vertices list.
 109                         vertices.insert (vertices.begin(), start);
 110
 111                         // append the end control ponit to the end of the vertices list.
 112                         vertices.push_back (end);
 113                 }
 114
 115                 // When looping, remember that each cycle requires 4 points, starting
 116                 // with i and ending with i+3.  So we don't loop through all the points.
 117
 118                 for (Points::size_type i = 0; i < vertices.size() - 3; i++) {
 119
 120                         // Actually calculate the Catmull-Rom curve for one segment.
 121                         Points r;
 122
 123                         _interpolate (vertices, i, points_per_segment, curve_type, r);
 124
 125                         // Since the middle points are added twice, once for each bordering
 126                         // segment, we only add the 0 index result point for the first
 127                         // segment.  Otherwise we will have duplicate points.
 128
 129                         if (results.size() > 0) {
 130                                 r.erase (r.begin());
 131                         }
 132
 133                         // Add the coordinates for the segment to the result list.
 134
 135                         results.insert (results.end(), r.begin(), r.end());
 136                 }
 137         }
 138
 139 private:
 140         /**
 141          * Calculate the same values but introduces the ability to "parameterize" the t
 142          * values used in the calculation. This is based on Figure 3 from
 143          * http://www.cemyuksel.com/research/catmullrom_param/catmullrom.pdf
 144          *
 145          * @param p An array of double values of length 4, where interpolation
 146          * occurs from p1 to p2.
 147          * @param time An array of time measures of length 4, corresponding to each
 148          * p value.
 149          * @param t the actual interpolation ratio from 0 to 1 representing the
 150          * position between p1 and p2 to interpolate the value.
 151          */
 152         static double
 153                 __interpolate (double p[4], double time[4], double t)
 154         {
 155                 const double L01 = p[0] * (time[1] - t) / (time[1] - time[0]) + p[1] * (t - time[0]) / (time[1] - time[0]);
 156                 const double L12 = p[1] * (time[2] - t) / (time[2] - time[1]) + p[2] * (t - time[1]) / (time[2] - time[1]);
 157                 const double L23 = p[2] * (time[3] - t) / (time[3] - time[2]) + p[3] * (t - time[2]) / (time[3] - time[2]);
 158                 const double L012 = L01 * (time[2] - t) / (time[2] - time[0]) + L12 * (t - time[0]) / (time[2] - time[0]);
 159                 const double L123 = L12 * (time[3] - t) / (time[3] - time[1]) + L23 * (t - time[1]) / (time[3] - time[1]);
 160                 const double C12 = L012 * (time[2] - t) / (time[2] - time[1]) + L123 * (t - time[1]) / (time[2] - time[1]);
 161                 return C12;
 162         }
 163
 164         /**
 165          * Given a list of control points, this will create a list of points_per_segment
 166          * points spaced uniformly along the resulting Catmull-Rom curve.
 167          *
 168          * @param points The list of control points, leading and ending with a
 169          * coordinate that is only used for controling the spline and is not visualized.
 170          * @param index The index of control point p0, where p0, p1, p2, and p3 are
 171          * used in order to create a curve between p1 and p2.
 172          * @param points_per_segment The total number of uniformly spaced interpolated
 173          * points to calculate for each segment. The larger this number, the
 174          * smoother the resulting curve.
 175          * @param curve_type Clarifies whether the curve should use uniform, chordal
 176          * or centripetal curve types. Uniform can produce loops, chordal can
 177          * produce large distortions from the original lines, and centripetal is an
 178          * optimal balance without spaces.
 179          * @return the list of coordinates that define the CatmullRom curve
 180          * between the points defined by index+1 and index+2.
 181          */
 182         static void
 183                 _interpolate (const Points& points, Points::size_type index, int points_per_segment, SplineType curve_type, Points& results)
 184         {
 185                 double x[4];
 186                 double y[4];
 187                 double time[4];
 188
 189                 for (int i = 0; i < 4; i++) {
 190                         x[i] = points[index + i].x;
 191                         y[i] = points[index + i].y;
 192                         time[i] = i;
 193                 }
 194
 195                 double tstart = 1;
 196                 double tend = 2;
 197
 198                 if (curve_type != CatmullRomUniform) {
 199                         double total = 0;
 200                         for (int i = 1; i < 4; i++) {
 201                                 double dx = x[i] - x[i - 1];
 202                                 double dy = y[i] - y[i - 1];
 203                                 if (curve_type == CatmullRomCentripetal) {
 204                                         total += pow (dx * dx + dy * dy, .25);
 205                                 } else {
 206                                         total += pow (dx * dx + dy * dy, .5);
 207                                 }
 208                                 time[i] = total;
 209                         }
 210                         tstart = time[1];
 211                         tend = time[2];
 212                 }
 213
 214                 int segments = points_per_segment - 1;
 215                 results.push_back (points[index + 1]);
 216
 217                 for (int i = 1; i < segments; i++) {
 218                         double xi = __interpolate (x, time, tstart + (i * (tend - tstart)) / segments);
 219                         double yi = __interpolate (y, time, tstart + (i * (tend - tstart)) / segments);
 220                         results.push_back (Duple (xi, yi));
 221                 }
 222
 223                 results.push_back (points[index + 2]);
 224         }
 225 };
 226
 227 }
 228
 229 #endif