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