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, std::vector<double> &x)
76 :
77 : {
78 : int nplusc, /* nplus2, */ i;
79 :
80 0 : nplusc = n + c;
81 : /* nplus2 = n + 2; */
82 :
83 : #if defined(__GNUC__)
84 : #pragma GCC diagnostic push
85 : #pragma GCC diagnostic ignored "-Wnull-dereference"
86 : #endif
87 :
88 0 : x[1] = 0;
89 0 : for (i = 2; i <= nplusc; i++)
90 : {
91 0 : x[i] = i - 1;
92 : }
93 :
94 : #if defined(__GNUC__)
95 : #pragma GCC diagnostic pop
96 : #endif
97 0 : }
98 :
99 : /************************************************************************/
100 : /* knot() */
101 : /************************************************************************/
102 :
103 : /*
104 : Subroutine to generate a B-spline open knot vector with multiplicity
105 : equal to the order at the ends.
106 :
107 : c = order of the basis function
108 : n = the number of defining polygon vertices
109 : nplus2 = index of x() for the first occurrence of the maximum knot
110 : vector value nplusc = maximum value of the knot vector -- $n + c$ x()
111 : = array containing the knot vector
112 : */
113 :
114 1 : static void knot(int n, int c, double x[])
115 :
116 : {
117 1 : const int nplusc = n + c;
118 1 : const int nplus2 = n + 2;
119 :
120 1 : x[1] = 0.0;
121 11 : for (int i = 2; i <= nplusc; i++)
122 : {
123 10 : if ((i > c) && (i < nplus2))
124 4 : x[i] = x[i - 1] + 1.0;
125 : else
126 6 : x[i] = x[i - 1];
127 : }
128 1 : }
129 :
130 : /************************************************************************/
131 : /* basis() */
132 : /************************************************************************/
133 :
134 : /* Subroutine to generate B-spline basis functions--open knot vector
135 :
136 : C code for An Introduction to NURBS
137 : by David F. Rogers. Copyright (C) 2000 David F. Rogers,
138 : All rights reserved.
139 :
140 : Name: rbasis (split into rbasis and basis by AT)
141 : Language: C
142 : Subroutines called: none
143 : Book reference: Chapter 4, Sec. 4. , p 296
144 :
145 : c = order of the B-spline basis function
146 : d = first term of the basis function recursion relation
147 : e = second term of the basis function recursion relation
148 : h[] = array containing the homogeneous weights
149 : npts = number of defining polygon vertices
150 : nplusc = constant -- npts + c -- maximum number of knot values
151 : r[] = array containing the rationalbasis functions
152 : r[1] contains the basis function associated with B1 etc.
153 : t = parameter value
154 : N[] = array output from bbasis containing the values of the
155 : basis functions; should have capacity npts + c + 1
156 : x[] = knot vector
157 : */
158 :
159 1322 : void basis(int c, double t, int npts, double x[], double N[])
160 :
161 : {
162 1322 : const int nplusc = npts + c;
163 :
164 : /* calculate the first order nonrational basis functions n[i] */
165 :
166 13082 : for (int i = 1; i <= nplusc - 1; i++)
167 : {
168 11760 : if ((t >= x[i]) && (t < x[i + 1]))
169 1292 : N[i] = 1.0;
170 : else
171 10468 : N[i] = 0.0;
172 : }
173 :
174 : /* calculate the higher order nonrational basis functions */
175 :
176 5288 : for (int k = 2; k <= c; k++)
177 : {
178 31314 : for (int i = 1; i <= nplusc - k; i++)
179 : {
180 27348 : double d = 0.0;
181 27348 : double e = 0.0;
182 27348 : if (N[i] != 0) /* if the lower order basis function is zero skip the
183 : calculation */
184 : {
185 7656 : double denom = x[i + k - 1] - x[i];
186 7656 : if (denom != 0)
187 7656 : d = ((t - x[i]) * N[i]) / denom;
188 : }
189 : // else
190 : // d = 0.0 ;
191 :
192 27348 : if (N[i + 1] != 0) /* if the lower order basis function is zero skip
193 : the calculation */
194 : {
195 7656 : double denom = x[i + k] - x[i + 1];
196 7656 : if (denom != 0)
197 7656 : e = ((x[i + k] - t) * N[i + 1]) / denom;
198 : }
199 : // else
200 : // e = 0.0;
201 :
202 27348 : N[i] = d + e;
203 : }
204 : }
205 :
206 1322 : if (t == (double)x[nplusc])
207 : { /* pick up last point */
208 30 : N[npts] = 1;
209 : }
210 1322 : }
211 :
212 : /************************************************************************/
213 : /* rbasis() */
214 : /* */
215 : /* Subroutine to generate rational B-spline basis functions. */
216 : /* See the comment for basis() above. */
217 : /************************************************************************/
218 :
219 1320 : static void rbasis(int c, double t, int npts, double x[], double h[],
220 : double r[])
221 :
222 : {
223 1320 : const int nplusc = npts + c;
224 :
225 2640 : std::vector<double> temp;
226 1320 : temp.resize(nplusc + 1);
227 :
228 1320 : basis(c, t, npts, x, &temp[0]);
229 :
230 : /* calculate sum for denominator of rational basis functions */
231 :
232 1320 : double sum = 0.0;
233 9104 : for (int i = 1; i <= npts; i++)
234 : {
235 7784 : sum = sum + temp[i] * h[i];
236 : }
237 :
238 : /* form rational basis functions and put in r vector */
239 :
240 9104 : for (int i = 1; i <= npts; i++)
241 : {
242 7784 : if (sum != 0)
243 : {
244 7784 : r[i] = (temp[i] * h[i]) / (sum);
245 : }
246 : else
247 : {
248 0 : r[i] = 0;
249 : }
250 : }
251 1320 : }
252 :
253 : /************************************************************************/
254 : /* rbspline2() */
255 : /************************************************************************/
256 :
257 : /* Subroutine to generate a rational B-spline curve.
258 :
259 : C code for An Introduction to NURBS
260 : by David F. Rogers. Copyright (C) 2000 David F. Rogers,
261 : All rights reserved.
262 :
263 : Name: rbspline.c
264 : Language: C
265 : Subroutines called: knot.c, rbasis.c, fmtmul.c
266 : Book reference: Chapter 4, Alg. p. 297
267 :
268 : b[] = array containing the defining polygon vertices
269 : b[1] contains the x-component of the vertex
270 : b[2] contains the y-component of the vertex
271 : b[3] contains the z-component of the vertex
272 : h[] = array containing the homogeneous weighting factors
273 : k = order of the B-spline basis function
274 : nbasis = array containing the basis functions for a single value of t
275 : nplusc = number of knot values
276 : npts = number of defining polygon vertices
277 : p[,] = array containing the curve points
278 : p[1] contains the x-component of the point
279 : p[2] contains the y-component of the point
280 : p[3] contains the z-component of the point
281 : p1 = number of points to be calculated on the curve
282 : t = parameter value 0 <= t <= npts - k + 1
283 : bCalculateKnots = when set to true, x will be filled with the knot()
284 : routine, otherwise its content will be used. knots[] = array containing
285 : the knot vector (must be npts + k + 1 large)
286 : */
287 :
288 30 : void rbspline2(int npts, int k, int p1, double b[], double h[],
289 : bool bCalculateKnots, double knots[], double p[])
290 :
291 : {
292 30 : const int nplusc = npts + k;
293 :
294 60 : std::vector<double> nbasis;
295 30 : nbasis.resize(npts + 1);
296 :
297 : /* generate the uniform open knot vector */
298 :
299 30 : if (bCalculateKnots == true)
300 1 : knot(npts, k, knots);
301 :
302 30 : int icount = 0;
303 :
304 : /* calculate the points on the rational B-spline curve */
305 :
306 30 : double t = knots[1];
307 30 : const double step = (knots[nplusc] - knots[1]) / ((double)(p1 - 1));
308 :
309 30 : const double eps = 5e-6 * (knots[nplusc] - knots[1]);
310 :
311 1350 : for (int i1 = 1; i1 <= p1; i1++)
312 : {
313 : /* avoid undershooting the final knot */
314 1320 : if ((double)knots[nplusc] - t < eps)
315 : {
316 30 : t = (double)knots[nplusc];
317 : }
318 :
319 : /* generate the basis function for this value of t */
320 1320 : rbasis(k, t, npts, knots, h, &(nbasis[0]));
321 5280 : for (int j = 1; j <= 3; j++)
322 : { /* generate a point on the curve */
323 3960 : int jcount = j;
324 3960 : p[icount + j] = 0.;
325 :
326 27312 : for (int i = 1; i <= npts; i++)
327 : { /* Do local matrix multiplication */
328 23352 : const double temp = nbasis[i] * b[jcount];
329 23352 : p[icount + j] = p[icount + j] + temp;
330 23352 : jcount = jcount + 3;
331 : }
332 : }
333 1320 : icount = icount + 3;
334 1320 : t = t + step;
335 : }
336 30 : }
337 :
338 : /************************************************************************/
339 : /* rbspline() */
340 : /************************************************************************/
341 :
342 : /* Subroutine to generate a rational B-spline curve using an uniform open knot
343 : vector
344 :
345 : C code for An Introduction to NURBS
346 : by David F. Rogers. Copyright (C) 2000 David F. Rogers,
347 : All rights reserved.
348 :
349 : Name: rbspline.c
350 : Language: C
351 : Subroutines called: knot.c, rbasis.c, fmtmul.c
352 : Book reference: Chapter 4, Alg. p. 297
353 :
354 : b[] = array containing the defining polygon vertices
355 : b[1] contains the x-component of the vertex
356 : b[2] contains the y-component of the vertex
357 : b[3] contains the z-component of the vertex
358 : h[] = array containing the homogeneous weighting factors
359 : k = order of the B-spline basis function
360 : nbasis = array containing the basis functions for a single value of t
361 : nplusc = number of knot values
362 : npts = number of defining polygon vertices
363 : p[,] = array containing the curve points
364 : p[1] contains the x-component of the point
365 : p[2] contains the y-component of the point
366 : p[3] contains the z-component of the point
367 : p1 = number of points to be calculated on the curve
368 : t = parameter value 0 <= t <= npts - k + 1
369 : x[] = array containing the knot vector
370 : */
371 :
372 0 : void rbspline(int npts, int k, int p1, double b[], double h[], double p[])
373 :
374 : {
375 0 : std::vector<double> x(npts + k + 1, 0.0);
376 0 : rbspline2(npts, k, p1, b, h, true, &x[0], p);
377 0 : }
378 :
379 : /************************************************************************/
380 : /* rbsplinu() */
381 : /************************************************************************/
382 :
383 : /* Subroutine to generate a rational B-spline curve using an uniform periodic
384 : knot vector
385 :
386 : C code for An Introduction to NURBS
387 : by David F. Rogers. Copyright (C) 2000 David F. Rogers,
388 : All rights reserved.
389 :
390 : Name: rbsplinu.c
391 : Language: C
392 : Subroutines called: knotu.c, rbasis.c, fmtmul.c
393 : Book reference: Chapter 4, Alg. p. 298
394 :
395 : b[] = array containing the defining polygon vertices
396 : b[1] contains the x-component of the vertex
397 : b[2] contains the y-component of the vertex
398 : b[3] contains the z-component of the vertex
399 : h[] = array containing the homogeneous weighting factors
400 : k = order of the B-spline basis function
401 : nbasis = array containing the basis functions for a single value of t
402 : nplusc = number of knot values
403 : npts = number of defining polygon vertices
404 : p[,] = array containing the curve points
405 : p[1] contains the x-component of the point
406 : p[2] contains the y-component of the point
407 : p[3] contains the z-component of the point
408 : p1 = number of points to be calculated on the curve
409 : t = parameter value 0 <= t <= npts - k + 1
410 : x[] = array containing the knot vector
411 : */
412 :
413 0 : void rbsplinu(int npts, int k, int p1, double b[], double h[], double p[])
414 :
415 : {
416 : int i, j, icount, jcount;
417 : int i1;
418 : int nplusc;
419 :
420 : double step;
421 : double t;
422 : double temp;
423 0 : std::vector<double> nbasis;
424 0 : std::vector<double> x;
425 :
426 0 : nplusc = npts + k;
427 :
428 0 : x.resize(nplusc + 1);
429 0 : nbasis.resize(npts + 1);
430 :
431 : /* zero and redimension the knot vector and the basis array */
432 :
433 0 : for (i = 0; i <= npts; i++)
434 : {
435 0 : nbasis[i] = 0.;
436 : }
437 :
438 0 : for (i = 0; i <= nplusc; i++)
439 : {
440 0 : x[i] = 0;
441 : }
442 :
443 : /* generate the uniform periodic knot vector */
444 :
445 0 : knotu(npts, k, x);
446 :
447 0 : icount = 0;
448 :
449 : /* calculate the points on the rational B-spline curve */
450 :
451 0 : t = k - 1;
452 0 : step = ((double)((npts) - (k - 1))) / ((double)(p1 - 1));
453 :
454 0 : for (i1 = 1; i1 <= p1; i1++)
455 : {
456 :
457 0 : if (x[nplusc] - t < 5e-6)
458 : {
459 0 : t = x[nplusc];
460 : }
461 :
462 0 : rbasis(
463 0 : k, t, npts, &(x[0]), h,
464 0 : &(nbasis[0])); /* generate the basis function for this value of t */
465 :
466 0 : for (j = 1; j <= 3; j++)
467 : { /* generate a point on the curve */
468 0 : jcount = j;
469 0 : p[icount + j] = 0.;
470 :
471 0 : for (i = 1; i <= npts; i++)
472 : { /* Do local matrix multiplication */
473 0 : temp = nbasis[i] * b[jcount];
474 0 : p[icount + j] = p[icount + j] + temp;
475 0 : jcount = jcount + 3;
476 : }
477 : }
478 0 : icount = icount + 3;
479 0 : t = t + step;
480 : }
481 0 : }
|
