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: s1354.c,v 1.3 2001-03-19 15:58:46 afr Exp $
45 *
46 */
47
48 #define S1354
49
50 #include "sislP.h"
51
52 #if defined(SISLNEEDPROTOTYPES)
s1354(SISLCurve * oldcurve,SISLCurve * rankcurve,rank_info * ranking,double eps[],double epsco[],int startfix,int endfix,int mini,int maxi,SISLCurve ** newcurve,double maxerr[],int * stat)53 void s1354(SISLCurve *oldcurve,SISLCurve *rankcurve,rank_info *ranking,
54 double eps[],double epsco[],int startfix,int endfix,int mini,
55 int maxi,SISLCurve **newcurve,double maxerr[],int *stat)
56 #else
57 void s1354(oldcurve, rankcurve, ranking, eps, epsco,
58 startfix, endfix, mini, maxi, newcurve, maxerr, stat)
59 SISLCurve *oldcurve;
60 SISLCurve *rankcurve;
61 rank_info *ranking;
62 double eps[];
63 double epsco[];
64 int startfix;
65 int endfix;
66 int mini;
67 int maxi;
68 SISLCurve **newcurve;
69 double maxerr[];
70 int *stat;
71 #endif
72 /*
73 *********************************************************************
74 *
75 *********************************************************************
76 *
77 * Purpose: To remove as many knots as possible from oldcurve
78 * without perturbing this spline more than eps, using the ranking
79 * information in rank_info which was obtained from the spline
80 * rankcurve; and compute an approximation to oldcurve on the
81 * knot vector of rankcurve. The knot vector of rankcurve must be
82 * a subsequence of the knot vector of oldcurve and the knots are
83 * always removed from rankcurve so that the final reduced knot
84 * vector will always be a subsequence of both the knot vectors of
85 * oldcurve and rankcurve. The final approximation is given in
86 * newcurve.
87 *
88 *
89 *
90 *
91 * Input:
92 * oldcurve - pointer to the original spline from which the knots
93 * are to be removed.
94 *
95 * rankcurve - pointer to the curve that the ranking information
96 * is based on. It is assumed the knot vector of
97 * rankcurve is a subsequence of the knot vector of
98 * oldcurve.
99 *
100 * ranking - a pointer to a rank_info object containing the
101 * result of the ranking computations. The struct
102 * rank_info contains four variables listed below
103 * (t, k, and n refer to the knot vector of rankcurve,
104 * its polynomial order and the number of B-spline
105 * coefficients):
106 * ranking->prio - integer array of dimension (n-k) containing the
107 * ranking of the interior knots. The knots are listed
108 * in order of increasing ranking number, cf. the
109 * second reference, with knots with the same ranking
110 * number listed in the order in which they occur
111 * in t. To differentiate between the different
112 * ranking numbers the array ranking->groups is used.
113 * ranking->groups - integer array of dimension (n-k) used to partition
114 * ranking->prio into blocks of knots with the same
115 * ranking number. ranking->groups[0] contains the
116 * number of knots with the smallest ranking number,
117 * ranking->groups[1] contains the number of knots with
118 * the two smallest ranking numbers etc.. Only the first
119 * ranking->antgr elements of ranking->groups are used.
120 * ranking->antgr - the number of distinct ranking numbers among the
121 * interior knots of t.
122 * ranking->antrem - an estimate of the number of knots that can be
123 * removed from curve.
124 *
125 * eps - double array containing the tolerance. The
126 * approximation must at each parameter value deviate
127 * less than eps from oldcurve (in each component).
128 *
129 * epsco - double array of containing a tolerance which
130 * is used when the number of coefficients in the
131 * approximating spline becomes very small. More
132 * specifically, if so many knots are removed that
133 * the corresponding spline approximation has fewer
134 * knots than startfix+endfix (the total number of
135 * constraints), then it is required that the error
136 * is less than epsco in each component.
137 *
138 * startfix - the number of derivatives to be kept fixed at the
139 * left end of the parameter interval.
140 *
141 * endfix - the number of derivatives to be kept fixed at the
142 * right end of the parameter interval.
143 *
144 * mini - a lower bound on how many knots that can be removed.
145 *
146 * maxi - an upper bound on how many knots that can be removed.
147 *
148 *
149 * Output:
150 * newcurve - pointer to the spline approximation on the reduced
151 * knot vector.
152 * maxerr - double array containing an upper bound for the
153 * pointwise error in each of the components of the
154 * spline approximation. The two curves oldcurve and
155 * newcurve are compared at the same parameter value,
156 * i.e., if oldcurve is f and newcurve is g, then
157 * |f(t)-g(t)| <= eps
158 * in each of the components.
159 *
160 * stat - status message
161 * > 0 : warning
162 * = 0 : ok
163 * < 0 : error
164 *
165 *
166 * Method:
167 * The routine determines the maximum number of knots that can be removed
168 * by binary search starting by removing ranking->antrem (antrem below)
169 * knots. In other words, it first tries to remove antrem knots from
170 * oldcurve. If this gives an error less than the tolerance, it tries
171 * to remove (antrem+maxi)/2 knots, while if the error is greater than
172 * the tolerance, the routine will try to remove (mini+antrem)/2 of the
173 * knots. This is continued until the number of knots that can be removed
174 * has been determined exactly, together with the corresponding
175 * approximation. The acceptable approximation with the fewest knots is
176 * kept in the variables igtau,igc,ign, and this is
177 * initialized to the spline given by etprio,edprio,imprio.
178 *
179 * During each iteration it must be determined which knots to remove,
180 * given the number of knots to remove. The basis for this choice is
181 * the ranking information (below groups and prio denote the arrays
182 * ranking->groups and ranking->prio and t and k denote the knot vector
183 * and order of rankcurve). Suppose we are to remove r knots.
184 * First the smallest integer i satisying groups(i) >= r is determined.
185 * In prio all the interior knots of t are listed (or rather their
186 * index in t, t[k-1+5] has index 5), and since we are to remove
187 * all knots in the first ranking groups, we remove the first groups[i-1]
188 * knots of prio (assuming for simplicity that groups[-1]=0).
189 * In the last group we reach into (which consists of knots no.
190 * prio[groups[i-1]+1], prio[groups[i-1]+2],...,prio[groups[i]])
191 * knots are removed uniformly on index, meaning that if we are to remove
192 * half the knots in this group, we pick every other knot as they are
193 * listed in prio, if we are to remove 1/3 of the knots, we remove every
194 * three knots and so on. In general, it is difficult to say what
195 * uniformly on subscripts should mean (how do you choose one knot from
196 * six in a symmetric way) but the approach used here seems to work
197 * reasonably well (see the program text below).
198 * When it is known which knots to remove, the corresponding knot vector
199 * can be constructed and an approximation computed.
200 * the error is then checked as indicated above, and the iterations proceed.
201 *
202 * For more information cf. the references:
203 *
204 *
205 * References:
206 * 1. A Data-Reduction Strategy for Splines with Applications to the
207 * Approximation of Functions and data, IMA J. of Num. Anal. 8 (1988),
208 * pp. 185-208.
209 *
210 * 2. Knot Removal for Parametric B-spline Curves and Surfaces,
211 * CAGD 4 (1987), pp. 217-230.
212 *
213 *
214 * Calls: sh1365, s6err
215 *
216 * Written by : Knut Moerken, University of Oslo, July 1992, based on an
217 * earlier Fortran version.
218 * Changed by: Paal Fugelli, SINTEF, 1994-07.
219 * Changed to fix several major memory leakage problems.
220 *
221 *********************************************************************
222 */
223 {
224 int k = oldcurve->ik; /* Make some parameters more easily */
225 int dim = oldcurve->idim; /* accessible. */
226 int mprio = rankcurve->in;
227 int *prio = ranking->prio;
228 int *group = ranking->groups;
229 int antrem = ranking->antrem;
230 int antgr = ranking->antgr;
231 char *del_array; /* Boolean array that is used for
232 marking what knots to include. */
233 char big, bigco; /* Boolean variables that are used
234 for checking whether the error
235 is small enough. */
236 int lstat=0; /* Local status variable. */
237 int pos=0; /* Parameter to s6err. */
238 int nlim = MAX(k, startfix+endfix); /* The minimum number of B-spline
239 coefficients that should be
240 left in newcurve. */
241
242 /* For the use of the other variables,
243 see the code below. */
244 int i, start, stop, indx, count, r, p, hn;
245 SISLCurve *hcurve = SISL_NULL;
246 double h;
247 double *local_err = SISL_NULL, *l2_err = SISL_NULL, *ltau = SISL_NULL;
248
249 /* Allocate memory for the local arrays. */
250
251 del_array = newarray(mprio-k, char);
252 if (del_array == SISL_NULL) goto err101;
253
254 /* In case we do not enter the while loop at all we must give newcurve
255 a value. */
256
257 *newcurve = newCurve(mprio, k, rankcurve->et, rankcurve->ecoef, 1, dim, 1);
258 if (newcurve == SISL_NULL) goto err101;
259
260 /* Iterate by binary search until the lower and upper bound on how many
261 knots to include are essentially equal. */
262
263 while (mini+1 < maxi)
264 {
265
266 /* To start with none of the knots of rankspline are marked
267 for removal. */
268
269 for (i=0; i<mprio-k; i++) del_array[i] = 0;
270
271 /* We then have to find out which knots to remove. We remove knots
272 group by group, with start pointing (into prio) to the first knot
273 of the current group and stop to the one following the last of this
274 group. */
275
276 start = 0;
277 stop = group[0];
278 count = 0;
279
280 /* We remove all knots of each group and mark them in del_array
281 until this would mean that we have removed too many. */
282
283 while (stop <= antrem)
284 {
285 for (i=start; i<stop; i++) del_array[prio[i]] = 1;
286
287 count++;
288
289 if (count < antgr)
290 {
291 start = stop;
292 stop = group[count];
293 }
294 else
295 {
296
297 /* start=stop signifies that we have reached the end. */
298
299 stop = stop + 1;
300 start = stop + 1;
301 }
302 }
303
304 /* Now there are p more knots to remove from a group containing
305 p knots. These p knots are removed "uniformly on index". */
306
307 r = stop - start;
308 p = antrem - start;
309
310 if (p > 0)
311 {
312 h = (double) (r+1) / (double) p;
313 for (i=0; i<p; i++)
314 {
315 indx = start - 1 + (int) floor( h*(i+0.5)+0.5 );
316 del_array[prio[indx]] = 1;
317 }
318 }
319
320 /* Gives the number of coefficients in the new spline to be computed. */
321
322 hn = mprio - antrem ;
323
324 /* The new knot vector is stored in ltau. It might already be allocated
325 from the last iteration, so free it first if required. */
326
327 if (ltau != SISL_NULL) freearray(ltau);
328 ltau = newarray(hn+k, double);
329 if (ltau == SISL_NULL) goto err101;
330
331 /* Set the first and last k knots. */
332
333 for (i=0; i<k; i++)
334 {
335 ltau[i] = rankcurve->et[i];
336 ltau[i+hn] = rankcurve->et[i+mprio];
337 }
338
339 /* Set the remaining knots as indicated by del_array. */
340
341 for (indx=k, i=0; i<mprio-k; i++)
342 if (!del_array[i]) ltau[indx++] = rankcurve->et[i+k];
343
344 /* Compute an approximation on the new knot vector and store it
345 in hcurve.
346 sh1365() will allocate space for hcurve (with icopy==1), local_err
347 and l2_err.
348 Need local_err to store the error of the approximation until we know
349 that we have an approximation with error smaller than the tolerance,
350 then we can transfer the error to maxerr.
351 Must remember to free local_err and l2_err since there will be
352 a memory leak if they were allocated in the previous iteration. */
353
354 if (local_err != SISL_NULL) freearray(local_err);
355 if (l2_err != SISL_NULL) freearray(l2_err);
356
357 sh1365(oldcurve, ltau, k, hn, startfix, endfix,
358 &hcurve, &local_err, &l2_err, &lstat);
359 if (lstat < 0) goto err;
360
361 /* Check the error. big signifies that the error in one of the
362 components is larger than the tolerance. bigco signifies that
363 the error in one of the components is larger than the (usually
364 much smaller) tolerance epsco. */
365 big = 0;
366 bigco = 0;
367
368 for (i=0; i<dim; i++)
369 {
370 big = big || (local_err[i] > eps[i]);
371 bigco = bigco || (local_err[i] > epsco[i]);
372 }
373
374 /* The error is too large if it is big or if it is bigco and the
375 number of coefficients is smaller than nlim. The latter test is
376 important when the sum of the number of derivatives to be kept
377 fixed at the two ends is larger than k. */
378
379 big = big || (bigco && hn < nlim);
380
381 if (big)
382 {
383
384 /* If the error was too big, we just throw away hcurve and
385 indicate that an upper bound for the number of knots to
386 remove is antrem. */
387
388 if (hcurve != SISL_NULL) freeCurve(hcurve);
389 hcurve = SISL_NULL;
390 maxi = antrem;
391 }
392 else
393 {
394
395 /* If the error is acceptable we know that we can remove at least
396 antrem knots and store hcurve in newcurve. We also save the
397 error in maxerr. */
398
399 mini = antrem;
400 if (*newcurve != SISL_NULL) freeCurve(*newcurve);
401 *newcurve = hcurve;
402 hcurve = SISL_NULL;
403 for (i=0; i<dim; i++) maxerr[i] = local_err[i];
404 }
405
406 /* The number of knots to be removed next time is half way between
407 mini and maxi. */
408
409 antrem = mini + (maxi-mini)/2;
410 }
411
412 *stat = 0;
413
414 goto out;
415
416
417
418 /* Error in memory allocation. */
419
420 err101:
421 *stat = -101;
422 goto out;
423
424 /* Error in lower level routine. */
425
426 err:
427 *stat = lstat;
428 s6err("s1354", *stat, pos);
429
430 /* Clean up before exit. */
431
432 out:
433 if (hcurve != SISL_NULL) freeCurve(hcurve);
434 if (del_array != SISL_NULL) freearray(del_array);
435 if (local_err != SISL_NULL) freearray(local_err);
436 if (l2_err != SISL_NULL) freearray(l2_err);
437 if (ltau != SISL_NULL) freearray(ltau);
438
439 return;
440 }
441