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