1 //---------------------------------------------------------------------------- 2 // Anti-Grain Geometry - Version 2.4 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 #include "agg_bspline.h" 21 22 namespace agg 23 { 24 //------------------------------------------------------------------------ bspline()25 bspline::bspline() : 26 m_max(0), 27 m_num(0), 28 m_x(0), 29 m_y(0), 30 m_last_idx(-1) 31 { 32 } 33 34 //------------------------------------------------------------------------ bspline(int num)35 bspline::bspline(int num) : 36 m_max(0), 37 m_num(0), 38 m_x(0), 39 m_y(0), 40 m_last_idx(-1) 41 { 42 init(num); 43 } 44 45 //------------------------------------------------------------------------ bspline(int num,const double * x,const double * y)46 bspline::bspline(int num, const double* x, const double* y) : 47 m_max(0), 48 m_num(0), 49 m_x(0), 50 m_y(0), 51 m_last_idx(-1) 52 { 53 init(num, x, y); 54 } 55 56 57 //------------------------------------------------------------------------ init(int max)58 void bspline::init(int max) 59 { 60 if(max > 2 && max > m_max) 61 { 62 m_am.resize(max * 3); 63 m_max = max; 64 m_x = &m_am[m_max]; 65 m_y = &m_am[m_max * 2]; 66 } 67 m_num = 0; 68 m_last_idx = -1; 69 } 70 71 72 //------------------------------------------------------------------------ add_point(double x,double y)73 void bspline::add_point(double x, double y) 74 { 75 if(m_num < m_max) 76 { 77 m_x[m_num] = x; 78 m_y[m_num] = y; 79 ++m_num; 80 } 81 } 82 83 84 //------------------------------------------------------------------------ prepare()85 void bspline::prepare() 86 { 87 if(m_num > 2) 88 { 89 int i, k, n1; 90 double* temp; 91 double* r; 92 double* s; 93 double h, p, d, f, e; 94 95 for(k = 0; k < m_num; k++) 96 { 97 m_am[k] = 0.0; 98 } 99 100 n1 = 3 * m_num; 101 102 pod_array<double> al(n1); 103 temp = &al[0]; 104 105 for(k = 0; k < n1; k++) 106 { 107 temp[k] = 0.0; 108 } 109 110 r = temp + m_num; 111 s = temp + m_num * 2; 112 113 n1 = m_num - 1; 114 d = m_x[1] - m_x[0]; 115 e = (m_y[1] - m_y[0]) / d; 116 117 for(k = 1; k < n1; k++) 118 { 119 h = d; 120 d = m_x[k + 1] - m_x[k]; 121 f = e; 122 e = (m_y[k + 1] - m_y[k]) / d; 123 al[k] = d / (d + h); 124 r[k] = 1.0 - al[k]; 125 s[k] = 6.0 * (e - f) / (h + d); 126 } 127 128 for(k = 1; k < n1; k++) 129 { 130 p = 1.0 / (r[k] * al[k - 1] + 2.0); 131 al[k] *= -p; 132 s[k] = (s[k] - r[k] * s[k - 1]) * p; 133 } 134 135 m_am[n1] = 0.0; 136 al[n1 - 1] = s[n1 - 1]; 137 m_am[n1 - 1] = al[n1 - 1]; 138 139 for(k = n1 - 2, i = 0; i < m_num - 2; i++, k--) 140 { 141 al[k] = al[k] * al[k + 1] + s[k]; 142 m_am[k] = al[k]; 143 } 144 } 145 m_last_idx = -1; 146 } 147 148 149 150 //------------------------------------------------------------------------ init(int num,const double * x,const double * y)151 void bspline::init(int num, const double* x, const double* y) 152 { 153 if(num > 2) 154 { 155 init(num); 156 int i; 157 for(i = 0; i < num; i++) 158 { 159 add_point(*x++, *y++); 160 } 161 prepare(); 162 } 163 m_last_idx = -1; 164 } 165 166 167 //------------------------------------------------------------------------ bsearch(int n,const double * x,double x0,int * i)168 void bspline::bsearch(int n, const double *x, double x0, int *i) 169 { 170 int j = n - 1; 171 int k; 172 173 for(*i = 0; (j - *i) > 1; ) 174 { 175 if(x0 < x[k = (*i + j) >> 1]) j = k; 176 else *i = k; 177 } 178 } 179 180 181 182 //------------------------------------------------------------------------ interpolation(double x,int i) const183 double bspline::interpolation(double x, int i) const 184 { 185 int j = i + 1; 186 double d = m_x[i] - m_x[j]; 187 double h = x - m_x[j]; 188 double r = m_x[i] - x; 189 double p = d * d / 6.0; 190 return (m_am[j] * r * r * r + m_am[i] * h * h * h) / 6.0 / d + 191 ((m_y[j] - m_am[j] * p) * r + (m_y[i] - m_am[i] * p) * h) / d; 192 } 193 194 195 //------------------------------------------------------------------------ extrapolation_left(double x) const196 double bspline::extrapolation_left(double x) const 197 { 198 double d = m_x[1] - m_x[0]; 199 return (-d * m_am[1] / 6 + (m_y[1] - m_y[0]) / d) * 200 (x - m_x[0]) + 201 m_y[0]; 202 } 203 204 //------------------------------------------------------------------------ extrapolation_right(double x) const205 double bspline::extrapolation_right(double x) const 206 { 207 double d = m_x[m_num - 1] - m_x[m_num - 2]; 208 return (d * m_am[m_num - 2] / 6 + (m_y[m_num - 1] - m_y[m_num - 2]) / d) * 209 (x - m_x[m_num - 1]) + 210 m_y[m_num - 1]; 211 } 212 213 //------------------------------------------------------------------------ get(double x) const214 double bspline::get(double x) const 215 { 216 if(m_num > 2) 217 { 218 int i; 219 220 // Extrapolation on the left 221 if(x < m_x[0]) return extrapolation_left(x); 222 223 // Extrapolation on the right 224 if(x >= m_x[m_num - 1]) return extrapolation_right(x); 225 226 // Interpolation 227 bsearch(m_num, m_x, x, &i); 228 return interpolation(x, i); 229 } 230 return 0.0; 231 } 232 233 234 //------------------------------------------------------------------------ get_stateful(double x) const235 double bspline::get_stateful(double x) const 236 { 237 if(m_num > 2) 238 { 239 // Extrapolation on the left 240 if(x < m_x[0]) return extrapolation_left(x); 241 242 // Extrapolation on the right 243 if(x >= m_x[m_num - 1]) return extrapolation_right(x); 244 245 if(m_last_idx >= 0) 246 { 247 // Check if x is not in current range 248 if(x < m_x[m_last_idx] || x > m_x[m_last_idx + 1]) 249 { 250 // Check if x between next points (most probably) 251 if(m_last_idx < m_num - 2 && 252 x >= m_x[m_last_idx + 1] && 253 x <= m_x[m_last_idx + 2]) 254 { 255 ++m_last_idx; 256 } 257 else 258 if(m_last_idx > 0 && 259 x >= m_x[m_last_idx - 1] && 260 x <= m_x[m_last_idx]) 261 { 262 // x is between pevious points 263 --m_last_idx; 264 } 265 else 266 { 267 // Else perform full search 268 bsearch(m_num, m_x, x, &m_last_idx); 269 } 270 } 271 return interpolation(x, m_last_idx); 272 } 273 else 274 { 275 // Interpolation 276 bsearch(m_num, m_x, x, &m_last_idx); 277 return interpolation(x, m_last_idx); 278 } 279 } 280 return 0.0; 281 } 282 283 } 284 285