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 #define S1961
43
44 #include "sislP.h"
45
46 #if defined(SISLNEEDPROTOTYPES)
s1961(double ep[],int im,int idim,int ipar,double epar[],double eeps[],int ilend,int irend,int iopen,double afctol,int itmax,int ik,SISLCurve ** rc,double emxerr[],int * jstat)47 void s1961(double ep[],int im,int idim,int ipar,double epar[],
48 double eeps[],int ilend,int irend,int iopen,double afctol,
49 int itmax,int ik,SISLCurve **rc,double emxerr[],
50 int *jstat)
51 #else
52 void s1961(ep,im,idim,ipar,epar,eeps,ilend,irend,iopen,afctol,
53 itmax,ik,rc,emxerr,jstat)
54 double ep[];
55 int im;
56 int idim;
57 int ipar;
58 double epar[];
59 double eeps[];
60 int ilend;
61 int irend;
62 int iopen;
63 double afctol;
64 int itmax;
65 int ik;
66 SISLCurve **rc;
67 double emxerr[];
68 int *jstat;
69 #endif
70 /*
71 ************************************************************
72 *
73 * Purpose: To compute a spline-approximation to the data given by the
74 * points ep, and represent it as a B-spline curve with
75 * parameterization determined by the parameter ipar.
76 * The approximation is determined by first forming the piecewise
77 * linear interpolant to the data, and then performing knot
78 * removal on this initial approximation.
79 *
80 * Input :
81 * Ep - Array (length idim*im) containing the points to
82 * be approximated.
83 * Im - The no. of data points.
84 * Idim - The dimension of the euclidean space in which the data
85 * points lie, i.e. the number of components of each data point.
86 * Ipar - Flag indicating the type of parameterization to be used:
87 * = 1 : Paramterize by accumulated cord length.
88 * (Arc length parametrization for the piecewise
89 * linear interpolant.)
90 * = 2 : Uniform parameterization.
91 * = 3 : Parametrization given by epar.
92 * If ipar<1 or ipar>3, it will be set to 1.
93 * Epar - Array (length im) containing a parametrization
94 * of the given data.
95 * Eeps - Array (length idim) containing the tolerance to be
96 * used during the data reduction stage. The final
97 * approximation to the data will deviate less than eeps
98 * from the piecewise linear interpolant in each of the
99 * idim components.
100 * Ilend - The no. of derivatives that are not allowed to change
101 * at the left end of the curve.
102 * The (0 - (ilend-1)) derivatives will be kept fixed.
103 * If ilend <0, this routine will set it to 0.
104 * If ilend<ik, this routine will set it to ik.
105 * Irend - The no. of derivatives that are not allowed to change
106 * at the right end of the curve.
107 * The (0 - (irend-1)) derivatives will be kept fixed.
108 * If irend <0, this routine will set it to 0.
109 * If irend<ik, this routine will set it to ik.
110 * iopen - Open/closed parameter.
111 * = 1 : Produce open curve.
112 * = 0 : Produce closed, non-periodic curve.
113 * = -1 : Produce closed, periodic curve.
114 * If a closed or periodic curve is to be produced and the
115 * start- and endpoint is more distant than the length of
116 * the tolerance, a new point is added. Note that if the
117 * parametrization is given as input, the parametrization
118 * if the last point will be arbitrary.
119 * Afctol = Number indicating how the tolerance is to be shared
120 * between the two data reduction stages. For the linear
121 * reduction, a tolerance of afctol*eeps will be used,
122 * while a tolerance of (1-afctol)*eeps will be used
123 * during the final data reduction. (Similarly for edgeps.)
124 * Itmax - Max. no. of iterations in the data-reduction routine.
125 * Ik - The polynomial order of the approximation.
126 *
127 * Output:
128 * Jstat - Output staus:
129 * < 0 : Error.
130 * = 0 : Ok.
131 * > o : Warning.
132 * Rc - Pointer to curve.
133 * Emxerr - Array (length idim) (allocated outside this routine.)
134 * containing for each component an upper bound on the
135 * max. deviation of the final approximation from the
136 * initial piecewise linear interpolant.
137 *
138 * Method: First the piecewise linear interpoolant is computed by a call
139 * to s1350 or s1351, and then knots are removed from this
140 * interpolant. The degree of the resulting linear spline is then
141 * raised to ik, and a second data-reduction is performed.
142 *
143 * Written by: C.R.Birkeland Si Oslo,Norway April 1993.
144 * Changed and renamed by : Vibeke Skytt, SINTEF Oslo, 02.95. Introduced
145 * periodicity.
146 ********************************************************************
147 */
148 {
149 double *maxerr = SISL_NULL; /* Arrays used to store error estimates */
150 double *error1 = SISL_NULL;
151 int i;
152 int stat = 0; /* Loop control variables */
153 int kpos = 0;
154 SISLCurve *ocurve = SISL_NULL; /* Local spline curve */
155 double *sp = SISL_NULL; /* Extended data points in closed/periodic case. */
156 double *spar = SISL_NULL; /* Extended par. values. */
157
158 /* Check Input */
159
160 if (im < 2 || idim < 1) goto err103;
161 if (ipar < 1 || ipar > 3) ipar = 1;
162
163 if ((iopen == 0 || iopen == -1) &&
164 s6dist(ep, ep+(im-1)*idim, idim) > s6length(eeps, idim, &stat))
165 {
166 /* Add an extra point to the input points. First allocated scratch
167 for extended arrays. */
168
169 if ((sp = newarray((im+1)*idim, DOUBLE)) == SISL_NULL) goto err101;
170 memcopy(sp, ep, im*idim, DOUBLE);
171 memcopy(sp+im*idim, ep, idim, DOUBLE);
172 if (ipar == 3)
173 {
174 if ((spar = newarray(im+1, DOUBLE)) == SISL_NULL) goto err101;
175 memcopy(spar, epar, im, DOUBLE);
176 spar[im] = spar[im-1] + s6dist(sp+(im-1)*idim, sp+im*idim, idim);
177 }
178 }
179 else
180 {
181 sp = ep;
182 spar = epar;
183 }
184
185 /* Set default value for afctol */
186
187 if(afctol < 0 || afctol > 1) afctol =0;
188
189 /* Piecewise linear interpolant to the data */
190
191 if (ipar == 3)
192 s1350(sp, spar, im, idim, 2, &ocurve, &stat);
193 else
194 s1351(sp, ipar, im, idim, 2, &ocurve, &stat);
195 if (stat<0) goto error;
196
197 /* Compute tolerance for linear reduction */
198
199 error1 = newarray(idim, DOUBLE);
200 maxerr = newarray(idim, DOUBLE);
201 if (error1 == SISL_NULL || maxerr == SISL_NULL) goto err101;
202 for (i=0; i<idim; i++)
203 maxerr[i] = eeps[i]*afctol;
204
205 /* Data reduction on piecewise linear interpolant */
206
207 s1940(ocurve, maxerr, ilend, irend, iopen, itmax, rc,
208 error1, &stat);
209 if (stat<0) goto error;
210 freeCurve(ocurve);
211
212 /* Piecewise linear interpolant to the reduced data
213 (coefficients of linear spline curve rc)
214 expressed as B-spline of order ik */
215
216 s1350((*rc)->ecoef, &((*rc)->et)[1], (*rc)->in,
217 idim, ik, &ocurve, &stat);
218 if (stat<0) goto error;
219 freeCurve(*rc);
220
221 /* Compute tolerance for last reduction */
222
223 for (i=0; i<idim; i++)
224 maxerr[i] = eeps[i] - error1[i];
225
226 /* Perform data reduction on the hermite interpolant */
227
228 s1940(ocurve, maxerr, ilend, irend, iopen, itmax, rc,
229 emxerr, &stat);
230 if (stat<0) goto error;
231
232 /* Compute total error after reduction */
233
234 for (i=0; i<idim; i++)
235 emxerr[i] += error1[i];
236
237 *jstat = 0;
238 goto out;
239
240 /* Error in scratch allocation. */
241
242 err101 :
243 *jstat = -101;
244 s6err("s1961",*jstat,kpos);
245 goto out;
246
247 /* Error in input */
248
249 err103:
250 *jstat = -103;
251 s6err("s1961",*jstat,kpos);
252 goto out;
253
254 /* Error in lower level routine. */
255
256 error:
257 *jstat = stat;
258 s6err("s1961",*jstat,kpos);
259 goto out;
260
261 /* Exit */
262
263 out:
264 if (maxerr != SISL_NULL) freearray(maxerr);
265 if (error1 != SISL_NULL) freearray(error1);
266 if (ocurve != SISL_NULL) freeCurve(ocurve);
267 if (spar != SISL_NULL && spar != epar) freearray(spar);
268 if (sp != SISL_NULL && sp != ep) freearray(sp);
269 return;
270 }
271