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: s1227.c,v 1.2 2001-03-19 15:58:42 afr Exp $
45 *
46 */
47
48
49 #define S1227
50
51 #include "sislP.h"
52
53 #if defined(SISLNEEDPROTOTYPES)
54 void
s1227(SISLCurve * pc1,int ider,double ax,int * ileft,double eder[],int * jstat)55 s1227(SISLCurve *pc1,int ider,double ax,int *ileft,double eder[],int *jstat)
56 #else
57 void s1227(pc1,ider,ax,ileft,eder,jstat)
58 SISLCurve *pc1;
59 int ider;
60 double ax;
61 int *ileft;
62 double eder[];
63 int *jstat;
64 #endif
65 /*
66 *********************************************************************
67 *
68 *********************************************************************
69 *
70 * PURPOSE : To compute the value and ider first derivatives of the
71 * B-spline curve pointed to by pc1, at the point with
72 * parameter value ax from the left hand side.
73 *
74 *
75 *
76 * INPUT : pc1 - Pointer to the curve for which position
77 * and derivatives are to be computed.
78 * ider - The number of derivatives to compute.
79 * < 0 : Error.
80 * = 0 : Compute position.
81 * = 1 : Compute position and first derivative.
82 * etc.
83 * ax - The parameter value at which to compute
84 * position and derivatives.
85 *
86 *
87 *
88 * INPUT/OUTPUT : ileft - Pointer to the interval in the knot vector
89 * where ax is located. If et is the knot vector,
90 * the relation
91 *
92 * et[ileft] < ax <= et[ileft+1]
93 *
94 * should hold. (If ax == et[ik-1] then ileft should
95 * be ik-1. Here in is the number of B-spline
96 * coefficients.)
97 * If ileft does not have the right value upon
98 * entry to the routine, its value will be changed
99 * to the value satisfying the above condition.
100 *
101 *
102 *
103 * OUTPUT : eder - Double array of dimension [(ider+1)*idim]
104 * containing the position and derivative vectors.
105 * (idim is the number of components of each B-spline
106 * coefficient, i.e. the dimension of the Euclidean
107 * space in which the curve lies.)
108 * These vectors are stored in the following order:
109 * First the idim components of the position vector,
110 * then the idim components of the tangent vector,
111 * then the idim components of the second derivative
112 * vector, and so on.
113 * (The C declaration of eder as a two dimensional array
114 * would therefore be eder[ider+1,idim].)
115 * jstat - Status messages
116 * > 0 : Warning.
117 * = 0 : Ok.
118 * < 0 : Error.
119 *
120 *
121 * METHOD : First we find if the parameter ax is at a knot value, if
122 * so we artificially shortens the curve to end at this value.
123 * Then the traditional calculation from the right hand side
124 * can be used.
125 *
126 * Suppose that the given curve is of the form
127 *
128 * f(x) = sum(i) c(i) B(i,k)(x)
129 *
130 * where c are the B-spline coefficients and B(i,k)(x) the
131 * B-splines accociated with the knot vector et.
132 * (For idim > 1 this is a vector equation with idim
133 * components; however, the B-spline is not a vector.)
134 * It is known that for a given value of x there are at most
135 * k (the order of the splines) nonzero B-splines.
136 * If ileft has the correct value these B-splines will be
137 *
138 * B(ileft-k+1,k),B(ileft+k+2,k),...,B(ileft,k).
139 *
140 * The position and derivatives are computed by
141 * first computing the values and derivatives of all the
142 * B-splines at x, and then multiplying and summing.
143 *
144 * REFERENCES :
145 *
146 *-
147 * CALLS : compbder - Computes B-spline values and derivatives at
148 * a given point.
149 *
150 * WRITTEN BY : Knut Moerken, University of Oslo, August 1988.
151 * MODIFIED BY : Mike Floater, SI, Oslo, Norway, April 1991 for rational case.
152 * REVISED BY : Johannes Kaasa, SI, Aug. 92 (In case of NURBS the maximum
153 * derivative is not set equal order).
154 *
155 *********************************************************************
156 */
157 {
158 int kind; /* Type of curve */
159 int kstat=0; /* Local status variable. */
160 int kpos=0; /* The position of the error. */
161 int kn; /* The number of B-splines, i.e., the dimension of
162 the spline space associated with the knot
163 vector. */
164 int kk; /* The polynomial order of the curve. */
165 int kmult; /* Multiplicity of knot value */
166 int kdim; /* The dimension of the space in which the curve
167 lies. Equivalently, the number of components
168 of each B-spline coefficient. */
169 int kleft=0; /* Local version of ileft which is used in order to
170 avoid the pointer. */
171 int kder; /* Local version of ider. Since derivatives of order
172 higher than kk-1 are all zero, we set
173 kder = min(kk-1,ider). */
174 int ki,kj,kih,kjh; /* Control variables in for loops and for stepping
175 through arrays. */
176 int kl,kl1,kl2; /* Control variables in for loops and for stepping
177 through arrays. */
178 double *st; /* Pointer to the first element of the knot vector
179 of the curve. The knot vector has [kn+kk]
180 elements. */
181 double *scoef; /* Pointer to the first element of the curve's
182 B-spline coefficients. This is assumed to be an
183 array with [kn*kdim] elements stored in the
184 following order:
185 First the kdim components of the first B-spline
186 coefficient, then the kdim components of the
187 second B-spline coefficient and so on. */
188 double tt; /* Dummy variable used for holding an array element
189 in a for loop. */
190 double *ebder=SISL_NULL; /* Pointer to an array of dimension [kk*(ider+1)]
191 which will contain the values and ider first derivatives
192 of the kk nonzero B-splines at ax.
193 These are stored in the following order:
194 First the value, 1. derivative etc. of the
195 first nonzero B-spline, then the same for the
196 second nonzero B-spline and so on. */
197 double *sder=SISL_NULL; /* Pointer to array used for storage of points, if
198 non rational sder points to eder, if rational sder
199 has to be allocated to make room for the homogenous
200 coordinate */
201
202 /* Copy curve attributes to local parameters. */
203
204 kn = pc1 -> in;
205 kk = pc1 -> ik;
206 st = pc1 -> et;
207 scoef = pc1 -> ecoef;
208 kdim = pc1 -> idim;
209 kind = pc1 ->ikind;
210
211 if (kind == 2 || kind == 4)
212 {
213 scoef = pc1 -> rcoef;
214 kdim +=1;
215 sder = newarray(kdim*(ider+1),DOUBLE);
216 if (sder==SISL_NULL) goto err101;
217 }
218 else
219 {
220 scoef = pc1 -> ecoef;
221 sder = eder;
222 }
223
224 /* Check the input. */
225
226 if (kdim < 1) goto err102;
227
228 if (kk < 1) goto err110;
229
230 if (kn < kk) goto err111;
231
232 /* Find in which interval ax lies */
233
234 s1219(st,kk,kn,&kleft,ax,&kstat);
235 if (kstat<0) goto error;
236
237 /* To force the derivative to be taken from the left we artificially
238 shorten the curve if ax==st[kleft] */
239
240 kmult = s6knotmult(st,kk,kn,&kleft,ax,&kstat);
241 if (kstat < 0) goto error;
242
243
244 if (ax == st[kleft] && kleft > kk-1) kn = kleft-kmult+1;
245 if (st[kk-1] == st[kk] || st[kn-1] == st[kn]) goto err112;
246
247 if (ider < 0) goto err178;
248
249 if (pc1->ikind == 1 || pc1->ikind == 3)
250 kder = min(kk-1,ider);
251 else
252 kder = ider;
253
254 /* Allocate space for B-spline values and derivatives. */
255
256 ebder = newarray(kk*(kder+1),DOUBLE);
257 if (ebder == SISL_NULL) goto err101;
258
259 /* Set all the elements of sder to 0. */
260
261 for (ki=0; ki<(ider+1)*kdim; ki++) sder[ki] = (double)0.0;
262
263 /* Compute the values and derivatives of the nonzero B-splines and
264 update ileft if necessary. */
265
266 s1220(st,kk,kn,ileft,ax,kder,ebder,&kstat);
267
268 if (kstat < 0) goto error;
269
270 kleft = *ileft;
271
272 /* Multiply together as indicated above. */
273
274 /* ki steps through the appropriate kk B-spline coefficients while kih steps
275 through the B-spline value and derivatives for the B-spline given by ki.*/
276
277 kih = 0;
278 for (ki=kleft-kk+1; ki<=kleft; ki++)
279 {
280
281 /* kj counts through the kder+1 derivatives to be computed.
282 kjh steps through sder once for each ki to accumulate the contribution
283 from the different B-splines.
284 kl1 points to the first component of B-spline coefficient no. ki. */
285
286 kjh = 0; kl1 = kdim*ki;
287 for (kj=0; kj<=kder; kj++)
288 {
289
290 /* The value of the B-spline derivative is stored in tt while
291 kl2 steps through the idim components of this B-spline
292 coefficient. */
293
294 tt = ebder[kih++]; kl2 = kl1;
295 for (kl=0; kl<kdim; kl++,kjh++,kl2++)
296 {
297 sder[kjh] += scoef[kl2]*tt;
298 }
299 }
300 }
301
302 /* Free memory. */
303
304 /* If rational curve calculate the derivatives based on derivatives in
305 homogenous coordinates */
306
307 if (kind == 2 || kind == 4)
308 {
309 s6ratder(sder,pc1->idim,ider,eder,&kstat);
310 if (kstat<0) goto error;
311 freearray(sder);
312 }
313
314 freearray(ebder);
315
316 /* Successful computations. */
317
318 *jstat = 0;
319 goto out;
320
321 /* Not enough memory. */
322 err101: *jstat = -101;
323 s6err("S1227",*jstat,kpos);
324 goto out;
325
326 /* kdim less than 1. */
327 err102: *jstat = -102;
328 s6err("S1227",*jstat,kpos);
329 goto out;
330
331 /* Polynomial order less than 1. */
332 err110: *jstat = -110;
333 s6err("S1227",*jstat,kpos);
334 goto out;
335
336 /* Fewer B-splines than the order. */
337 err111: *jstat = -111;
338 s6err("S1227",*jstat,kpos);
339 goto out;
340
341 /* Error in knot vector.
342 (The first or last interval of the knot vector is empty.) */
343 err112: *jstat = -112;
344 s6err("S1227",*jstat,kpos);
345 goto out;
346
347 /* Illegal derivative requested. */
348 err178: *jstat = -178;
349 s6err("S1227",*jstat,kpos);
350 goto out;
351
352 /* Error in lower level routine. */
353
354 error: *jstat = kstat;
355 s6err("S1227",*jstat,kpos);
356 goto out;
357
358 out: return;
359 }
360