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