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: s1358.c,v 1.2 2001-03-19 15:58:47 afr Exp $
45 *
46 */
47
48
49 #define S1358
50
51 #include "sislP.h"
52
53 #if defined(SISLNEEDPROTOTYPES)
54 void
s1358(double epoint[],int inbpnt,int idim,double ntype[],double epar[],int icnsta,int icnend,int iopen,int ik,double astpar,double * cendpar,SISLCurve ** rc,double ** gpar,int * jnbpar,int * jstat)55 s1358(double epoint[],int inbpnt,int idim,double ntype[],double epar[],
56 int icnsta,int icnend,int iopen,int ik,double astpar,
57 double *cendpar,SISLCurve **rc,double **gpar,int *jnbpar,int *jstat)
58 #else
59 void s1358(epoint,inbpnt,idim,ntype,epar,icnsta,icnend,iopen,ik,astpar,
60 cendpar,rc,gpar,jnbpar,jstat)
61 double epoint[];
62 int inbpnt;
63 int idim;
64 double ntype[];
65 double epar[];
66 int icnsta;
67 int icnend;
68 int iopen;
69 int ik;
70 double astpar;
71 double *cendpar;
72 SISLCurve **rc;
73 double **gpar;
74 int *jnbpar;
75 int *jstat;
76 #endif
77 /*
78 *************************************************************************
79 *
80 * Purpose: To calculate a B-spline curve interpolating a set of points.
81 * The points can be assigned a tangent (derivative).
82 * The curve can be closed or open. If end-conditions are
83 * conflicting, the condition closed curve rules out other
84 * end conditions.
85 * The parametrization is given by the array Epar.
86 *
87 * Input:
88 * Epoint - Array (length idim*inbpnt) containing the points/
89 * derivatives to be interpolated.
90 * Inbpnt - No. of points/derivatives in the epoint array.
91 * Idim - The dimension of the space in which the points lie.
92 * ntype - Array (length inbpnt) containing type indicator for
93 * points/derivatives/second-derivatives:
94 * 1 - Ordinary point.
95 * 2 - Knuckle point. (Is treated as an ordinary point.)
96 * 3 - Derivative to next point.
97 * 4 - Derivative to prior point.
98 * ( 5 - Second derivative to next point. )
99 * ( 6 - Second derivative to prior point. )
100 * 13 - Start-point of tangent to next point.
101 * 14 - End-point of tangent to prior point.
102 * Epar - Array containing the wanted parametrization. Only parameter
103 * values corresponding to position points are given.
104 * For closed curves, one additional parameter value
105 * must be spesified. The last entry contains
106 * the parametrization of the repeted start point.
107 * (if the endpoint is equal to the startpoint of
108 * the interpolation the lenght of the array should
109 * be equal to inpt1 also in the closed case).
110 * Icnsta - Additional condition at the start of the curve:
111 * 0 : No additional condition.
112 * 1 : Zero curvature at start.
113 * Icnend - Additional condition at the end of the curve:
114 * 0 : No additional condition.
115 * 1 : Zero curvature at end.
116 * Iopen - Flag telling if the curve should be open or closed:
117 * Ik - The order of the B-spline curve to be produced.
118 * Astpar - Parameter-value to be used at the start of the curve.
119 *
120 * Output:
121 * Jstat - status variable:
122 * < 0 : Error.
123 * = 0 : Ok.
124 * > 0 : Warning.
125 * cendpar - Parameter-value used at the end of the curve.
126 * Rc - Pointer to output-curve.
127 * Gpar - Pointer to the parameter-values of the points in
128 * the curve. Represented only once, although derivatives
129 * and second-derivatives will have the same parameter-
130 * value as the points.
131 * Jnbpar - No. of different parameter-values.
132 *
133 * Method: First the parametrization of the curve and the type
134 * specification of points/derivatives are checked and/or
135 * corrected. Then the knots are calculated, and the
136 * interpolation is performed.
137 *
138 * Calls: s1907,s1908,s1902,s1891,s1713,s1750,s6err.
139 *
140 * Written by: Trond Vidar Stensby, SI, 1991-07
141 * The Fortran routine, p19539, is written by: T. Dokken SI.
142 * Bug fix: Michael Floater, 82.9.93. The construction of gpar and *jnbpar
143 * at the end of the routine was incorrect.
144 *****************************************************************
145 */
146 {
147 int kpos = 0;
148 SISLCurve *qc = SISL_NULL; /* Temporary SISLCurves. */
149 SISLCurve *qc2 = SISL_NULL;
150 int *ltype = SISL_NULL; /* Type of interpolation condition
151 (temporary) */
152 int *ltype2 = SISL_NULL; /* Type of interpolation condition (finial) */
153 int *sder = SISL_NULL; /* Vector of derivative indicators. */
154 int knpt; /* Number of interpolation conditions. */
155 int kordr; /* Local order. */
156 int kstat; /* Status variable. */
157 int kn; /* Number of coefficients of B-spline curve.*/
158 int ki; /* Loop control variable. */
159 int kright = 1; /* One equation system to solve in
160 interpolation. */
161 int knlr = 0; /* Indicates shape of interpolation matrix. */
162 int knrc = 0; /* Indicates shape of interpolation matrix. */
163 double *lpar = SISL_NULL; /* Parameter values. (temporary) */
164 double *lcond = SISL_NULL; /* Interpolation conditions. (temporary) */
165 double *sknot = SISL_NULL; /* Knot vector. */
166 double *spar = SISL_NULL; /* Parameter valued. (finial) */
167 double *scond = SISL_NULL; /* Interpolation conditions. (finial) */
168 double *scoef = SISL_NULL; /* Coefficients of curve. */
169 int *ityp = SISL_NULL;
170
171 /* -> new statement guen & ujk Thu Jul 2 14:59:05 MESZ 1992 */
172
173 /* make compatible to old use of s1604 */
174
175 if (iopen==SISL_CRV_CLOSED) iopen = SISL_CRV_PERIODIC;
176
177 /* a new version with input-iopen == rc->cuopen
178 should be made, with a name different from any
179 other.
180 NOTE: There is an error in this function when iopen = 0
181 qc as input to s1713 (and s1750) then has wrong flag !!*/
182
183 /* <- new statement guen & ujk Thu Jul 2 14:59:05 MESZ 1992 */
184
185 ityp = newarray (inbpnt, INT);
186 if (ityp == SISL_NULL) goto err101;
187
188 for (ki=0; ki < inbpnt; ki++) ityp[ki] = (int)ntype[ki];
189 *jstat = 0;
190
191 /* Transform interpolation conditions. */
192
193 s1907 (epoint, ityp, epar, iopen, icnsta, icnend, inbpnt,
194 idim, &lcond, <ype, &lpar,&knpt, &kstat);
195 if (kstat < 0) goto error;
196
197 /* Test interpolation conditions, and adjust the input conditions
198 if necessary. */
199
200 s1908 (lcond, ltype, lpar, knpt, ik, idim, iopen, &scond, <ype2,
201 &spar, &knpt, &kstat);
202 if (kstat < 0) goto error;
203
204 /* Allocate scratch for derivative indicator. */
205
206 sder = newarray (knpt, INT);
207 if (sder == SISL_NULL) goto err101;
208
209 for (ki = 0; ki < knpt; ki++)
210 sder[ki] = abs (ltype2[ki]);
211
212 kordr = MIN (ik, knpt);
213
214 if (!(iopen == SISL_CRV_OPEN))
215 {
216 knlr = kordr / 2;
217 knrc = kordr - knlr - 1;
218 knpt--;
219 }
220
221 /* Produce knot vector. */
222
223 s1902 (lpar, knpt, kordr, iopen, &sknot, &kstat);
224 if (kstat < 0) goto error;
225
226 /* Perform interpolation. */
227
228 s1891 (spar, scond, idim, knpt, kright, sder, iopen, sknot,
229 &scoef, &kn, kordr, knlr, knrc, &kstat);
230 if (kstat < 0) goto error;
231
232 /* Express the curve as a curve object. */
233
234 qc = newCurve (kn, kordr, sknot, scoef, 1, idim, 1);
235 if (qc == SISL_NULL) goto err101;
236 qc->cuopen = (iopen == SISL_CRV_OPEN) ? iopen : SISL_CRV_PERIODIC;
237
238 if (iopen == SISL_CRV_CLOSED)
239 {
240 /* A closed, non-periodic curve is expected. Pick the part of the
241 interpolation curve that has got a full basis. */
242
243 s1713 (qc, sknot[kordr - 1], sknot[kn], &qc2, &kstat);
244 if (kstat < 0) goto error;
245
246 if (qc != SISL_NULL) freeCurve (qc);
247 qc = qc2;
248 }
249
250 if (kordr < ik)
251 {
252 /* The order of the curve is less than expected. Increase the order. */
253
254 qc2 = SISL_NULL;
255 s1750 (qc, ik, &qc2, &kstat);
256 if (kstat < 0) goto error;
257
258 if (qc != SISL_NULL) freeCurve (qc);
259 qc = qc2;
260 }
261
262 /* Set open/closed parameter of curve. */
263
264 qc->cuopen = iopen;
265
266 /* Set end of parameter interval. */
267
268 *cendpar = *(qc->et + qc->in);
269
270 /* Interpolation performed. */
271
272 /* Find distinct parameter values. */
273
274 *gpar = lpar;
275
276 /* Bug fix, MSF, 28.9.93. Change
277 *jnbpar = 1;
278 to: */
279
280 *jnbpar = 0;
281
282 /*called from s1333: s1358(epoint,21,21,only 1,epar,0,0,1,3,,
283 cendpar,rc,gpar,jnbpar,jstat) */
284 /* Error: iopen = 1, order=3, hp has a case where the size of gpar is 21 and
285 when leaving this loop *jnbpar (and knpt and ki) is also 21 !!!!!!!!!!!!? */
286 /* The above error is now fixed, MSF.
287 The following chunk of code merely removes double knots
288 from the knot vector spar (length knpt) producing
289 both the number *jnbpar of different knots and
290 an array gpar with the different knots of length *jnbpar. */
291
292 for (ki = 1; ki<knpt; ki++)
293 {
294 if (spar[ki - 1] < spar[ki])
295
296 /* Bug fix. MSF, 28.9.93 Change (*gpar)[(*jnbpar)++] = spar[ki]; to: */
297
298 (*gpar)[(*jnbpar)++] = spar[ki-1];
299 }
300
301 /* Bug fix. MSF, 28.9.93 Change (*gpar)[(*jnbpar)++] = spar[ki]; to: */
302
303 (*gpar)[(*jnbpar)++] = spar[ki-1];
304
305 *rc = qc;
306 goto out;
307
308
309 /* Error in scratch allocation. */
310
311 err101:
312 *jstat = -101;
313 s6err ("s1358", *jstat, kpos);
314 goto out;
315
316 /* Error in lower level routine. */
317
318 error:
319 *jstat = kstat;
320 s6err ("s1358", *jstat, kpos);
321 goto out;
322
323 out:
324 /* Free scratch occupied by local arrays. */
325
326 if (scond != SISL_NULL) freearray (scond);
327 if (scoef != SISL_NULL) freearray (scoef);
328 if (sder != SISL_NULL) freearray (sder);
329 if (ltype != SISL_NULL) freearray (ltype);
330
331 /* -> added, guen & ujk Wed Jul 1 18:48:27 MESZ 1992 */
332 if (ityp != SISL_NULL) freearray (ityp);
333 /* <- added, guen & ujk Wed Jul 1 18:48:27 MESZ 1992 */
334
335 if (ltype2 != SISL_NULL) freearray (ltype2);
336 if (lcond != SISL_NULL) freearray (lcond);
337 if (sknot != SISL_NULL) freearray (sknot);
338 if (spar != SISL_NULL) freearray (spar);
339
340 return;
341 }
342
343