agg/source/agg_bspline.cpp

   1 //----------------------------------------------------------------------------
   2 // Anti-Grain Geometry - Version 2.3
   3 // Copyright (C) 2002-2005 Maxim Shemanarev (http://www.antigrain.com)
   4 //
   5 // Permission to copy, use, modify, sell and distribute this software
   6 // is granted provided this copyright notice appears in all copies.
   7 // This software is provided "as is" without express or implied
   8 // warranty, and with no claim as to its suitability for any purpose.
   9 //
  10 //----------------------------------------------------------------------------
  11 // Contact: mcseem@antigrain.com
  12 //          mcseemagg@yahoo.com
  13 //          http://www.antigrain.com
  14 //----------------------------------------------------------------------------
  15 //
  16 // class bspline
  17 //
  18 //----------------------------------------------------------------------------
  19
  20
  21 #include "agg_bspline.h"
  22
  23 namespace agg
  24 {
  25
  26     //------------------------------------------------------------------------
  27     bspline::~bspline()
  28     {
  29         delete [] m_am;
  30     }
  31
  32
  33     //------------------------------------------------------------------------
  34     bspline::bspline() :
  35         m_max(0),
  36         m_num(0),
  37         m_x(0),
  38         m_y(0),
  39         m_am(0),
  40         m_last_idx(-1)
  41     {
  42     }
  43
  44     //------------------------------------------------------------------------
  45     bspline::bspline(int num) :
  46         m_max(0),
  47         m_num(0),
  48         m_x(0),
  49         m_y(0),
  50         m_am(0),
  51         m_last_idx(-1)
  52     {
  53         init(num);
  54     }
  55
  56     //------------------------------------------------------------------------
  57     bspline::bspline(int num, const double* x, const double* y) :
  58         m_max(0),
  59         m_num(0),
  60         m_x(0),
  61         m_y(0),
  62         m_am(0),
  63         m_last_idx(-1)
  64     {
  65         init(num, x, y);
  66     }
  67
  68
  69     //------------------------------------------------------------------------
  70     void bspline::init(int max)
  71     {
  72         if(max > 2 && max > m_max)
  73         {
  74             delete [] m_am;
  75             m_am = new double[max * 3];
  76             m_max = max;
  77             m_x = m_am + m_max;
  78             m_y = m_am + m_max * 2;
  79         }
  80         m_num = 0;
  81         m_last_idx = -1;
  82     }
  83
  84
  85     //------------------------------------------------------------------------
  86     void bspline::add_point(double x, double y)
  87     {
  88         if(m_num < m_max)
  89         {
  90             m_x[m_num] = x;
  91             m_y[m_num] = y;
  92             ++m_num;
  93         }
  94     }
  95
  96
  97     //------------------------------------------------------------------------
  98     void bspline::prepare()
  99     {
 100         if(m_num > 2)
 101         {
 102             int i, k, n1;
 103             double* temp;
 104             double* r;
 105             double* s;
 106             double* al;
 107             double h, p, d, f, e;
 108
 109             for(k = 0; k < m_num; k++)
 110             {
 111                 m_am[k] = 0.0;
 112             }
 113
 114             n1 = 3 * m_num;
 115
 116             al = new double[n1];
 117             temp = al;
 118
 119             for(k = 0; k < n1; k++)
 120             {
 121                 temp[k] = 0.0;
 122             }
 123
 124             r = temp + m_num;
 125             s = temp + m_num * 2;
 126
 127             n1 = m_num - 1;
 128             d = m_x[1] - m_x[0];
 129             e = (m_y[1] - m_y[0]) / d;
 130
 131             for(k = 1; k < n1; k++)
 132             {
 133                 h     = d;
 134                 d     = m_x[k + 1] - m_x[k];
 135                 f     = e;
 136                 e     = (m_y[k + 1] - m_y[k]) / d;
 137                 al[k] = d / (d + h);
 138                 r[k]  = 1.0 - al[k];
 139                 s[k]  = 6.0 * (e - f) / (h + d);
 140             }
 141
 142             for(k = 1; k < n1; k++)
 143             {
 144                 p = 1.0 / (r[k] * al[k - 1] + 2.0);
 145                 al[k] *= -p;
 146                 s[k] = (s[k] - r[k] * s[k - 1]) * p;
 147             }
 148
 149             m_am[n1]     = 0.0;
 150             al[n1 - 1]   = s[n1 - 1];
 151             m_am[n1 - 1] = al[n1 - 1];
 152
 153             for(k = n1 - 2, i = 0; i < m_num - 2; i++, k--)
 154             {
 155                 al[k]   = al[k] * al[k + 1] + s[k];
 156                 m_am[k] = al[k];
 157             }
 158             delete [] al;
 159         }
 160         m_last_idx = -1;
 161     }
 162
 163
 164
 165     //------------------------------------------------------------------------
 166     void bspline::init(int num, const double* x, const double* y)
 167     {
 168         if(num > 2)
 169         {
 170             init(num);
 171             int i;
 172             for(i = 0; i < num; i++)
 173             {
 174                 add_point(*x++, *y++);
 175             }
 176             prepare();
 177         }
 178         m_last_idx = -1;
 179     }
 180
 181
 182     //------------------------------------------------------------------------
 183     void bspline::bsearch(int n, const double *x, double x0, int *i)
 184     {
 185         int j = n - 1;
 186         int k;
 187
 188         for(*i = 0; (j - *i) > 1; )
 189         {
 190             if(x0 < x[k = (*i + j) >> 1]) j = k;
 191             else                         *i = k;
 192         }
 193     }
 194
 195
 196
 197     //------------------------------------------------------------------------
 198     double bspline::interpolation(double x, int i) const
 199     {
 200         int j = i + 1;
 201         double d = m_x[i] - m_x[j];
 202         double h = x - m_x[j];
 203         double r = m_x[i] - x;
 204         double p = d * d / 6.0;
 205         return (m_am[j] * r * r * r + m_am[i] * h * h * h) / 6.0 / d +
 206                ((m_y[j] - m_am[j] * p) * r + (m_y[i] - m_am[i] * p) * h) / d;
 207     }
 208
 209
 210     //------------------------------------------------------------------------
 211     double bspline::extrapolation_left(double x) const
 212     {
 213         double d = m_x[1] - m_x[0];
 214         return (-d * m_am[1] / 6 + (m_y[1] - m_y[0]) / d) *
 215                (x - m_x[0]) +
 216                m_y[0];
 217     }
 218
 219     //------------------------------------------------------------------------
 220     double bspline::extrapolation_right(double x) const
 221     {
 222         double d = m_x[m_num - 1] - m_x[m_num - 2];
 223         return (d * m_am[m_num - 2] / 6 + (m_y[m_num - 1] - m_y[m_num - 2]) / d) *
 224                (x - m_x[m_num - 1]) +
 225                m_y[m_num - 1];
 226     }
 227
 228     //------------------------------------------------------------------------
 229     double bspline::get(double x) const
 230     {
 231         if(m_num > 2)
 232         {
 233             int i;
 234
 235             // Extrapolation on the left
 236             if(x < m_x[0]) return extrapolation_left(x);
 237
 238             // Extrapolation on the right
 239             if(x >= m_x[m_num - 1]) return extrapolation_right(x);
 240
 241             // Interpolation
 242             bsearch(m_num, m_x, x, &i);
 243             return interpolation(x, i);
 244         }
 245         return 0.0;
 246     }
 247
 248
 249     //------------------------------------------------------------------------
 250     double bspline::get_stateful(double x) const
 251     {
 252         if(m_num > 2)
 253         {
 254             // Extrapolation on the left
 255             if(x < m_x[0]) return extrapolation_left(x);
 256
 257             // Extrapolation on the right
 258             if(x >= m_x[m_num - 1]) return extrapolation_right(x);
 259
 260             if(m_last_idx >= 0)
 261             {
 262                 // Check if x is not in current range
 263                 if(x < m_x[m_last_idx] || x > m_x[m_last_idx + 1])
 264                 {
 265                     // Check if x between next points (most probably)
 266                     if(m_last_idx < m_num - 2 &&
 267                        x >= m_x[m_last_idx + 1] &&
 268                        x <= m_x[m_last_idx + 2])
 269                     {
 270                         ++m_last_idx;
 271                     }
 272                     else
 273                     if(m_last_idx > 0 &&
 274                        x >= m_x[m_last_idx - 1] &&
 275                        x <= m_x[m_last_idx])
 276                     {
 277                         // x is between pevious points
 278                         --m_last_idx;
 279                     }
 280                     else
 281                     {
 282                         // Else perform full search
 283                         bsearch(m_num, m_x, x, &m_last_idx);
 284                     }
 285                 }
 286                 return interpolation(x, m_last_idx);
 287             }
 288             else
 289             {
 290                 // Interpolation
 291                 bsearch(m_num, m_x, x, &m_last_idx);
 292                 return interpolation(x, m_last_idx);
 293             }
 294         }
 295         return 0.0;
 296     }
 297
 298 }
 299