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
22namespace agg
23{
24    //------------------------------------------------------------------------
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    //------------------------------------------------------------------------
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    //------------------------------------------------------------------------
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    //------------------------------------------------------------------------
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    //------------------------------------------------------------------------
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    //------------------------------------------------------------------------
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    //------------------------------------------------------------------------
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    //------------------------------------------------------------------------
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    //------------------------------------------------------------------------
183    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    //------------------------------------------------------------------------
196    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    //------------------------------------------------------------------------
205    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    //------------------------------------------------------------------------
214    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    //------------------------------------------------------------------------
235    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