1 /*
2 * Copyright (C) 1998, 2000-2007, 2010, 2011, 2012, 2013 SINTEF ICT,
3 * Applied Mathematics, Norway.
4 *
5 * Contact information: E-mail: tor.dokken@sintef.no
6 * SINTEF ICT, Department of Applied Mathematics,
7 * P.O. Box 124 Blindern,
8 * 0314 Oslo, Norway.
9 *
10 * This file is part of SISL.
11 *
12 * SISL is free software: you can redistribute it and/or modify
13 * it under the terms of the GNU Affero General Public License as
14 * published by the Free Software Foundation, either version 3 of the
15 * License, or (at your option) any later version.
16 *
17 * SISL is distributed in the hope that it will be useful,
18 * but WITHOUT ANY WARRANTY; without even the implied warranty of
19 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
20 * GNU Affero General Public License for more details.
21 *
22 * You should have received a copy of the GNU Affero General Public
23 * License along with SISL. If not, see
24 * <http://www.gnu.org/licenses/>.
25 *
26 * In accordance with Section 7(b) of the GNU Affero General Public
27 * License, a covered work must retain the producer line in every data
28 * file that is created or manipulated using SISL.
29 *
30 * Other Usage
31 * You can be released from the requirements of the license by purchasing
32 * a commercial license. Buying such a license is mandatory as soon as you
33 * develop commercial activities involving the SISL library without
34 * disclosing the source code of your own applications.
35 *
36 * This file may be used in accordance with the terms contained in a
37 * written agreement between you and SINTEF ICT.
38 */
39
40 #include "sisl-copyright.h"
41
42 /*
43 *
44 * $Id: s1938.c,v 1.2 2001-03-19 15:58:56 afr Exp $
45 *
46 */
47
48
49 #define S1938
50
51 #include "sislP.h"
52 #define MAX_SIZE 50
53
54
55 #if defined(SISLNEEDPROTOTYPES)
56 void
s1938(SISLSurf * srf,double etr1[],int inr1,double etr2[],int inr2,double ** surfr,int * jstat)57 s1938 (SISLSurf * srf, double etr1[], int inr1, double etr2[], int inr2,
58 double **surfr, int *jstat)
59 #else
60 void
61 s1938 (srf, etr1, inr1, etr2, inr2, surfr, jstat)
62 SISLSurf *srf;
63 double etr1[];
64 int inr1;
65 double etr2[];
66 int inr2;
67 double **surfr;
68 int *jstat;
69
70 #endif
71 /*
72 *********************************************************************
73 *
74 *********************************************************************
75 *
76 * PURPOSE: To express a B-spline surface curve using a refined basis.
77 *
78 *
79 * INPUT: srf - The original surface.
80 * etr1 - The refined knot vector in the first parameter
81 * direction.
82 * inr1 - The number of vertices of the subdivided surface
83 * in first parameter direction.
84 * etr2 - The refined knot vector in the second parameter
85 * direction.
86 * inr2 - The number of vertices of the subdivided surface
87 * in the second parameter direction.
88 *
89 * OUTPUT: surfr - Array containing the vertices of the refined
90 * surface.
91 * jstat - Status variable.
92 * < 0: Error.
93 * = 0: Ok.
94 * > 0: Warning.
95 *
96 * METHOD: The vertices are calculated using the "Oslo"-algorithm
97 * developped by Cohen, Lyche and Riesenfeld.
98 *
99 *
100 * REFERENCES: Cohen, Lyche, Riesenfeld: Discrete B-splines and subdivision
101 * techniques in computer aided geometric design, computer
102 * graphics and image processing, vol 14, no.2 (1980)
103 *
104 * CALLS: s1937, s6err.
105 *
106 *
107 * WRITTEN BY: Christophe R. Birkeland, SI, 1991-08
108 * REVISED BY: Johannes Kaasa, SI, May 1992 (Introduced NURBS)
109 *
110 *********************************************************************
111 */
112 {
113 int kpos = 0;
114 int ki, kj, kl, kr, low; /* Loop control parameters */
115 int start;
116 int rem; /* Used to store array index */
117 int nu; /* Pointer into knot vector:
118 knt[nu-1]<=etd[kj]<knt[nu] */
119 int ins1; /* Number of vertices of curve in first
120 parameter direction */
121 int iordr1; /* Order of curve in first
122 parameter direction */
123 int ins2; /* Number of vertices of curve in second
124 parameter direction */
125 int iordr2; /* Order of curve in second
126 parameter direction */
127 double *knt1 = SISL_NULL; /* Original knot-vector of surface. */
128 double *knt2 = SISL_NULL; /* Original knot vector of surface. */
129 double *coef = SISL_NULL; /* Pointer to array of coefficients of
130 the surface */
131 int idim; /* Dimension of space where the
132 curve lies */
133 double sum; /* Used to store vertices of
134 new curve */
135 double sarray[MAX_SIZE];
136 int alloc_needed=FALSE;
137 double *alfa = SISL_NULL; /* Array needed in subroutine
138 s1937 (Oslo-algorithm) */
139 double *ktsurf = SISL_NULL; /* Array for internal use only */
140
141 *jstat = 0;
142
143
144 /* Initialization. */
145
146 knt1 = srf->et1;
147 knt2 = srf->et2;
148 ins1 = srf->in1;
149 ins2 = srf->in2;
150 iordr1 = srf->ik1;
151 iordr2 = srf->ik2;
152 if (srf->ikind == 2 || srf->ikind == 4)
153 {
154 idim = srf->idim + 1;
155 coef = srf->rcoef;
156 }
157 else
158 {
159 idim = srf->idim;
160 coef = srf->ecoef;
161 }
162
163 /* Test if legal input. */
164
165 if (iordr1 < 1 || iordr2 < 1)
166 goto err115;
167 if (ins1 < iordr1 || ins2 < iordr2)
168 goto err116;
169 if (idim < 1)
170 goto err102;
171
172
173 /* Allocate array for internal use only. */
174
175 if (MAX(iordr1,iordr2) > MAX_SIZE)
176 {
177 if ((alfa = newarray(MAX(iordr1,iordr2),DOUBLE)) == SISL_NULL)
178 goto err101;
179 alloc_needed = TRUE;
180 }
181 else
182 alfa = sarray;
183
184 ktsurf = newarray (inr1 * inr2 * idim, DOUBLE);
185 if (ktsurf == SISL_NULL)
186 goto err101;
187
188
189 /* Allocate array surfr for output. */
190
191 *surfr = newarray (inr1 * inr2 * idim, DOUBLE);
192 if (*surfr == SISL_NULL)
193 goto err101;
194
195
196 /* Find if etr1 is a refinement of the original
197 knot vector knt1 (srf->et1). */
198
199 kj = iordr1 - 1;
200
201 for (ki = 0; kj < ins1; ki++)
202 {
203 if (ki >= inr1)
204 goto err116;
205 if (knt1[kj] > etr1[ki])
206 continue;
207 if (knt1[kj] < etr1[ki])
208 goto err117;
209 kj++;
210 }
211
212 /* etr1 is a refinement of original knot vector knt1
213 * Produce surface refined in one direction. */
214
215 nu = 1;
216 for (kj = 0; kj < inr1; kj++)
217 {
218 /* We want to find knt1[nu-1] <= etr1[kj] < knt1[nu]
219 The copying of knots guarantees the nu-value to be found.
220 Since kj is increasing, the nu-values will be increasing
221 due to copying of knots. */
222
223 for (; (((knt1[nu - 1] > etr1[kj]) || (etr1[kj] >= knt1[nu]))
224 && (nu != ins1)); nu++) ;
225
226
227 /* Now we have knt1[nu-1] <= etr1[kj] < knt1[nu],
228 so the discrete B-splines can be calculated. */
229
230 s1937 (knt1, iordr1, kj + 1, nu, alfa, etr1);
231
232
233 /* Compute the temporary surface. */
234
235 low = nu - iordr1 + 1;
236 for (ki = 0; ki < ins2; ki++)
237 {
238 rem = idim * kj + idim * inr1 * ki;
239
240 for (kl = 0; kl < idim; kl++)
241 {
242 sum = (double) 0.0;
243 start = nu - iordr1 + 1;
244 if (start < 1)
245 start = 1;
246
247 for (kr = start; kr <= nu; kr++)
248 sum += alfa[kr - low] * coef[ki * ins1 * idim + (kr - 1) * idim + kl];
249 ktsurf[rem + kl] = sum;
250 }
251 }
252 }
253
254 /* Find if etr2 is a refinement of the original
255 * knot vector knt2 (srf->et2). */
256
257 kj = iordr2 - 1;
258
259 for (ki = 0; kj < ins2; ki++)
260 {
261 if (ki >= inr2)
262 goto err116;
263 if (knt2[kj] > etr2[ki])
264 continue;
265 if (knt2[kj] < etr2[ki])
266 goto err117;
267 kj++;
268 }
269
270 /* etr2 is a refinement of original knot vector knt2
271 * Produce surface refined in one direction. */
272
273 nu = 1;
274 for (ki = 0; ki < inr2; ki++)
275 {
276 /* We want to find knt2[nu-1] <= etr2[ki] < knt2[nu]
277 The copying of knots guarantees the nu-value to be found.
278 Since kj is increasing, the nu-values will be increasing
279 due to copying of knots. */
280
281 for (; (((knt2[nu - 1] > etr2[ki]) || (etr2[ki] >= knt2[nu]))
282 && (nu != ins2)); nu++) ;
283
284
285 /* Now we have knt2[nu-1] <= etr2[kj] < knt2[nu],
286 so the discrete B-splines can be calculated */
287
288 s1937 (knt2, iordr2, ki + 1, nu, alfa, etr2);
289
290
291 /* Compute the temporary surface. */
292
293 low = nu - iordr2 + 1;
294 for (kj = 0; kj < inr1; kj++)
295 for (kl = 0; kl < idim; kl++)
296 {
297 sum = (double) 0.0;
298 start = nu - iordr2 + 1;
299 if (start < 1)
300 start = 1;
301 for (kr = start; kr <= nu; kr++)
302 sum += alfa[kr - low] * ktsurf[kj * idim + (kr - 1) * idim * inr1 + kl];
303 (*surfr)[ki * idim * inr1 + kj * idim + kl] = sum;
304 }
305 }
306
307
308 /* OK. */
309
310 goto out;
311
312
313 /* Memory error */
314
315 err101:
316 *jstat = -101;
317 s6err ("s1938", *jstat, kpos);
318 goto out;
319
320 /* Error in B-spline surface description:
321
322 Dimension less than 1. */
323
324
325 err102:
326 *jstat = -102;
327 s6err ("s1938", *jstat, kpos);
328 goto out;
329
330 /* Order less than 1. */
331
332 err115:
333 *jstat = -115;
334 s6err ("s1938", *jstat, kpos);
335 goto out;
336
337 /* No. of vertices less than order. */
338
339 err116:
340 *jstat = -116;
341 s6err ("s1938", *jstat, kpos);
342 goto out;
343
344 /* Error in knot vector */
345
346 err117:
347 *jstat = -117;
348 s6err ("s1938", *jstat, kpos);
349 goto out;
350
351 out:
352 if (alloc_needed)
353 freearray (alfa);
354 if (ktsurf != SISL_NULL)
355 freearray (ktsurf);
356
357 return;
358 }
359 #undef MAX_SIZE
360