Anti-Grain Geometry Tutorial
agg_bspline.cpp
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  //------------------------------------------------------------------------
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 
Definition: agg_arc.cpp:24