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: s1933.c,v 1.4 2001-03-19 15:58:56 afr Exp $
45 *
46 */
47
48
49 #define S1933
50
51 #include "sislP.h"
52
53
54 #if defined(SISLNEEDPROTOTYPES)
55 void
s1933(int inbcrv,SISLCurve * crvarr[],double start,double stop,double ** it,int * in,int * iordr,int * jstat)56 s1933 (int inbcrv, SISLCurve * crvarr[], double start, double stop,
57 double **it, int *in, int *iordr, int *jstat)
58 #else
59 void
60 s1933 (inbcrv, crvarr, start, stop, it, in, iordr, jstat)
61 int inbcrv;
62 SISLCurve *crvarr[];
63 double start;
64 double stop;
65 double **it;
66 int *in;
67 int *iordr;
68 int *jstat;
69
70 #endif
71 /*
72 *********************************************************************
73 *
74 *********************************************************************
75 *
76 * PURPOSE : To map a number of knot vectors of different lengths
77 * into the same interval, and to produce a union of the
78 * knot vectors taking care to use the greatest of the
79 * orders of the input knot vectors and to reflect the
80 * continuity requirements reflected in all input knot
81 * vectors. To avoid too many new knots, knots lying
82 * closer than 0.0001 times the parameter interval
83 * (relative distance), might be moved.
84 *
85 *
86 * INPUT : inbcrv - Number of curves to be interpolated by a
87 * spline lofted surface.
88 * crvarr - Array (length inbcrv) of pointers to the curves
89 * in the curve-set.
90 * start - Start parameter value of B-spline basis to
91 * be made.
92 * stop - End parameter value of B-spline basis to
93 * be made.
94 *
95 *
96 * OUTPUT : it - The new knot vector (union of the old ones).
97 * in - The number of B-spline basis functions in basis
98 * produced.
99 * iordr - The order of the B-spline basis produced.
100 * jstat - Output status:
101 * < 0: Error.
102 * = 0: Ok.
103 * > 0: Warning.
104 *
105 * METHOD : Phase 1:
106 * Find highest order of all curves.
107 * Map all knot-vectors into [astart,astop].
108 * Lift all knot-vectors to the highest order and make union.
109 * Phase 2:
110 * Run through the union vector and detect knots lying
111 * closer to each other than 10E-4 times the parameter
112 * interval.
113 * Find values to be used for the close knot-values.
114 * Move knots in original knot-vectors to match these values.
115 * Phase 3:
116 * Lift all knot-vectors (now adjusted) to the highest
117 * order and make union.
118 *
119 * REFERENCES : Fortran version:
120 * T.Dokken, SI, 1981-12
121 *
122 * CALLS : s1934,s1754,s1935,s6err.
123 *
124 *
125 * WRITTEN BY : Christophe R. Birkeland, SI, 1991-07
126 * REVISED BY : Paal Fugelli, SINTEF, Oslo 18/07-1994. Removed memory
127 * leaks.
128 *
129 *********************************************************************
130 */
131 {
132 int ki, kj, kl, kr, kp; /* Loop control variables */
133 int ktell;
134 int kfirst; /* Pointer to first knot value in
135 * interval [tinmin,tinmax] */
136 int klast; /* Pointer to last knot value in
137 * interval [tinmin,tinmax] */
138 int kn; /* Number of vertices in curve */
139 int kbgn; /* Pointer to first knot value
140 * satisfying kt[kbgn] > kt[*iordr-1] */
141 int kend; /* Pointer to last knot value
142 * satisfying kt[kend] < kt[kn] */
143 int kmax; /* Number of distinct knot values to be
144 * kept in the interval [tinmin,tinmax] */
145 int kdum1; /* Dummy variables for use in algorithm */
146 int kdum2;
147 int kdiff;
148 int kdf; /* Equals (kdiff div 2) */
149 int knumb; /* Number of vertices in a curve */
150 int kpos = 0; /* Error position indicator */
151 int kstat;
152 double tepsco; /* Greatest distance for two knots to be
153 * regarded as the same knot */
154 double dum;
155 double tlast; /* Value of knot vector */
156 double tval;
157 double tincre;
158 double tincr2; /* Equals tincre / 2.0 */
159 double tmin, tmax;
160 double tinmin, tinmax;
161 double tbgn, tend; /* tbgn=kt[*iordr-1] & tend = kt[kn]
162 * Used to find kbgn and kend */
163 double *kt = SISL_NULL; /* New knot vector */
164 double *knot = SISL_NULL; /* Knot vector */
165 double *incknt = SISL_NULL; /* Used to store the union of two
166 * knot vectors */
167 SISLCurve *curve = SISL_NULL;
168
169 *jstat = 0;
170
171
172 /* Initailzation of variables */
173
174 *in = 0;
175 *iordr = 0;
176
177
178 /* Test if legal input */
179
180 if (inbcrv < 2)
181 goto err179;
182 for (ki = 0; ki < inbcrv; ki++)
183 if (crvarr[ki]->in <crvarr[ki]->ik || crvarr[ki]->ik < 1)
184 goto err112;
185
186
187 /* Find highest order used in description */
188
189 for (ki = 0; ki < inbcrv; ki++)
190 *iordr = MAX (*iordr, crvarr[ki]->ik);
191
192 kn = *iordr;
193
194
195 /* Normalize the knot vectors */
196
197 for (ki = 0; ki < inbcrv; ki++)
198 {
199 curve = crvarr[ki];
200 s1934 (curve->et, curve->in, curve->ik, start, stop, &kstat);
201 if (kstat < 0)
202 goto error;
203 }
204
205 /* Make start knot vector */
206
207 kt = newarray (*iordr * 2, double);
208 if (kt == SISL_NULL)
209 goto err101;
210 for (ki = 0; ki < *iordr; ki++)
211 {
212 kt[ki] = start;
213 kt[*iordr + ki] = stop;
214 }
215
216 /* Run through all knot vectors, lift the order, map into right
217 * interval and make the union with the candidates already existing */
218
219 for (ki = 0; ki < inbcrv; ki++)
220 {
221 /* Increase order of basis */
222
223 curve = crvarr[ki];
224 s1754 (curve->et, curve->in, curve->ik, *iordr,
225 &incknt, &knumb, &kstat);
226 if (kstat < 0)
227 goto error;
228
229
230 /* Make union of old union and new knot vector */
231
232 s1935 (kt, kn, incknt, knumb, &knot, &kn, *iordr, &kstat);
233 if (kstat < 0)
234 goto error;
235
236 if (incknt != SISL_NULL) freearray(incknt); /* PFU 18/07-94. */
237
238 if (kt != SISL_NULL)
239 freearray (kt);
240 kt = knot;
241 knot = SISL_NULL; /* PFU 18/07-94 */
242 }
243
244 /* The knot vector produced might contain knots originating
245 * from different curves that are located very close to each
246 * other. We will move internal knots when they are closer
247 * to each other than tepsco.
248 *
249 * Find first knot bigger than kt[*iordr-1] and last
250 * knot less than kt[kn]. */
251
252 kend = kn - 1;
253 tbgn = kt[*iordr - 1];
254 tend = kt[kn];
255
256 /* Find first knot kt[kbgn] > tbgn */
257
258 for (kbgn = *iordr; (tbgn >= kt[kbgn]) && (kbgn < kend); kbgn++) ;
259
260 /* Find last knot kt[kend] < tend */
261
262 for (; (tend <= kt[kend]) && (kbgn < kend); kend--) ;
263
264 /* Set greatest distance for two knots to be regarded
265 * as the same knot */
266
267 tepsco = (double) 0.00000000000001 * (stop - start);
268
269 /* Loop runing through all the knots of the union knot
270 * vector and moving knots when possible */
271
272 for (ki = kbgn; ki < kend; ki = kl + 1)
273 {
274 /* Find knots closer to knot ki than tepsco in
275 * positive direction */
276
277 kl = ki - 1;
278 do
279 {
280 kl++;
281 dum = MAX (fabs(kt[kl + 1]), fabs(kt[ki]))/stop;
282 if (dum == (double) 0.0)
283 dum = (double) 1.0;
284 }
285 while ((kl < kend) && ((fabs (kt[kl + 1] - kt[ki]) / dum) <= tepsco));
286
287
288 /* If only one value found, then shifting is not necessary */
289
290 if (ki == kl)
291 continue;
292
293
294 /* The knots close to kt[ki-1] are kt[ki]...kt[kl-1] */
295
296 tmin = kt[ki];
297 tmax = MIN (kt[kl] + tepsco*dum, kt[kend]);
298
299
300 /* If interval is degenerate we have finished the moving */
301
302 if (tmin >= tmax)
303 break;
304
305
306 /* For each curve, count the number of distinct knot
307 values in the interval [tmin,tmax].
308 Accumulate max and min knot values in the interval
309 where we know that tmin is a knot. */
310
311 tinmin = tmin;
312 tinmax = tmax;
313 kmax = 0;
314 for (kj = 0; kj < inbcrv; kj++)
315 {
316 curve = crvarr[kj];
317 ktell = 0;
318 tlast = (double) 0.0;
319
320 for (kr = *iordr - 1; kr <= curve->in; kr++)
321 if ((tmin <= curve->et[kr]) && (curve->et[kr] <= tmax))
322 {
323 if (ktell == 0)
324 ktell = 1;
325 else if (curve->et[kr] > tlast)
326 ktell++;
327 tlast = curve->et[kr];
328 tinmin = MIN (tinmin, tlast);
329 tinmax = MAX (tinmax, tlast);
330 }
331 kmax = MAX (kmax, ktell);
332 }
333
334 /* kmax now contains the number of distinct knot values
335 * to be kept in the interval [tinmin,tinmax]. */
336
337 if (kmax > 1)
338 tincre = (tinmax - tinmin) / (double) (kmax - 1);
339 else
340 tincre = (double) 0.0;
341 tincr2 = tincre / (double) 2.0;
342
343 /* The values used will be
344 tinmin,tinmin+tincre,...,tinmin+(kmax-1)*tincre.
345 We want to use the values closest to the actual
346 knot values when possible.
347 Run through each curve and move knots if tinmin<tinmax.
348 If they are equal, knots should not be moved. */
349
350 if (tinmin >= tinmax)
351 continue;
352
353 for (kj = 0; kj < inbcrv; kj++)
354 {
355 curve = crvarr[kj];
356 ktell = 0;
357 for (kr = *iordr - 1; kr <= curve->in; kr++)
358 if ((tinmin <= curve->et[kr]) && (curve->et[kr] <= tinmax))
359 {
360 if (ktell == 0)
361 {
362 ktell = 1;
363 kfirst = kr;
364 }
365 else if (curve->et[kr] >= tlast)
366 {
367 ktell++;
368 klast = kr;
369 }
370 tlast = curve->et[kr];
371 }
372
373 /* ktell contains the number of knots in [tinmin,tinmax]
374 on curve kj. kfirst is a pointer to first knot value,
375 klast is the pointer to last knot value */
376
377 if (ktell == 1)
378 {
379
380 /* Only one knot value in interval, move to closest
381 * legal value */
382
383 if ((kmax == 1) || (tincre == (double) 0.0))
384 tval = tinmin;
385 else
386 {
387 kdum1 = (int)((curve->et[kfirst] - tinmin + tincr2) / tincre);
388 tval = tinmin + (double)kdum1 * tincre;
389 }
390 tlast = curve->et[kfirst];
391 do
392 {
393 curve->et[kfirst] = tval;
394 kfirst++;
395 } while (curve->et[kfirst] == tlast);
396 continue;
397 }
398
399 if ((ktell <= 1) || (tincre <= (double) 0.0))
400 continue;
401
402
403 /* More than one point found. */
404
405 kdum1 = (int)((curve->et[kfirst] - tinmin + tincr2) / tincre);
406 kdum2 = (int)((curve->et[klast] - tinmin + tincr2) / tincre);
407 kdiff = ktell - kdum2 + kdum1 - 1;
408 if (kdiff > 0)
409 {
410 /* Change kdum1 and kdum2 to allow for ktell
411 different values. */
412
413 kdf = kdiff / 2;
414 kdum1 = MAX (0, kdum1 - kdf);
415 kdum2 = MIN (kmax - 1, kdum2 + kdiff - kdf);
416 kdiff = ktell - kdum2 + kdum1 - 1;
417
418 if (kdiff > 0)
419 {
420 if (kdum1 == 0)
421 kdum2 += kdiff;
422 else
423 {
424 if (kdum2 != ktell)
425 goto err170;
426 kdum1 -= kdiff;
427 }
428 }
429 }
430
431 if ((kdum1 < 0) || (kdum2 > ktell))
432 goto err170;
433
434
435 /* Use kdum1 as start value, move knots. */
436
437 tval = tinmin + kdum1 * tincre;
438 kr = kfirst;
439 tlast = curve->et[kr];
440 for (kp = 0; kp < ktell; kp++, kr++)
441 if (curve->et[kr] == tlast)
442 curve->et[kr] = tval;
443 else
444 {
445 if (curve->et[kr] <= tlast)
446 goto err170;
447 tval += tincre;
448 tlast = curve->et[kr];
449 curve->et[kr] = tval;
450 }
451 }
452 }
453
454 /* Make start knot vector */
455
456 if (kt != SISL_NULL) freearray(kt); /* PFU 18/07-94 */
457
458 kt = newarray (*iordr * 2, double);
459 if (kt == SISL_NULL)
460 goto err101;
461
462 for (ki = 0; ki < *iordr; ki++)
463 {
464 kt[ki] = start;
465 kt[*iordr + ki] = stop;
466 }
467
468 kn = *iordr;
469
470 /* Run through all knot vectors, lift the order, map into right
471 interval and make the union with the candidates already existing. */
472
473 for (ki = 0; ki < inbcrv; ki++)
474 {
475 /* Increase order of basis. */
476
477 curve = crvarr[ki];
478 s1754 (curve->et, curve->in, curve->ik, *iordr, &incknt, &knumb, &kstat);
479 if (kstat < 0)
480 goto error;
481
482
483 /* Make union of old union and new knot vector */
484
485 s1935 (kt, kn, incknt, knumb, &knot, &kn, *iordr, &kstat);
486 if (kstat < 0)
487 goto error;
488
489 if (incknt != SISL_NULL) freearray(incknt); /* PFU 18/07-94 */
490
491 if (kt != SISL_NULL)
492 freearray (kt);
493 kt = knot;
494 knot = SISL_NULL; /* PFU 18/07-94 */
495 }
496
497 /* No errors */
498
499 *in = kn;
500 *it = kt;
501 kt = SISL_NULL; /* PFU 18/07-94 */
502 goto out;
503
504
505 /* Memory error */
506
507 err101:
508 *jstat = -101;
509 s6err ("s1933", *jstat, kpos);
510 goto out;
511
512 /* Error in input */
513
514 err112:
515 *jstat = -112;
516 s6err ("s1933", *jstat, kpos);
517 goto out;
518
519 /* Too few curves for spline lofted surface */
520
521 err179:
522 *jstat = -179;
523 s6err ("s1933", *jstat, kpos);
524 goto out;
525
526 /* Error in lower level routines */
527
528 error:
529 *jstat = kstat;
530 s6err ("s1933", *jstat, kpos);
531 goto out;
532
533 /* Special error in moving of knot values */
534
535 err170:
536 *jstat = -170;
537 s6err ("s1933", *jstat, kpos);
538 goto out;
539
540 out:
541 if (kt != SISL_NULL) freearray(kt); /* PFU 18/07-94 */
542 if (knot != SISL_NULL) freearray(knot); /* PFU 18/07-94 */
543 if (incknt != SISL_NULL) freearray(incknt); /* PFU 18/07-94 */
544 return;
545 }
546