1 //---------------------------------------------------------------------------- 2 // Anti-Grain Geometry (AGG) - Version 2.5 3 // A high quality rendering engine for C++ 4 // Copyright (C) 2002-2006 Maxim Shemanarev 5 // Contact: mcseem@antigrain.com 6 // mcseemagg@yahoo.com 7 // http://antigrain.com 8 // 9 // AGG is free software; you can redistribute it and/or 10 // modify it under the terms of the GNU General Public License 11 // as published by the Free Software Foundation; either version 2 12 // of the License, or (at your option) any later version. 13 // 14 // AGG is distributed in the hope that it will be useful, 15 // but WITHOUT ANY WARRANTY; without even the implied warranty of 16 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 17 // GNU General Public License for more details. 18 // 19 // You should have received a copy of the GNU General Public License 20 // along with AGG; if not, write to the Free Software 21 // Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, 22 // MA 02110-1301, USA. 23 //---------------------------------------------------------------------------- 24 25 #include "agg_bspline.h" 26 27 namespace agg 28 { 29 //------------------------------------------------------------------------ bspline()30 bspline::bspline() : 31 m_max(0), 32 m_num(0), 33 m_x(0), 34 m_y(0), 35 m_last_idx(-1) 36 { 37 } 38 39 //------------------------------------------------------------------------ bspline(int num)40 bspline::bspline(int num) : 41 m_max(0), 42 m_num(0), 43 m_x(0), 44 m_y(0), 45 m_last_idx(-1) 46 { 47 init(num); 48 } 49 50 //------------------------------------------------------------------------ bspline(int num,const double * x,const double * y)51 bspline::bspline(int num, const double* x, const double* y) : 52 m_max(0), 53 m_num(0), 54 m_x(0), 55 m_y(0), 56 m_last_idx(-1) 57 { 58 init(num, x, y); 59 } 60 61 62 //------------------------------------------------------------------------ init(int max)63 void bspline::init(int max) 64 { 65 if(max > 2 && max > m_max) 66 { 67 m_am.resize(max * 3); 68 m_max = max; 69 m_x = &m_am[m_max]; 70 m_y = &m_am[m_max * 2]; 71 } 72 m_num = 0; 73 m_last_idx = -1; 74 } 75 76 77 //------------------------------------------------------------------------ add_point(double x,double y)78 void bspline::add_point(double x, double y) 79 { 80 if(m_num < m_max) 81 { 82 m_x[m_num] = x; 83 m_y[m_num] = y; 84 ++m_num; 85 } 86 } 87 88 89 //------------------------------------------------------------------------ prepare()90 void bspline::prepare() 91 { 92 if(m_num > 2) 93 { 94 int i, k, n1; 95 double* temp; 96 double* r; 97 double* s; 98 double h, p, d, f, e; 99 100 for(k = 0; k < m_num; k++) 101 { 102 m_am[k] = 0.0; 103 } 104 105 n1 = 3 * m_num; 106 107 pod_array<double> al(n1); 108 temp = &al[0]; 109 110 for(k = 0; k < n1; k++) 111 { 112 temp[k] = 0.0; 113 } 114 115 r = temp + m_num; 116 s = temp + m_num * 2; 117 118 n1 = m_num - 1; 119 d = m_x[1] - m_x[0]; 120 e = (m_y[1] - m_y[0]) / d; 121 122 for(k = 1; k < n1; k++) 123 { 124 h = d; 125 d = m_x[k + 1] - m_x[k]; 126 f = e; 127 e = (m_y[k + 1] - m_y[k]) / d; 128 al[k] = d / (d + h); 129 r[k] = 1.0 - al[k]; 130 s[k] = 6.0 * (e - f) / (h + d); 131 } 132 133 for(k = 1; k < n1; k++) 134 { 135 p = 1.0 / (r[k] * al[k - 1] + 2.0); 136 al[k] *= -p; 137 s[k] = (s[k] - r[k] * s[k - 1]) * p; 138 } 139 140 m_am[n1] = 0.0; 141 al[n1 - 1] = s[n1 - 1]; 142 m_am[n1 - 1] = al[n1 - 1]; 143 144 for(k = n1 - 2, i = 0; i < m_num - 2; i++, k--) 145 { 146 al[k] = al[k] * al[k + 1] + s[k]; 147 m_am[k] = al[k]; 148 } 149 } 150 m_last_idx = -1; 151 } 152 153 154 155 //------------------------------------------------------------------------ init(int num,const double * x,const double * y)156 void bspline::init(int num, const double* x, const double* y) 157 { 158 if(num > 2) 159 { 160 init(num); 161 int i; 162 for(i = 0; i < num; i++) 163 { 164 add_point(*x++, *y++); 165 } 166 prepare(); 167 } 168 m_last_idx = -1; 169 } 170 171 172 //------------------------------------------------------------------------ bsearch(int n,const double * x,double x0,int * i)173 void bspline::bsearch(int n, const double *x, double x0, int *i) 174 { 175 int j = n - 1; 176 int k; 177 178 for(*i = 0; (j - *i) > 1; ) 179 { 180 if(x0 < x[k = (*i + j) >> 1]) j = k; 181 else *i = k; 182 } 183 } 184 185 186 187 //------------------------------------------------------------------------ interpolation(double x,int i) const188 double bspline::interpolation(double x, int i) const 189 { 190 int j = i + 1; 191 double d = m_x[i] - m_x[j]; 192 double h = x - m_x[j]; 193 double r = m_x[i] - x; 194 double p = d * d / 6.0; 195 return (m_am[j] * r * r * r + m_am[i] * h * h * h) / 6.0 / d + 196 ((m_y[j] - m_am[j] * p) * r + (m_y[i] - m_am[i] * p) * h) / d; 197 } 198 199 200 //------------------------------------------------------------------------ extrapolation_left(double x) const201 double bspline::extrapolation_left(double x) const 202 { 203 double d = m_x[1] - m_x[0]; 204 return (-d * m_am[1] / 6 + (m_y[1] - m_y[0]) / d) * 205 (x - m_x[0]) + 206 m_y[0]; 207 } 208 209 //------------------------------------------------------------------------ extrapolation_right(double x) const210 double bspline::extrapolation_right(double x) const 211 { 212 double d = m_x[m_num - 1] - m_x[m_num - 2]; 213 return (d * m_am[m_num - 2] / 6 + (m_y[m_num - 1] - m_y[m_num - 2]) / d) * 214 (x - m_x[m_num - 1]) + 215 m_y[m_num - 1]; 216 } 217 218 //------------------------------------------------------------------------ get(double x) const219 double bspline::get(double x) const 220 { 221 if(m_num > 2) 222 { 223 int i; 224 225 // Extrapolation on the left 226 if(x < m_x[0]) return extrapolation_left(x); 227 228 // Extrapolation on the right 229 if(x >= m_x[m_num - 1]) return extrapolation_right(x); 230 231 // Interpolation 232 bsearch(m_num, m_x, x, &i); 233 return interpolation(x, i); 234 } 235 return 0.0; 236 } 237 238 239 //------------------------------------------------------------------------ get_stateful(double x) const240 double bspline::get_stateful(double x) const 241 { 242 if(m_num > 2) 243 { 244 // Extrapolation on the left 245 if(x < m_x[0]) return extrapolation_left(x); 246 247 // Extrapolation on the right 248 if(x >= m_x[m_num - 1]) return extrapolation_right(x); 249 250 if(m_last_idx >= 0) 251 { 252 // Check if x is not in current range 253 if(x < m_x[m_last_idx] || x > m_x[m_last_idx + 1]) 254 { 255 // Check if x between next points (most probably) 256 if(m_last_idx < m_num - 2 && 257 x >= m_x[m_last_idx + 1] && 258 x <= m_x[m_last_idx + 2]) 259 { 260 ++m_last_idx; 261 } 262 else 263 if(m_last_idx > 0 && 264 x >= m_x[m_last_idx - 1] && 265 x <= m_x[m_last_idx]) 266 { 267 // x is between pevious points 268 --m_last_idx; 269 } 270 else 271 { 272 // Else perform full search 273 bsearch(m_num, m_x, x, &m_last_idx); 274 } 275 } 276 return interpolation(x, m_last_idx); 277 } 278 else 279 { 280 // Interpolation 281 bsearch(m_num, m_x, x, &m_last_idx); 282 return interpolation(x, m_last_idx); 283 } 284 } 285 return 0.0; 286 } 287 288 } 289 290