Line data Source code
1 : /******************************************************************************
2 : *
3 : * Project: DXF Translator
4 : * Purpose: Low level spline interpolation.
5 : * Author: David F. Rogers
6 : *
7 : ******************************************************************************
8 :
9 : This code is derived from the code associated with the book "An Introduction
10 : to NURBS" by David F. Rogers. More information on the book and the code is
11 : available at:
12 :
13 : http://www.nar-associates.com/nurbs/
14 :
15 : David F. Rogers agreed explicitly that the subset of his code in this file to
16 : be used under the following license (cf
17 : https://github.com/OSGeo/gdal/issues/6524)
18 :
19 : Copyright (c) 2009, David F. Rogers
20 : All rights reserved.
21 :
22 : Redistribution and use in source and binary forms, with or without
23 : modification, are permitted provided that the following conditions are met:
24 :
25 : * Redistributions of source code must retain the above copyright notice,
26 : this list of conditions and the following disclaimer.
27 : * Redistributions in binary form must reproduce the above copyright notice,
28 : this list of conditions and the following disclaimer in the documentation
29 : and/or other materials provided with the distribution.
30 : * Neither the name of David F. Rogers nor the names of its contributors
31 : may be used to endorse or promote products derived from this software
32 : without specific prior written permission.
33 :
34 : THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
35 : AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
36 : IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
37 : ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
38 : LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
39 : CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
40 : SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
41 : INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
42 : CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
43 : ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
44 : POSSIBILITY OF SUCH DAMAGE.
45 : */
46 :
47 : #include <stdio.h>
48 : #include <vector>
49 : #include "cpl_port.h"
50 :
51 : /* used by ogrdxflayer.cpp */
52 : void rbspline2(int npts, int k, int p1, double b[], double h[],
53 : bool bCalculateKnots, double x[], double p[]);
54 :
55 : /* used by ogrdxf_leader.cpp */
56 : void basis(int c, double t, int npts, double x[], double N[]);
57 :
58 : /* used by DWG driver */
59 : void rbspline(int npts, int k, int p1, double b[], double h[], double p[]);
60 : void rbsplinu(int npts, int k, int p1, double b[], double h[], double p[]);
61 :
62 : /************************************************************************/
63 : /* knotu() */
64 : /************************************************************************/
65 :
66 : /* Subroutine to generate a B-spline uniform (periodic) knot vector.
67 :
68 : c = order of the basis function
69 : n = the number of defining polygon vertices
70 : nplus2 = index of x() for the first occurrence of the maximum knot
71 : vector value nplusc = maximum value of the knot vector -- $n + c$ x[]
72 : = array containing the knot vector
73 : */
74 :
75 0 : static void knotu(int n, int c, double x[])
76 :
77 : {
78 : int nplusc, /* nplus2, */ i;
79 :
80 0 : nplusc = n + c;
81 : /* nplus2 = n + 2; */
82 :
83 0 : x[1] = 0;
84 0 : for (i = 2; i <= nplusc; i++)
85 : {
86 0 : x[i] = i - 1;
87 : }
88 0 : }
89 :
90 : /************************************************************************/
91 : /* knot() */
92 : /************************************************************************/
93 :
94 : /*
95 : Subroutine to generate a B-spline open knot vector with multiplicity
96 : equal to the order at the ends.
97 :
98 : c = order of the basis function
99 : n = the number of defining polygon vertices
100 : nplus2 = index of x() for the first occurrence of the maximum knot
101 : vector value nplusc = maximum value of the knot vector -- $n + c$ x()
102 : = array containing the knot vector
103 : */
104 :
105 1 : static void knot(int n, int c, double x[])
106 :
107 : {
108 1 : const int nplusc = n + c;
109 1 : const int nplus2 = n + 2;
110 :
111 1 : x[1] = 0.0;
112 11 : for (int i = 2; i <= nplusc; i++)
113 : {
114 10 : if ((i > c) && (i < nplus2))
115 4 : x[i] = x[i - 1] + 1.0;
116 : else
117 6 : x[i] = x[i - 1];
118 : }
119 1 : }
120 :
121 : /************************************************************************/
122 : /* basis() */
123 : /************************************************************************/
124 :
125 : /* Subroutine to generate B-spline basis functions--open knot vector
126 :
127 : C code for An Introduction to NURBS
128 : by David F. Rogers. Copyright (C) 2000 David F. Rogers,
129 : All rights reserved.
130 :
131 : Name: rbasis (split into rbasis and basis by AT)
132 : Language: C
133 : Subroutines called: none
134 : Book reference: Chapter 4, Sec. 4. , p 296
135 :
136 : c = order of the B-spline basis function
137 : d = first term of the basis function recursion relation
138 : e = second term of the basis function recursion relation
139 : h[] = array containing the homogeneous weights
140 : npts = number of defining polygon vertices
141 : nplusc = constant -- npts + c -- maximum number of knot values
142 : r[] = array containing the rationalbasis functions
143 : r[1] contains the basis function associated with B1 etc.
144 : t = parameter value
145 : N[] = array output from bbasis containing the values of the
146 : basis functions; should have capacity npts + c + 1
147 : x[] = knot vector
148 : */
149 :
150 1322 : void basis(int c, double t, int npts, double x[], double N[])
151 :
152 : {
153 1322 : const int nplusc = npts + c;
154 :
155 : /* calculate the first order nonrational basis functions n[i] */
156 :
157 13082 : for (int i = 1; i <= nplusc - 1; i++)
158 : {
159 11760 : if ((t >= x[i]) && (t < x[i + 1]))
160 1292 : N[i] = 1.0;
161 : else
162 10468 : N[i] = 0.0;
163 : }
164 :
165 : /* calculate the higher order nonrational basis functions */
166 :
167 5288 : for (int k = 2; k <= c; k++)
168 : {
169 31314 : for (int i = 1; i <= nplusc - k; i++)
170 : {
171 27348 : double d = 0.0;
172 27348 : double e = 0.0;
173 27348 : if (N[i] != 0) /* if the lower order basis function is zero skip the
174 : calculation */
175 : {
176 7656 : double denom = x[i + k - 1] - x[i];
177 7656 : if (denom != 0)
178 7656 : d = ((t - x[i]) * N[i]) / denom;
179 : }
180 : // else
181 : // d = 0.0 ;
182 :
183 27348 : if (N[i + 1] != 0) /* if the lower order basis function is zero skip
184 : the calculation */
185 : {
186 7656 : double denom = x[i + k] - x[i + 1];
187 7656 : if (denom != 0)
188 7656 : e = ((x[i + k] - t) * N[i + 1]) / denom;
189 : }
190 : // else
191 : // e = 0.0;
192 :
193 27348 : N[i] = d + e;
194 : }
195 : }
196 :
197 1322 : if (t == (double)x[nplusc])
198 : { /* pick up last point */
199 30 : N[npts] = 1;
200 : }
201 1322 : }
202 :
203 : /************************************************************************/
204 : /* rbasis() */
205 : /* */
206 : /* Subroutine to generate rational B-spline basis functions. */
207 : /* See the comment for basis() above. */
208 : /************************************************************************/
209 :
210 1320 : static void rbasis(int c, double t, int npts, double x[], double h[],
211 : double r[])
212 :
213 : {
214 1320 : const int nplusc = npts + c;
215 :
216 2640 : std::vector<double> temp;
217 1320 : temp.resize(nplusc + 1);
218 :
219 1320 : basis(c, t, npts, x, &temp[0]);
220 :
221 : /* calculate sum for denominator of rational basis functions */
222 :
223 1320 : double sum = 0.0;
224 9104 : for (int i = 1; i <= npts; i++)
225 : {
226 7784 : sum = sum + temp[i] * h[i];
227 : }
228 :
229 : /* form rational basis functions and put in r vector */
230 :
231 9104 : for (int i = 1; i <= npts; i++)
232 : {
233 7784 : if (sum != 0)
234 : {
235 7784 : r[i] = (temp[i] * h[i]) / (sum);
236 : }
237 : else
238 : {
239 0 : r[i] = 0;
240 : }
241 : }
242 1320 : }
243 :
244 : /************************************************************************/
245 : /* rbspline2() */
246 : /************************************************************************/
247 :
248 : /* Subroutine to generate a rational B-spline curve.
249 :
250 : C code for An Introduction to NURBS
251 : by David F. Rogers. Copyright (C) 2000 David F. Rogers,
252 : All rights reserved.
253 :
254 : Name: rbspline.c
255 : Language: C
256 : Subroutines called: knot.c, rbasis.c, fmtmul.c
257 : Book reference: Chapter 4, Alg. p. 297
258 :
259 : b[] = array containing the defining polygon vertices
260 : b[1] contains the x-component of the vertex
261 : b[2] contains the y-component of the vertex
262 : b[3] contains the z-component of the vertex
263 : h[] = array containing the homogeneous weighting factors
264 : k = order of the B-spline basis function
265 : nbasis = array containing the basis functions for a single value of t
266 : nplusc = number of knot values
267 : npts = number of defining polygon vertices
268 : p[,] = array containing the curve points
269 : p[1] contains the x-component of the point
270 : p[2] contains the y-component of the point
271 : p[3] contains the z-component of the point
272 : p1 = number of points to be calculated on the curve
273 : t = parameter value 0 <= t <= npts - k + 1
274 : bCalculateKnots = when set to true, x will be filled with the knot()
275 : routine, otherwise its content will be used. knots[] = array containing
276 : the knot vector (must be npts + k + 1 large)
277 : */
278 :
279 30 : void rbspline2(int npts, int k, int p1, double b[], double h[],
280 : bool bCalculateKnots, double knots[], double p[])
281 :
282 : {
283 30 : const int nplusc = npts + k;
284 :
285 60 : std::vector<double> nbasis;
286 30 : nbasis.resize(npts + 1);
287 :
288 : /* generate the uniform open knot vector */
289 :
290 30 : if (bCalculateKnots == true)
291 1 : knot(npts, k, knots);
292 :
293 30 : int icount = 0;
294 :
295 : /* calculate the points on the rational B-spline curve */
296 :
297 30 : double t = knots[1];
298 30 : const double step = (knots[nplusc] - knots[1]) / ((double)(p1 - 1));
299 :
300 30 : const double eps = 5e-6 * (knots[nplusc] - knots[1]);
301 :
302 1350 : for (int i1 = 1; i1 <= p1; i1++)
303 : {
304 : /* avoid undershooting the final knot */
305 1320 : if ((double)knots[nplusc] - t < eps)
306 : {
307 30 : t = (double)knots[nplusc];
308 : }
309 :
310 : /* generate the basis function for this value of t */
311 1320 : rbasis(k, t, npts, knots, h, &(nbasis[0]));
312 5280 : for (int j = 1; j <= 3; j++)
313 : { /* generate a point on the curve */
314 3960 : int jcount = j;
315 3960 : p[icount + j] = 0.;
316 :
317 27312 : for (int i = 1; i <= npts; i++)
318 : { /* Do local matrix multiplication */
319 23352 : const double temp = nbasis[i] * b[jcount];
320 23352 : p[icount + j] = p[icount + j] + temp;
321 23352 : jcount = jcount + 3;
322 : }
323 : }
324 1320 : icount = icount + 3;
325 1320 : t = t + step;
326 : }
327 30 : }
328 :
329 : /************************************************************************/
330 : /* rbspline() */
331 : /************************************************************************/
332 :
333 : /* Subroutine to generate a rational B-spline curve using an uniform open knot
334 : vector
335 :
336 : C code for An Introduction to NURBS
337 : by David F. Rogers. Copyright (C) 2000 David F. Rogers,
338 : All rights reserved.
339 :
340 : Name: rbspline.c
341 : Language: C
342 : Subroutines called: knot.c, rbasis.c, fmtmul.c
343 : Book reference: Chapter 4, Alg. p. 297
344 :
345 : b[] = array containing the defining polygon vertices
346 : b[1] contains the x-component of the vertex
347 : b[2] contains the y-component of the vertex
348 : b[3] contains the z-component of the vertex
349 : h[] = array containing the homogeneous weighting factors
350 : k = order of the B-spline basis function
351 : nbasis = array containing the basis functions for a single value of t
352 : nplusc = number of knot values
353 : npts = number of defining polygon vertices
354 : p[,] = array containing the curve points
355 : p[1] contains the x-component of the point
356 : p[2] contains the y-component of the point
357 : p[3] contains the z-component of the point
358 : p1 = number of points to be calculated on the curve
359 : t = parameter value 0 <= t <= npts - k + 1
360 : x[] = array containing the knot vector
361 : */
362 :
363 0 : void rbspline(int npts, int k, int p1, double b[], double h[], double p[])
364 :
365 : {
366 0 : std::vector<double> x(npts + k + 1, 0.0);
367 0 : rbspline2(npts, k, p1, b, h, true, &x[0], p);
368 0 : }
369 :
370 : /************************************************************************/
371 : /* rbsplinu() */
372 : /************************************************************************/
373 :
374 : /* Subroutine to generate a rational B-spline curve using an uniform periodic
375 : knot vector
376 :
377 : C code for An Introduction to NURBS
378 : by David F. Rogers. Copyright (C) 2000 David F. Rogers,
379 : All rights reserved.
380 :
381 : Name: rbsplinu.c
382 : Language: C
383 : Subroutines called: knotu.c, rbasis.c, fmtmul.c
384 : Book reference: Chapter 4, Alg. p. 298
385 :
386 : b[] = array containing the defining polygon vertices
387 : b[1] contains the x-component of the vertex
388 : b[2] contains the y-component of the vertex
389 : b[3] contains the z-component of the vertex
390 : h[] = array containing the homogeneous weighting factors
391 : k = order of the B-spline basis function
392 : nbasis = array containing the basis functions for a single value of t
393 : nplusc = number of knot values
394 : npts = number of defining polygon vertices
395 : p[,] = array containing the curve points
396 : p[1] contains the x-component of the point
397 : p[2] contains the y-component of the point
398 : p[3] contains the z-component of the point
399 : p1 = number of points to be calculated on the curve
400 : t = parameter value 0 <= t <= npts - k + 1
401 : x[] = array containing the knot vector
402 : */
403 :
404 0 : void rbsplinu(int npts, int k, int p1, double b[], double h[], double p[])
405 :
406 : {
407 : int i, j, icount, jcount;
408 : int i1;
409 : int nplusc;
410 :
411 : double step;
412 : double t;
413 : double temp;
414 0 : std::vector<double> nbasis;
415 0 : std::vector<double> x;
416 :
417 0 : nplusc = npts + k;
418 :
419 0 : x.resize(nplusc + 1);
420 0 : nbasis.resize(npts + 1);
421 :
422 : /* zero and redimension the knot vector and the basis array */
423 :
424 0 : for (i = 0; i <= npts; i++)
425 : {
426 0 : nbasis[i] = 0.;
427 : }
428 :
429 0 : for (i = 0; i <= nplusc; i++)
430 : {
431 0 : x[i] = 0;
432 : }
433 :
434 : /* generate the uniform periodic knot vector */
435 :
436 0 : knotu(npts, k, &(x[0]));
437 :
438 0 : icount = 0;
439 :
440 : /* calculate the points on the rational B-spline curve */
441 :
442 0 : t = k - 1;
443 0 : step = ((double)((npts) - (k - 1))) / ((double)(p1 - 1));
444 :
445 0 : for (i1 = 1; i1 <= p1; i1++)
446 : {
447 :
448 0 : if (x[nplusc] - t < 5e-6)
449 : {
450 0 : t = x[nplusc];
451 : }
452 :
453 0 : rbasis(
454 0 : k, t, npts, &(x[0]), h,
455 0 : &(nbasis[0])); /* generate the basis function for this value of t */
456 :
457 0 : for (j = 1; j <= 3; j++)
458 : { /* generate a point on the curve */
459 0 : jcount = j;
460 0 : p[icount + j] = 0.;
461 :
462 0 : for (i = 1; i <= npts; i++)
463 : { /* Do local matrix multiplication */
464 0 : temp = nbasis[i] * b[jcount];
465 0 : p[icount + j] = p[icount + j] + temp;
466 0 : jcount = jcount + 3;
467 : }
468 : }
469 0 : icount = icount + 3;
470 0 : t = t + step;
471 : }
472 0 : }
|
