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: s1919.c,v 1.2 2001-03-19 15:58:56 afr Exp $
45 *
46 */
47
48
49 #define S1919
50
51 #include "sislP.h"
52
53 #if defined(SISLNEEDPROTOTYPES)
54 void
s1919(double et[],double prev[],double curr[],double deriv[],double follow[],int in,int ik,int idim,int iip,int iif,double ap,double ac,double af,int * jstat)55 s1919 (double et[], double prev[], double curr[], double deriv[],
56 double follow[], int in, int ik, int idim, int iip, int iif,
57 double ap, double ac, double af, int *jstat)
58 #else
59 void
60 s1919 (et, prev, curr, deriv, follow, in, ik, idim, iip, iif, ap, ac, af, jstat)
61 double et[];
62 double prev[];
63 double curr[];
64 double deriv[];
65 double follow[];
66 int in;
67 int ik;
68 int idim;
69 int iip;
70 int iif;
71 double ap;
72 double ac;
73 double af;
74 int *jstat;
75 #endif
76 /*
77 *********************************************************************
78 *
79 *********************************************************************
80 *
81 * PURPOSE : To make the product of the tangent legth and the parametrization
82 * of the curve epd as close as possible to the curve length.
83 *
84 * INPUT : et - The original knot vector.
85 *
86 * OUTPUT :jstat - Status variable:
87 * > 0 : warning
88 * = 0 : ok
89 * < 0 : error
90 *
91 * METHOD :
92 *
93 * REFERENCES : Fortran version by Tor Dokken, SI, 1991-07
94 *
95 * CALLS : s1890, s1221, s1891, s6err.
96 *
97 * WRITTEN BY : Trond Vidar Stensby, SI, 1991-07
98 *
99 *********************************************************************
100 */
101 {
102 SISLCurve *tcurve = SISL_NULL; /* Temporary curve. */
103 int knh; /* Local number of verticews. */
104 int left; /* Used when calling s1221 */
105 int pos, pos2; /* Used for efficent adressing of arrays. */
106 int ki, kj; /* Loop control variables. */
107 double tdist1; /* Distance between previous and current. */
108 double tdist2; /* Distance between current and following. */
109 double tlength; /* Length of tangemt. */
110 double tdiff; /* Dummy. */
111 double tfak;
112 double tval1; /* Adjusted parameter values. */
113 double tval2;
114 double *kpar = SISL_NULL; /* Used when calculating parametrization. */
115 int *kder = SISL_NULL;
116 double *kpc = SISL_NULL; /* Points on previous curve. */
117 double *kdc = SISL_NULL; /* Points on derivative curve. */
118 double *kcc = SISL_NULL; /* Points on current curve. */
119 double *kfc = SISL_NULL; /* Points on following curve. */
120 double *epd = SISL_NULL; /* Vertices of the derivative curve. */
121 int kstat = 0;
122 int kpos = 0;
123
124 *jstat = 0;
125
126
127 /* Test if legal input. */
128
129 if (ik <= 1 || in <ik)
130 goto err112;
131
132
133 /* Produce parameter values. */
134
135 s1890 (et, ik, in, &kpar, &kder, &kstat);
136 if (kstat < 0)
137 goto error;
138
139
140 /* Allocate temporary arrays. */
141
142 kpc = newarray (idim * in, DOUBLE);
143 if (kpc == SISL_NULL)
144 goto err101;
145 kcc = newarray (idim * in, DOUBLE);
146 if (kcc == SISL_NULL)
147 goto err101;
148 kdc = newarray (idim * in, DOUBLE);
149 if (kdc == SISL_NULL)
150 goto err101;
151 kfc = newarray (idim * in, DOUBLE);
152 if (kfc == SISL_NULL)
153 goto err101;
154
155
156 if (iip == 1)
157 {
158 /* Caculate interpolation points on previous curve. */
159
160 tcurve = newCurve (in, ik, et, prev, 1, idim, 1);
161 if (tcurve == SISL_NULL)
162 goto err101;
163
164 left = 0;
165 for (ki = 0; ki < in; ki++)
166 {
167 s1221 (tcurve, 0, kpar[ki], &left, &kpc[ki * idim], &kstat);
168 if (kstat < 0)
169 goto error;
170 }
171 if (tcurve != SISL_NULL)
172 freeCurve (tcurve);
173 }
174
175 /* Caculate interpolation points on current curve. */
176
177 tcurve = newCurve (in, ik, et, curr, 1, idim, 1);
178 if (tcurve == SISL_NULL)
179 goto err101;
180
181 left = 0;
182 for (ki = 0; ki < in; ki++)
183 {
184 s1221 (tcurve, 0, kpar[ki], &left, &kcc[ki * idim], &kstat);
185 if (kstat < 0)
186 goto error;
187 }
188 if (tcurve != SISL_NULL)
189 freeCurve (tcurve);
190
191
192 /* Caculate interpolation points on derivative curve. */
193
194 tcurve = newCurve (in, ik, et, deriv, 1, idim, 1);
195 if (tcurve == SISL_NULL)
196 goto err101;
197
198 left = 0;
199 for (ki = 0; ki < in; ki++)
200 {
201 s1221 (tcurve, 0, kpar[ki], &left, &kdc[ki * idim], &kstat);
202 if (kstat < 0)
203 goto error;
204 }
205 if (tcurve != SISL_NULL)
206 freeCurve (tcurve);
207
208
209 if (iif == 1)
210 {
211 /* Caculate interpolation points on following curve. */
212
213 tcurve = newCurve (in, ik, et, follow, 1, idim, 1);
214 if (tcurve == SISL_NULL)
215 goto err101;
216
217 left = 0;
218 for (ki = 0; ki < in; ki++)
219 {
220 s1221 (tcurve, 0, kpar[ki], &left, &kfc[ki * idim], &kstat);
221 if (kstat < 0)
222 goto error;
223 }
224 if (tcurve != SISL_NULL)
225 freeCurve (tcurve);
226 }
227
228 /* Adjust the points calculated on the derivative curve according to the
229 distances between previous, current and following curve. */
230
231 if (iip == 1)
232 tval1 = fabs (ac - ap);
233 if (iif == 1)
234 tval2 = fabs (af - ac);
235
236 pos = 0;
237 for (ki = 0; ki < in; ki++)
238 {
239 tdist1 = (double) 0.0;
240 tdist2 = (double) 0.0;
241 tlength = (double) 0.0;
242
243 for (kj = 0; kj < idim; kj++)
244 {
245 if (iip == 1)
246 {
247 tdiff = kcc[pos] - kpc[pos];
248 tdist1 += tdiff * tdiff;
249 }
250 if (iif == 1)
251 {
252 tdiff = kfc[pos] - kcc[pos];
253 tdist2 += tdiff * tdiff;
254 }
255 tlength += kdc[pos] * kdc[pos];
256 pos++;
257 }
258
259 tdist1 = sqrt (tdist1);
260 tdist2 = sqrt (tdist2);
261 tlength = sqrt (tlength);
262 if (tlength == (double) 0.0)
263 tlength = (double) 1.0;
264
265
266 /* Make scaling factor of tangent/derivative. */
267
268 tfak = (double) 1.0;
269 if (iip == 1 && iif != 1)
270 tfak = tdist1 / (tval1 * tlength);
271 else if (iip != 1 && iif == 1)
272 tfak = tdist2 / (tval2 * tlength);
273 else if (iip == 1 && iif == 1)
274 tfak = min (tdist1 / (tval1 * tlength),
275 tdist2 / (tval2 * tlength));
276
277
278 /* Find actual length of derivative curve at the parameter value.
279 Make new derivative length. */
280
281 pos2 = pos - idim;
282 for (kj = 0; kj < idim; kj++)
283 {
284 kdc[pos2 + kj] *= tfak;
285 }
286 }
287
288 /* Calculate new curve description. */
289
290 s1891 (kpar, kdc, idim, in, 1, kder, 1, et, &epd, &knh, ik, 0, 0, &kstat);
291 if (kstat < 0)
292 goto error;
293
294 memcopy (deriv, epd, idim * in, DOUBLE);
295
296
297 /* OK */
298
299 goto out;
300
301
302 /* Allocation error. */
303
304 err101:
305 *jstat = -101;
306 s6err ("s1919", *jstat, kpos);
307 goto out;
308
309 /* Error in description of B-spline. */
310
311 err112:
312 *jstat = -112;
313 s6err ("s1919", *jstat, kpos);
314 goto out;
315
316 /* Error in lower level routine. */
317
318 error:
319 *jstat = kstat;
320 s6err ("s1919", *jstat, kpos);
321 goto out;
322
323 out:
324 if (epd != SISL_NULL)
325 freearray (epd);
326 if (kpc != SISL_NULL)
327 freearray (kpc);
328 if (kcc != SISL_NULL)
329 freearray (kcc);
330 if (kdc != SISL_NULL)
331 freearray (kdc);
332 if (kfc != SISL_NULL)
333 freearray (kfc);
334 if (kpar != SISL_NULL)
335 freearray (kpar);
336 if (kder != SISL_NULL)
337 freearray (kder);
338
339 return;
340 }
341