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: s1340.c,v 1.5 2001-03-19 15:58:46 afr Exp $
45 *
46 */
47
48 #define S1340
49
50 #include "sislP.h"
51
52 #if defined(SISLNEEDPROTOTYPES)
53 void
s1340(SISLCurve * oldcurve,double eps[],int startfix,int endfix,double epsco,int itmax,SISLCurve ** newcurve,double maxerr[],int * stat)54 s1340(SISLCurve *oldcurve, double eps[], int startfix, int endfix,
55 double epsco, int itmax, SISLCurve **newcurve, double maxerr[],
56 int *stat)
57 #else
58 void
59 s1340(oldcurve, eps, startfix, endfix, epsco, itmax, newcurve,
60 maxerr, stat)
61 SISLCurve *oldcurve;
62 double eps[];
63 int startfix;
64 int endfix;
65 double epsco;
66 int itmax;
67 SISLCurve **newcurve;
68 double maxerr[];
69 int *stat;
70 #endif
71 /*
72 *********************************************************************
73 *
74 *********************************************************************
75 *
76 * Purpose: To remove as many knots as possible from the spline curve
77 * "oldcurve" without perturbing the curve more than a given
78 * tolerance. The tolerance is given by "eps" and the
79 * approximation by "newcurve".
80 *
81 *
82 * Input:
83 * oldcurve - pointer to the original spline curve. Note that
84 * if the polynomial order of the curve is k, then
85 * the knot vector of oldcurve is assumed to have
86 * knots of multiplicity k at the beginning and end.
87 *
88 * eps - double array giving the desired absolute accuracy
89 * of the final approximation as compared to oldcurve.
90 * If oldcurve is a spline curve in a space of
91 * dimension dim, then eps must have length dim.
92 * Note that it is not relative, but absolute accuracy
93 * that is being used. This means that the difference
94 * in component i at any parameter value, between
95 * the given curve and the approximation, is to be
96 * less than eps[i]. Note that in such comparisons
97 * the same parametrization is used for both curves.
98 *
99 * startfix - the number of derivatives to be kept fixed at
100 * the beginning of the knot interval.
101 * The (0 - (startfix-1)) derivatives will be kept fixed.
102 * If startfix <0, this routine will set it to 0.
103 * If startfix< the order of the curve, this routine
104 * will set it to the order.
105 *
106 * endfix - the number of derivatives to be kept fixed at
107 * the end of the knot interval.
108 * The (0 - (endfix-1)) derivatives will be kept fixed.
109 * If endfix <0, this routine will set it to 0.
110 * If endfix< the order of the curve, this routine
111 * will set it to the order.
112 *
113 * epsco - real number used to check equality between real
114 * numbers in some routines. Two numbers differing
115 * by a relative amount less than epsco will in some
116 * cases be considered equal. A suitable value is just
117 * above the unit roundoff of the machine. In IEEE
118 * double precision arithmetic 1.0E-15 is a reasonable
119 * choice.
120 * NB! The computations are not guaranteed to have
121 * relative accuracy less than epsco.
122 *
123 * itmax - maximum number of iterations. The routine will
124 * follow an iterative procedure trying to remove
125 * more and more knots. The process will almost always
126 * stop after less than 10 iterations and it will often
127 * stop after less than 5 iterations. A suitable
128 * value for itmax is therefore usually in the region
129 * 3-10.
130 *
131 *
132 *
133 *
134 * Output:
135 * newcurve - the spline approximation on the reduced
136 * knot vector.
137 * maxerr - double array containing an upper bound for the
138 * pointwise error in each of the components of the
139 * spline approximation. The two curves oldcurve and
140 * newcurve are compared at the same parameter value,
141 * i.e., if oldcurve is f and newcurve is g, then
142 * |f(t)-g(t)| <= eps
143 * in each of the components.
144 *
145 * stat - status message
146 * > 0 : warning
147 * = 0 : ok
148 * < 0 : error
149 *
150 *
151 *
152 *
153 * Method:
154 * The method is described in two papers listed as references below.
155 * First, the knots of the input spline, s, are ranked according to their
156 * relative importance in the representation of s (s1353). Using this
157 * information, knots are removed from s and an approximation g is computed
158 * on the reduced knot vector (s1354) with the error guaranteed to be
159 * less than the tolerance in all components. In order to remove more knots,
160 * the knots of g are ranked, and then knots are removed from g, resulting
161 * in a spline h with even fewer knots. The knot vector of h is then
162 * clearly a subsequence of the knot vector of s and s is approximated
163 * on this knot vector (sh1365). If the error is less than the tolerance
164 * then h is accepted as the new best approximation and stored in g.
165 * If the error is greater than the tolerance, knots are removed directly
166 * from s based on the ranking of g resulting in a new g. In either case a
167 * new best approximation g is obtained, and we can continue the iteration.
168 * The iterations stop when no more knots can be removed or all the interior
169 * knots have been removed. (The names s, g, and h are not used in the
170 * code.)
171 *
172 * NB! It is assumed that the input spline has knots of multiplicity k (the
173 * polynomial order) at each end. This property of the knot vector is
174 * used extensively by the lower level routines.
175 *
176 *
177 * References:
178 * 1. A Data-Reduction Strategy for Splines with Applications to the
179 * Approximation of Functions and data, IMA J. of Num. Anal. 8 (1988),
180 * pp. 185-208.
181 *
182 * 2. Knot Removal for Parametric B-spline Curves and Surfaces,
183 * CAGD 4 (1987), pp. 217-230.
184 *
185 *
186 * Calls: s1353, s1354, sh1365, s6err
187 *
188 * Written by: Knut Moerken, University of Oslo, November 1992, based
189 * on an earlier Fortran version by the same author.
190 * CHANGED BY: Paal Fugelli, SINTEF, 1994-07. Initialized pointers (to SISL_NULL)
191 * in 'ranking' to avoid potential memory leak when exiting through 'out'.
192 * Removed several other memory leaks.
193 * CHANGED BY: Per OEyvind Hvidsten, SINTEF, 1994-11. Added a freeCurve
194 * call before overwriting the *newarray pointer (thus removing a
195 * memory leak.
196 *
197 *********************************************************************
198 */
199 {
200 char ready, big; /* Boolean variables used for computing
201 the stopping criteria. */
202 int k = oldcurve->ik; /* Unwrapped version of the given curve. */
203 int dim = oldcurve->idim;
204 double *d = oldcurve->ecoef;
205 int itcount=0; /* Iteration counter. */
206 int n1 = oldcurve->in; /* n1 and n2 are used for storing the number
207 coefficients in consecutive spline
208 approximations. The situation n1=n2
209 signifies that no more knots can be
210 removed. */
211 int n2;
212 int i, mini, maxi, indx; /* Various auxiliary integer variables. */
213 double *local_err = SISL_NULL; /* Variables used for storing local error
214 estimates. */
215 double *l2err = SISL_NULL;
216 double *lepsco = SISL_NULL;
217 double *temp_err = SISL_NULL;
218 SISLCurve *tempcurve = SISL_NULL; /* Variables that are used for storing
219 temporary curves. */
220 SISLCurve *helpcurve = SISL_NULL;
221 rank_info ranking; /* Variable used for holding ranking
222 information (from s1353). */
223 int lstat=0; /* Local status variable. */
224 int pos=0; /* Parameter to s6err. */
225 SISLCurve *qc_kreg = SISL_NULL; /* Non-periodic version of the input curve. */
226
227
228 /* Initialize ranking ptrs in case of early exit through 'out' (PFU 05/07-94) */
229 ranking.prio = SISL_NULL;
230 ranking.groups = SISL_NULL;
231
232 /* Initialize maxerr to zero. */
233
234 for (i=0; i<dim; i++) maxerr[i] = 0.;
235
236 /* Only interior knots may be removed so if n1==k we can stop
237 straight away. */
238
239 if (n1 == k)
240 {
241 *newcurve = newCurve(n1, k, oldcurve->et, oldcurve->ecoef,
242 oldcurve->ikind, oldcurve->idim, 1);
243 if (*newcurve == SISL_NULL) goto err101;
244
245 *stat = 0;
246 goto out;
247 }
248
249 /* Make sure that the input curve is non-periodic. */
250
251 if (oldcurve->cuopen == SISL_CRV_PERIODIC)
252 {
253 make_cv_kreg(oldcurve, &qc_kreg, &lstat);
254 if (lstat < 0) goto error;
255
256 /* The input curve is closed and periodic. Make sure that the
257 endpoints of the curve is still matching by fixing the
258 position and the derivative in the endpoints of the curve.
259 The change made to startfix and endfix is only locallly. */
260
261 startfix = MAX(startfix, 2);
262 endfix = MAX(endfix, 2);
263 }
264 else
265 qc_kreg = oldcurve;
266
267 /* Allocate space for some local arrays. */
268
269 temp_err = newarray(dim, double);
270 if (temp_err == SISL_NULL) goto err101;
271
272 lepsco = newarray(dim, double);
273 if (lepsco == SISL_NULL) goto err101;
274
275 /* ranking is of type rank_info which is a struct described in s1353. */
276
277 ranking.prio = newarray(MAX(n1-k,1), int);
278 if (ranking.prio == SISL_NULL) goto err101;
279
280 ranking.groups = newarray(MAX(n1-k,1), int);
281 if (ranking.groups == SISL_NULL) goto err101;
282
283 /* lespco is needed in s1354. In component i we first store the l1-norm
284 of the component i of the B-spline coefficients of oldcurve. */
285
286 for (i=0; i<dim; i++) lepsco[i] = fabs(d[i]);
287
288 indx = 0;
289 for (i=1; i<n1*dim; i++)
290 {
291 lepsco[indx++] += fabs(d[i]);
292 if (indx == dim) indx = 0;
293 }
294
295 /* We can now compute the final value of lepsco which will be used for
296 checking if two numbers are almost equal. */
297
298 for (i=0; i<dim; i++)
299 lepsco[i] = MIN(lepsco[i]*epsco/n1, eps[i]);
300
301 /* This is where the knot removal process starts. */
302
303 /* mini and maxi are lower and upper bounds on how many knots that can
304 be removed. */
305
306 mini = 0;
307 maxi = n1 - k + 1;
308
309 /* Start by ranking the knots of oldcurve. */
310
311 s1353(qc_kreg, eps, &ranking, &lstat);
312 if (lstat < 0) goto error;
313
314 /* Based on the computed ranking, we remove as close to maxi knots
315 from oldcurve as possible, but such that we can always compute an
316 approximation to oldcurve on the reduced knot vector, with error less
317 than eps. newcurve will be created with icopy==1. */
318
319 s1354(qc_kreg, qc_kreg, &ranking, eps, lepsco, startfix, endfix,
320 mini, maxi, newcurve, maxerr, &lstat);
321 if (lstat < 0) goto error;
322
323 /* The spline stored in newcurve is now an approximation to oldcurve
324 with error less than eps. We will now iterate and try to remove knots
325 from newcurve. The integers n1 and n2 will be the number of knots in
326 the two most recent approximations. */
327
328 n2 = (*newcurve)->in;
329
330 /* Start the iterations. We have already done one iteration and we have
331 not had any problems with the error getting too big. */
332
333 itcount = 1;
334 big = 0;
335
336 /* If n1=n2 we are unable to remove any knots, and if n2=k there are
337 no interior knots left so we can stop. Likewise if itmax was 1. */
338
339 ready = (n1 == n2) || (n2 == k) || (itcount == itmax);
340 while (!ready)
341 {
342
343 /* We start by ranking the knots of newcurve which is the approximation
344 with the fewest knots that we have found so far. */
345
346 s1353(*newcurve, eps, &ranking, &lstat);
347 if (lstat < 0) goto error;
348
349 if (!big)
350 {
351
352 /* If we have not had any problems with the error in approximation
353 becoming too big we take the risk of removing as many knots as
354 we can from newcurve and in this way determining an even shorter
355 knot vector.
356 First we iterate in s1354 to see how many knots we can remove
357 and store the approximation to newcurve in helpcurve (will be
358 created with icopy==1). */
359
360 mini = 0; maxi = (*newcurve)->in - k + 1;
361 s1354(*newcurve, *newcurve, &ranking, eps, lepsco, startfix,
362 endfix, mini, maxi, &helpcurve, temp_err, &lstat);
363 if (lstat < 0) goto error;
364
365 /* Now, we have no guarantee that helpcurve is closer to oldcurve
366 than the tolerance so we compute a new approximation to oldcurve
367 on the knot vector of helpcurve and store this in tempcurve (with
368 icopy==1).
369 Must make sure that local_err is free'ed if allocated from last
370 iteration. */
371
372 if (local_err != SISL_NULL) freearray(local_err);
373
374 sh1365(qc_kreg, helpcurve->et, k, helpcurve->in,
375 startfix, endfix, &tempcurve,
376 &local_err, &l2err, &lstat);
377 if (lstat < 0) goto error;
378
379 /* Don't need l2err or helpcurve anymore. */
380
381 freearray(l2err);
382 freeCurve(helpcurve);
383 helpcurve = SISL_NULL;
384
385 /* We must now check if the new tempcurve is within the tolerance. */
386 i = 0;
387 while (!big && i<dim)
388 {
389 big = local_err[i] > eps[i];
390 i++;
391 }
392 }
393 else
394
395 /* If we have had problems with the error becoming greater than the
396 tolerance, we simply store newcurve in tempcurve and proceed. */
397
398 tempcurve = *newcurve;
399
400 if (big)
401 {
402
403 /* If at some stage we have had problems with the error becoming
404 larger than the tolerance, we do not get involved in the risky
405 approach of removing knots from newcurve. Instead we throw away
406 tempcurve and remove knots directly from oldcurve, but remember
407 that the only knots left to remove are the knots of newcurve for
408 which we use the ranking computed above. The result is stored
409 in helpcurve and later transferred to newcurve. */
410
411 mini= 0; maxi = tempcurve->in - k + 1;
412 s1354(qc_kreg, *newcurve, &ranking, eps, lepsco,
413 startfix, endfix, mini, maxi, &helpcurve, maxerr, &lstat);
414 if (lstat < 0) goto error;
415
416 if (*newcurve != SISL_NULL && *newcurve != tempcurve)
417 {
418 freeCurve(*newcurve);
419 *newcurve = SISL_NULL;
420 }
421
422 if (tempcurve != SISL_NULL)
423 {
424 freeCurve(tempcurve);
425 tempcurve = SISL_NULL;
426 }
427 *newcurve = helpcurve;
428 helpcurve = SISL_NULL;
429 }
430 else
431 {
432
433 /* Since we now know that the difference between oldcurve and tempcurve
434 is within the tolerance, we just have to store the error in maxerr,
435 throw away the old newcurve, and store tempcurve in newcurve. */
436
437 for (i=0; i<dim; i++) maxerr[i] = local_err[i];
438
439 if (*newcurve != SISL_NULL) freeCurve(*newcurve);
440 *newcurve = tempcurve;
441 tempcurve = SISL_NULL;
442 }
443
444 /* Now we must check if it is time to stop. */
445
446 n1 = n2;
447 n2 = (*newcurve)->in;
448 ++itcount;
449 ready = (n1 == n2) || (n2 == k) || (itcount == itmax);
450 }
451
452 /* Set periodicity flag. */
453
454 if (oldcurve->cuopen == SISL_CRV_CLOSED ||
455 oldcurve->cuopen == SISL_CRV_PERIODIC)
456 (*newcurve)->cuopen = SISL_CRV_CLOSED;
457
458 /* Success */
459
460 *stat = 0;
461 goto out;
462
463 /* Error in allocation of memory. */
464
465 err101:
466 *stat = -101;
467 s6err("s1340", *stat, pos);
468 goto out;
469
470 /* Error in lower level routine. */
471
472 error:
473 *stat = lstat;
474 s6err("s1340", *stat, pos);
475 goto out;
476
477 /* Clear up and free memory before exit. */
478
479 out:
480 if (temp_err != SISL_NULL) freearray(temp_err);
481 if (local_err != SISL_NULL) freearray(local_err);
482 if (l2err != SISL_NULL) freearray(l2err);
483 if (lepsco != SISL_NULL) freearray(lepsco);
484 if (ranking.prio != SISL_NULL) freearray(ranking.prio);
485 if (ranking.groups != SISL_NULL) freearray(ranking.groups);
486 if (qc_kreg != SISL_NULL && qc_kreg != oldcurve) freeCurve(qc_kreg);
487 if (helpcurve != SISL_NULL && helpcurve != (*newcurve)) freeCurve(helpcurve);
488 if (tempcurve != SISL_NULL && tempcurve != (*newcurve)) freeCurve(tempcurve);
489
490 return;
491 }
492