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 #define S1965
44
45 #include "sislP.h"
46
47
48 #if defined(SISLNEEDPROTOTYPES)
s1965(SISLSurf * oldsurf,double eps[],int edgefix[4],int iopen1,int iopen2,double edgeps[],int opt,int itmax,SISLSurf ** newsurf,double maxerr[],int * stat)49 void s1965(SISLSurf *oldsurf,double eps[],int edgefix[4],int iopen1,
50 int iopen2,double edgeps[],int opt,int itmax,
51 SISLSurf **newsurf,double maxerr[],int *stat)
52 #else
53 void s1965(oldsurf, eps, edgefix, iopen1, iopen2, edgeps, opt,
54 itmax, newsurf,maxerr, stat)
55 SISLSurf *oldsurf;
56 double eps[];
57 int edgefix[4];
58 int iopen1;
59 int iopen2;
60 double edgeps[];
61 int opt;
62 int itmax;
63 SISLSurf **newsurf;
64 double maxerr[];
65 int *stat;
66 #endif
67 /*
68 *********************************************************************
69 *
70 *********************************************************************
71 *
72 * Purpose: To remove as many knots as possible from the surface oldsurf
73 * without perturbing the surface more than the given tolerance.
74 * The tolerances are given by eps and epsfix, while the
75 * approximation is given by newsurf.
76 * Mark that the error in continuity over the start and end of
77 * a closed or periodic surface is only guaranteed to be within
78 * edgeps.
79 *
80 *
81 *
82 * Input: oldsurf - pointer to the original spline surface. Note
83 * if the polynomial orders of the surface are
84 * k1 and k2, then the two knot vectors are
85 * assumed to have knots of multiplicity k1 and
86 * k2 at the ends.
87 *
88 * eps - double array of length dim (the number of
89 * components of the surface, typically three)
90 * giving the desired accuracy of the
91 * final approximation compared to oldcurve.
92 * Note that in such comparisons the two
93 * surfaces are not reparametrized in any way.
94 *
95 * edgefix - integer array of dimension (4) giving the number
96 * of derivatives to be kept fixed along each edge
97 * of the surface. The numbering of the edges is the
98 * same as for edgeps below. All the derivatives of
99 * order < nend(i)-1 will be kept fixed along
100 * edge i. Hence nend(i)=0 indicates that nothing is
101 * to be kept fixed along edge i.
102 *
103 * NB! TO BE KEPT FIXED HERE MEANS TO HAVE ERROR LESS THAN
104 * EDGEPS. IN GENERAL, IT IS IMPOSSIBLE TO REMOVE KNOTS
105 * AND KEEP AN EDGE COMPLETELY FIXED.
106 *
107 * iopen1 - Open/closed parameter in first parameter direction.
108 * = 1 : Produce open surface.
109 * = 0 : Produce closed, non-periodic surface if possible.
110 * = -1 : Produce closed, periodic surface
111 *
112 * iopen2 - Open/closed parameter in second parameter direction.
113 * = 1 : Produce open surface.
114 * = 0 : Produce closed, non-periodic surface if possible.
115 * = -1 : Produce closed, periodic surface
116 *
117 * edgeps - double array of length 4*dim ([4,dim]) (dim is
118 * the number of components of each coefficient)
119 * containing the maximum deviation which is
120 * acceptable along the edges of the surface.
121 * edgeps[0]-edgeps[dim-1] gives the tolerance along
122 * the edge corresponding to x1 (the first parameter)
123 * having it's minimum value.
124 * edgeps[dim]-edgeps[2*dim-1] gives the tolerance
125 * along the edge corresponding to x1 (the first
126 * parameter) having it's maximum value.
127 * edgeps[2*dim]-edgeps[3*dim-1] gives the tolerance
128 * along the edge corresponding to x2 (the second
129 * parameter) having it's minimum value.
130 * edgeps[3*dim]-edgeps[4*dim-1] gives the tolerance
131 * along the edge corresponding to x2 (the second
132 * parameter) having it's maximum value.
133 * NB! EDGEPS WILL ONLY HAVE ANY SIGNIFICANCE IF THE
134 * CORRESPONDING ELEMENT OF EDGEFIX IS POSITIVE.
135 *
136 *
137 * itmax - maximum number of iterations. The routine will
138 * follow an iterative procedure trying to remove
139 * more and more knots, one direction at a time.
140 * The process will almost always stop after less
141 * than 10 iterations and it will often stop after
142 * less than 5 iterations. A suitable value for itmax
143 * is therefore usually in the region 3-10.
144 *
145 *
146 * opt - integer indicating the order in which the
147 * knot removal is to be performed.
148 * = 1 - remove knots in parameter direction 1 only.
149 * = 2 - remove knots in parameter direction 2 only.
150 * = 3 - remove knots first in parameter direction
151 * 1 and then 2.
152 * = 4 - remove knots first in parameter direction
153 * 2 and then 1.
154 *
155 *
156 *
157 *
158 *
159 * Output:
160 * newsurf - the approximating surface on the reduced knot vectors.
161 *
162 * maxerr - double array of length dim containing an upper
163 * bound for the pointwise error in each of the
164 * components of the spline approximation. The two
165 * surfaces oldsurf and newsurf are compared at the
166 * same parameter vaues, i.e., if oldsurf is f and
167 * newsurf is g then
168 * |f(u,v)-g(u,v)| <= eps
169 * in each of the components.
170 *
171 * stat - status message
172 * > 0 : warning
173 * = 0 : ok
174 * < 0 : error
175 *
176 *
177 *
178 *
179 * Method:
180 * The data reduction is performed in one parameter direction at a time,
181 * using a subroutine for data reduction for spline curves (s1940).
182 * Consider for example knot removal in the second parameter direction
183 * (suppose oldcurve has m1xm2 coefficients, each with dim components
184 * and suppose also that the knot vectors are t1 and t2).
185 * This amounts to removing knots from each of the m1*dim curves, each
186 * of length m2 and with knot vector t2, obtained by considering the
187 * coefficients d as an array of dimension [m2][dim*m1] instead
188 * of [m2][m1][dim] (note that in the code the arrays are always
189 * treated as one dimensional).
190 * A similar approach is possible when removing knots in the first parameter
191 * direction except that then d has to be transposed in the first two
192 * dimensions first.
193 * The remaining part of the subroutine is more or less obvious except
194 * perhaps the fact that if knots are to be removed in both directions,
195 * the permissible error must be halved so that the total error committed
196 * is kept within the original tolerance.
197 *
198 *
199 * References:
200 * 1. A Data-Reduction Strategy for Splines with Applications to the
201 * Approximation of Functions and data, IMA J. of Num. Anal. 8 (1988),
202 * pp. 185-208.
203 *
204 * 2. Knot Removal for Parametric B-spline Curves and Surfaces,
205 * CAGD 4 (1987), pp. 217-230.
206 *
207 *-
208 * Calls: s1940, s1352, s6chpar, s6err
209 *
210 * Written by: Knut Moerken, University of Oslo, December 1992, based
211 * on an earlier Fortran version by the same author.
212 * Changed by: Paal Fugelli, SINTEF, 1994-07.
213 * Rewritten to fix several major memory leakage problems.
214 * Changed and renamed by : Vibeke Skytt, SINTEF Oslo, 02.95. Introduced
215 * periodicity.
216 *
217 *********************************************************************
218 */
219
220 {
221 int k1 = oldsurf->ik1; /* Unwrap some of the most commonly */
222 int k2 = oldsurf->ik2; /* used parameters. */
223 int m1 = oldsurf->in1;
224 int m2 = oldsurf->in2;
225 int dim = oldsurf->idim;
226 int kopen1, kopen2;
227 double *et1 = SISL_NULL;
228 double *et2 = SISL_NULL;
229
230 double *local_eps = SISL_NULL; /* Declaration of local variables, */
231 double *local_edge_eps = SISL_NULL; /* for usage, see the code. */
232 double *clocal_eps = SISL_NULL;
233 double *clocal_err = SISL_NULL;
234 double *harray = SISL_NULL;
235 double *tempcoef = SISL_NULL;
236
237 double *loct1 = SISL_NULL;
238 double *loct2 = SISL_NULL;
239 SISLCurve *local_curve = SISL_NULL;
240 SISLCurve *newlcurve = SISL_NULL;
241
242 int lopt, i, it, n1, n2, antit, j, jh, lstat;
243 double factor;
244 int pos = 0; /* Parameter in s6err. */
245
246
247 /* Initialize kopen1 and kopen2 to make sure that the returned
248 surface has legal open/closed parametere also if datareduction
249 only is performed in one direction. */
250
251 kopen1 = oldsurf->cuopen_1;
252 kopen2 = oldsurf->cuopen_2;
253
254 /* Make sure the returned surface has a valid value. */
255 (*newsurf) = SISL_NULL;
256
257 /* Unwrap the two (possibly changed) input knot vectors */
258
259 et1 = oldsurf->et1;
260 et2 = oldsurf->et2;
261
262
263 /* We start by allocating space for local arrays (any error exit must
264 now handle free'ing of allocated space). */
265
266 local_eps = newarray(dim, DOUBLE);
267 if (local_eps == SISL_NULL) goto err101;
268
269 local_edge_eps = newarray(4*dim, DOUBLE);
270 if (local_edge_eps == SISL_NULL) goto err101;
271
272 /* If we are going to remove knots in both directions, we use half the
273 tolerance each time, otherwise the whole tolerance once. This is
274 signified by factor. */
275
276 factor = (opt < 3) ? 1.0 : 0.5;
277
278 /* Store a local version of the tolerance and initialize the error
279 to zero. */
280
281 for (i=0; i<dim; i++)
282 {
283 local_eps[i] = factor*eps[i];
284 maxerr[i] = 0.0;
285 }
286
287 /* We also need a local version of the edge tolerance. */
288
289 /* Originally the edgetolerance was also multiplied by fac. It turns
290 out that this is unnecessary: Each edge is completely fixed in one
291 of the calls to s1940. */
292
293 for (i=0; i<4*dim; i++) local_edge_eps[i] = edgeps[i];
294
295 /* Now we are ready to prepare for the iterations. it counts the iterations
296 while lopt says which direction to start with and antit tells us how
297 many iterations there should be (one or two). n1 and n2 are used for
298 keeping track of the number of coefficients. */
299
300 it = 0; lopt = (opt-1) % 2 + 1;
301 n1 = m1; n2 = m2;
302
303 antit = (opt < 3) ? 1 : 2;
304 while (it < antit)
305 {
306 if (lopt == 1)
307 {
308
309 /* Here we remove knots along the first parameter direction. */
310
311 /* Since we are going to treat the surface as a collection of
312 n2 curves in the first direction each with dim components,
313 we need a tolerance array of length dim*n2 and a similar
314 array for storing the error. */
315
316 clocal_eps = newarray(dim*n2, DOUBLE);
317 if (clocal_eps == SISL_NULL) goto err101;
318
319 clocal_err = newarray(dim*n2, DOUBLE);
320 if (clocal_err == SISL_NULL) goto err101;
321
322 /* The tolerance will in general not be the same for all the
323 n2*dim curves. If edges 1 and/or 3 are to be kept fixed then
324 first and/or last tolerance must be smaller than the tolerance
325 for middle curves. This is accomplished by a call to s1352. */
326
327 s1352(et2, n2, k2, local_eps, (local_edge_eps+2*dim),
328 (local_edge_eps+3*dim), dim, edgefix[2], edgefix[3],
329 clocal_eps, &lstat);
330 if (lstat < 0) goto error;
331
332 /* First we have to pick out the correct coefficients.
333 Then we must transpose the coefficient with respect to
334 n1 and n2 due to the way we store the coefficients of
335 surfaces. */
336
337 harray = newarray(n1*n2*dim, DOUBLE);
338 if (harray == SISL_NULL) goto err101;
339
340 if (opt == 4)
341 {
342 /* Here tempcoef points to the coefs of the first direction. */
343
344 s6chpar(tempcoef, n1, n2, dim, harray);
345 freearray(tempcoef);
346 tempcoef = harray;
347 }
348 else
349 {
350 s6chpar(oldsurf->ecoef, n1, n2, dim, harray);
351 tempcoef = harray;
352 }
353 harray = SISL_NULL;
354
355 /* We then create a curve in which to store the high dimensional
356 curve. */
357
358 local_curve = newCurve(n1, k1, oldsurf->et1, tempcoef,
359 1, dim*n2, 0);
360 if (local_curve == SISL_NULL) goto err101;
361 local_curve -> cuopen = oldsurf->cuopen_1;
362
363 /* We can now perform knot removal on this curve. The result
364 is stored in newlcurve (will always have icopy==1, i.e. a
365 proper copy). */
366
367 s1940(local_curve, clocal_eps, edgefix[0], edgefix[1], iopen1,
368 itmax, &newlcurve, clocal_err, &lstat);
369 if(lstat<0) goto error;
370
371 kopen1 = newlcurve->cuopen;
372
373 /* Remember to update n1 and throw away local_curve (has icopy==0). */
374
375 n1 = newlcurve->in;
376 freeCurve(local_curve);
377 local_curve = SISL_NULL;
378
379 /* Now we must transpose back to return to surface coeffiecients. */
380
381 tempcoef = increasearray(tempcoef,n1*n2*dim, DOUBLE);
382 if (tempcoef == SISL_NULL) goto err101;
383
384 harray = newlcurve->ecoef;
385 s6chpar(harray, n2, n1, dim, tempcoef);
386 harray = SISL_NULL;
387
388 /* We save the surface as a curve and create the final surface
389 later. */
390
391 local_curve = newlcurve;
392 freearray(newlcurve->ecoef);
393 local_curve->ecoef = tempcoef;
394 newlcurve = SISL_NULL;
395
396 /* Now local_curve has icopy==1, so the knots and coefs are proper
397 copies. Must make it safe to free local_curve. */
398
399 local_curve->icopy = 0;
400
401 /* Now we can just save the new knot vector in loct1 and update the
402 pointer to the "input" knot vector et1. */
403
404 loct1 = local_curve->et;
405 et1 = loct1;
406
407 /* tempcoef and loct1 are now safe from free'ing of local_curve, so
408 throw it away. */
409
410 freeCurve(local_curve);
411 local_curve = SISL_NULL;
412
413 /* Calculate the error in the approximation. */
414
415 /* The error in component i in the approximation is simply
416 the largest error in component i in any of the curves. */
417
418 for (i=0; i<dim; i++) clocal_eps[i] = 0.0;
419 for (jh=0, i=0; i<n2; i++)
420 for (j=0; j<dim; j++, jh++)
421 clocal_eps[j] = MAX(clocal_eps[j], clocal_err[jh]);
422
423 /* Accumulate the error in maxerr and use the remaining tolerance
424 in the next direction. */
425
426 for (i=0; i<dim; i++)
427 {
428 maxerr[i] += clocal_eps[i];
429 local_eps[i] = eps[i] - maxerr[i];
430 }
431 freearray(clocal_eps);
432 freearray(clocal_err);
433 }
434 else
435 {
436
437
438 /* Now we are to remove knots in the second direction. This is
439 very similar to removing knots in the first direction but
440 we do not need to transpose. */
441
442 /* First some local arrays. */
443
444 clocal_eps = newarray(dim*n1, DOUBLE);
445 if (clocal_eps == SISL_NULL) goto err101;
446
447 clocal_err = newarray(dim*n1, DOUBLE);
448 if (clocal_err == SISL_NULL) goto err101;
449
450 /* Compute a tolerance along the second direction. */
451
452 s1352(et1, n1, k1, local_eps, local_edge_eps,
453 local_edge_eps+dim, dim, edgefix[0], edgefix[1],
454 clocal_eps, &lstat);
455 if (lstat < 0) goto error;
456
457 /* Generate a high dimensional curve along the second direction
458 with the right coefficients. */
459
460 if (opt == 3)
461 {
462 /* Here tempcoef points to the coefs of the first direction. */
463 local_curve = newCurve(n2, k2, oldsurf->et2, tempcoef,
464 1, dim*n1, 0);
465 }
466 else
467 local_curve = newCurve(n2, k2, oldsurf->et2, oldsurf->ecoef, 1,
468 dim*n1, 0);
469 if (local_curve == SISL_NULL) goto err101;
470 local_curve->cuopen = oldsurf->cuopen_2;
471
472 /* Remove knots and store in newlcurve (will always get icopy==1, i.e.
473 proper copy. */
474
475 s1940(local_curve, clocal_eps, edgefix[2], edgefix[3], iopen2,
476 itmax, &newlcurve, clocal_err, &lstat);
477 if(lstat<0) goto error;
478
479 kopen2 = newlcurve->cuopen;
480
481 /* Remember to update n2. */
482
483 n2 = newlcurve->in;
484
485 /* Throw away the old local_curve which has icopy==0 (must throw away
486 tempcoef too -- non-null if opt==3), but keep newlcurve. */
487
488 freeCurve(local_curve);
489 if (tempcoef != SISL_NULL) freearray(tempcoef);
490 local_curve = newlcurve;
491 newlcurve = SISL_NULL;
492
493 /* Now local_curve has icopy==1, so the knots and coefs are proper
494 copies. Must make it safe to free local_curve. */
495
496 local_curve->icopy = 0;
497 tempcoef = local_curve->ecoef;
498
499 /* The new knot vector in the second direction will be stored in
500 loct2 so we must update the pointer to the "input" knot vector et2. */
501
502 loct2 = local_curve->et;
503 et2 = loct2;
504
505 /* tempcoef and loct2 are now safe from free'ing of local_curve, so
506 throw it away. */
507
508 freeCurve(local_curve);
509 local_curve = SISL_NULL;
510
511 /* Calculate the error in the approximation as above. */
512
513 for (i=0; i<dim; i++) clocal_eps[i] = 0.0;
514 for (jh=0, i=0; i<n1; i++)
515 for (j=0; j<dim; j++, jh++)
516 clocal_eps[j] = MAX(clocal_eps[j], clocal_err[jh]);
517
518 for (i=0; i<dim; i++)
519 {
520 maxerr[i] += clocal_eps[i];
521 local_eps[i] = eps[i] - maxerr[i];
522 }
523 freearray(clocal_eps);
524 freearray(clocal_err);
525 }
526
527 /* Prepare for the next iteration. */
528
529 ++it;
530 lopt = 3 - lopt;
531
532 }
533
534 /* It remains to create a new surface object. If we only removed knots
535 in one direction, one of loct1 and loct2 will be SISL_NULL. */
536
537 if (loct1 == SISL_NULL)
538 {
539 loct1 = newarray(n1+k1, DOUBLE);
540 if (loct1 == SISL_NULL) goto err101;
541
542 harray = et1;
543 for (i=0; i<n1+k1; i++) loct1[i] = harray[i];
544 }
545
546 if (loct2 ==SISL_NULL)
547 {
548 loct2 = newarray(n2+k2, DOUBLE);
549 if (loct2 == SISL_NULL) goto err101;
550
551 harray = et2;
552 for (i=0; i<n2+k2; i++) loct2[i] = harray[i];
553 }
554
555 /* Generate the new surface object. */
556
557 *newsurf = newSurf(n1, n2, k1, k2, loct1, loct2, tempcoef,
558 1, dim, 2);
559 if (*newsurf == SISL_NULL) goto err101;
560
561 /* Avoid free'ing the referenced knots and coefs on exit. */
562
563 loct1 = SISL_NULL;
564 loct2 = SISL_NULL;
565 tempcoef = SISL_NULL;
566
567 /* Set periodicity flag. */
568
569 (*newsurf)->cuopen_1 = kopen1;
570 (*newsurf)->cuopen_2 = kopen2;
571
572 *stat = 0;
573 goto out;
574
575 /* Error in memory allocation. */
576
577 err101:
578 *stat = -101;
579 s6err("s1965", *stat, pos);
580 goto out;
581
582 /* Error in lower level routine. */
583
584 error:
585 *stat = lstat;
586 s6err("s1965", *stat, pos);
587
588 /* Clean up. */
589
590 out:
591 if (local_eps != SISL_NULL) freearray(local_eps);
592 if (local_edge_eps != SISL_NULL) freearray(local_edge_eps);
593 if (clocal_eps != SISL_NULL) freearray(clocal_eps);
594
595 /* Must remember to free everything to avoid memory leak. */
596
597 if (clocal_eps != SISL_NULL) freearray(clocal_eps);
598 if (clocal_err != SISL_NULL) freearray(clocal_err);
599 if (tempcoef != SISL_NULL) freearray(tempcoef);
600 if (loct1 != SISL_NULL) freearray(loct1);
601 if (loct2 != SISL_NULL) freearray(loct2);
602 if (local_curve != SISL_NULL) freeCurve(local_curve); /* icopy==0 */
603 if (newlcurve != SISL_NULL) freeCurve(newlcurve); /* icopy ==1 */
604
605
606 return;
607
608 }
609