20 #include "agg_bspline.h" 35 bspline::bspline(
int num) :
46 bspline::bspline(
int num,
const double* x,
const double* y) :
58 void bspline::init(
int max)
60 if(max > 2 && max > m_max)
65 m_y = &m_am[m_max * 2];
73 void bspline::add_point(
double x,
double y)
85 void bspline::prepare()
95 for(k = 0; k < m_num; k++)
102 pod_array<double> al(n1);
105 for(k = 0; k < n1; k++)
111 s = temp + m_num * 2;
115 e = (m_y[1] - m_y[0]) / d;
117 for(k = 1; k < n1; k++)
120 d = m_x[k + 1] - m_x[k];
122 e = (m_y[k + 1] - m_y[k]) / d;
125 s[k] = 6.0 * (e - f) / (h + d);
128 for(k = 1; k < n1; k++)
130 p = 1.0 / (r[k] * al[k - 1] + 2.0);
132 s[k] = (s[k] - r[k] * s[k - 1]) * p;
136 al[n1 - 1] = s[n1 - 1];
137 m_am[n1 - 1] = al[n1 - 1];
139 for(k = n1 - 2, i = 0; i < m_num - 2; i++, k--)
141 al[k] = al[k] * al[k + 1] + s[k];
151 void bspline::init(
int num,
const double* x,
const double* y)
157 for(i = 0; i < num; i++)
159 add_point(*x++, *y++);
168 void bspline::bsearch(
int n,
const double *x,
double x0,
int *i)
173 for(*i = 0; (j - *i) > 1; )
175 if(x0 < x[k = (*i + j) >> 1]) j = k;
183 double bspline::interpolation(
double x,
int i)
const 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;
196 double bspline::extrapolation_left(
double x)
const 198 double d = m_x[1] - m_x[0];
199 return (-d * m_am[1] / 6 + (m_y[1] - m_y[0]) / d) *
205 double bspline::extrapolation_right(
double x)
const 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]) +
214 double bspline::get(
double x)
const 221 if(x < m_x[0])
return extrapolation_left(x);
224 if(x >= m_x[m_num - 1])
return extrapolation_right(x);
227 bsearch(m_num, m_x, x, &i);
228 return interpolation(x, i);
235 double bspline::get_stateful(
double x)
const 240 if(x < m_x[0])
return extrapolation_left(x);
243 if(x >= m_x[m_num - 1])
return extrapolation_right(x);
248 if(x < m_x[m_last_idx] || x > m_x[m_last_idx + 1])
251 if(m_last_idx < m_num - 2 &&
252 x >= m_x[m_last_idx + 1] &&
253 x <= m_x[m_last_idx + 2])
259 x >= m_x[m_last_idx - 1] &&
260 x <= m_x[m_last_idx])
268 bsearch(m_num, m_x, x, &m_last_idx);
271 return interpolation(x, m_last_idx);
276 bsearch(m_num, m_x, x, &m_last_idx);
277 return interpolation(x, m_last_idx);