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: s1710.c,v 1.3 2001-03-19 15:58:52 afr Exp $
45 *
46 */
47
48
49 #define S1710
50
51 #include "sislP.h"
52
53 #if defined(SISLNEEDPROTOTYPES)
54 void
s1710(SISLCurve * pc1,double apar,SISLCurve ** rcnew1,SISLCurve ** rcnew2,int * jstat)55 s1710 (SISLCurve * pc1, double apar, SISLCurve ** rcnew1, SISLCurve ** rcnew2, int *jstat)
56 #else
57 void
58 s1710 (pc1, apar, rcnew1, rcnew2, jstat)
59 SISLCurve *pc1;
60 double apar;
61 SISLCurve **rcnew1;
62 SISLCurve **rcnew2;
63 int *jstat;
64 #endif
65 /*
66 *******************************************************************
67 *
68 *********************************************************************
69 *
70 * PURPOSE : Subdivide a B-spline curve at a given parameter-value.
71 * NOTE: When the curve is periodic (ie. cuopen is
72 * set to SISL_CRV_PERIODIC and open basis
73 * with order-mult repeated knots and coeffic.)
74 * this function will return only ONE curve
75 * through rcnew1. This curve is the same
76 * geometric curve as pc1, but is represented on
77 * a closed basis with start and end at pc(apar).
78 * Cuopen is set to SISL_CRV_CLOSED.
79 * jstat equals 2 when this occurs.
80 *
81 *
82 *
83 * INPUT : pc1 - SISLCurve to subdivide.
84 * apar - Parameter-value at which to subdivide.
85 *
86 *
87 *
88 * OUTPUT : rcnew1 - First part of the subdivided curve.
89 * rcnew2 - Second part of the subdivided curve.
90 * If the parameter value is at the end of a
91 * curve SISL_NULL pointers might be returned
92 * jstat - status messages
93 * = 2 : pc1 periodic, rcnew2=SISL_NULL
94 * = 5 : parameter value at end of
95 * curve, rcnew1=SISL_NULL or
96 * rcnew2=SISL_NULL.
97 * > 0 : warning
98 * = 0 : ok
99 * < 0 : error
100 *
101 *
102 * METHOD : Inserting knots at the subdividing-point until
103 * we have a ktuple-knot. Then we may separate the
104 * curve into two parts.
105 *
106 *
107 * REFERENCES :
108 *
109 *-
110 * CALLS : newCurve - Allocate space for a new curve-object.
111 * freeCurve - Free space occupied by given curve-object.
112 * S1700.C - Making the knot-inserten-transformation matrix.
113 *
114 * WRITTEN BY : Arne Laksaa, SI, 88-06.
115 * Revised by : Tor Dokken, SI, 89-03. Can be used for converting
116 * from closed to open description
117 * MODIFIED BY : Mike Floater, SI, 91-01. Subdivide rational curves.
118 * MODIFIED BY : Ulf J. Krystad, SI, 92-01. Subdivide periodic crvs.
119 * MODIFIED BY : Arne Laksaa, SI, 92-09. Move apar to closest knot
120 * if it is close to a knot. Using D(N)EQUAL().
121 * Revised by : Paal Fugelli, SINTEF, Oslo, Norway, 94-08. Added error propagation.
122 *
123 *
124 **********************************************************************/
125 {
126 int kind = pc1->ikind; /* Type of curve pc1 is. */
127 int kstat; /* Local status variable. */
128 int kpos = 0; /* Position of error. */
129 int kmy; /* An index to the knot-vector. */
130 int kv, kv1; /* Number of knots we have to insert. */
131 int kpl, kfi, kla; /* To posisjon elements in trans.-matrix. */
132 int kk = pc1->ik; /* Order of the input curve. */
133 int kn = pc1->in; /* Number of the vertices in input curves. */
134 int kdim = pc1->idim; /* Dimensjon of the space in whice
135 * the curve lies. */
136 int kn1, kn2; /* Number of vertices in the new curves. */
137 int knum; /* Number of knots less and equal than
138 the intersection point. */
139 int ki, ki1; /* Control variable in loop. */
140 int kj, kj1, kj2; /* Control variable in loop. */
141 int newkind = 1; /* Type of curve the subcurves are */
142 double *s1, *s2, *s3, *s4; /* Pointers used in loop. */
143 double *st1 = SISL_NULL; /* The first new knot-vector. */
144 double *st2 = SISL_NULL; /* The second new knot-vector. */
145 double *salfa = SISL_NULL; /* A line of the trans.-matrix. */
146 double *scoef; /* Pointer to vertices. */
147 double *scoef1 = SISL_NULL; /* The first new vertice. */
148 double *scoef2 = SISL_NULL; /* The second new vertice. */
149 SISLCurve *q1 = SISL_NULL; /* Pointer to new curve-object. */
150 SISLCurve *q2 = SISL_NULL; /* Pointer to new curve-object. */
151 int incr; /* Number of extra knots copied
152 * during periodicity */
153 int mu; /* Multiplisity at the k'th knot */
154 int kleft = kk-1; /* Knot navigator */
155 double delta; /* Period size in knot array. */
156 double salfa_local[5]; /* Local help array. */
157
158 *rcnew1 = SISL_NULL;
159 *rcnew2 = SISL_NULL;
160
161 /* if pc1 is rational, do subdivision in homogeneous coordinates */
162 /* just need to set up correct dim and kind for the new curves at end of routine */
163 if (kind == 2 || kind == 4)
164 {
165 scoef = pc1->rcoef;
166 kdim++;
167 newkind++;
168 }
169 else
170 {
171 scoef = pc1->ecoef;
172 }
173
174 /* Check that we have a curve to subdivide. */
175
176 if (!pc1)
177 goto err150;
178
179
180 /* Periodic curve treatment, UJK jan 92--------------------------------- */
181 if (pc1->cuopen == SISL_CRV_PERIODIC)
182 {
183 delta = (pc1->et[kn] - pc1->et[kk - 1]);
184
185 /* Check that the intersection point is an interior point. */
186 /*if (apar < *(pc1->et) || apar > *(pc1->et + kn + kk - 1))*/
187 if ((apar < pc1->et[0] && DNEQUAL(apar, pc1->et[0])) ||
188 (apar > pc1->et[kn+kk-1] && DNEQUAL(apar, pc1->et[kn+kk-1])))
189 goto err158;
190
191 /* If inside the knot vector, but outside well define
192 intervall, we shift the parameter value one period. */
193 if (apar < pc1->et[kk - 1] && DNEQUAL(apar, pc1->et[kk - 1]))
194 apar += delta;
195 if (apar > pc1->et[kn] || DEQUAL(apar, pc1->et[kn]))
196 apar -= delta;
197
198 /* Now we create a new curve that is a copy of pc1,
199 but with the period repeated once,
200 this allows us to pick a whole period. */
201
202 /* Get multiplisity at start of full basis interval */
203 mu = s6knotmult(pc1->et, kk, kn, &kleft, pc1->et[kk-1], &kstat);
204 if (kstat < 0) goto err153;
205 if (mu >= kk) goto errinp;
206
207 /* Copy ----------------------------------- */
208 incr = kn - kk + mu;
209 if ((scoef1 = newarray ((kn + incr) * kdim, double)) == SISL_NULL)
210 goto err101;
211 if ((st1 = newarray (kn + kk + incr, double)) == SISL_NULL)
212 goto err101;
213
214 memcopy (scoef1, scoef, kn * kdim, double);
215 memcopy (st1, pc1->et, kn + kk, double);
216 memcopy (scoef1 + kn * kdim, scoef + (kk - mu) * kdim,
217 incr * kdim, double);
218
219
220 for (ki = 0; ki < incr; ki++)
221 st1[kn + kk + ki] = st1[kn + kk + ki - 1] +
222 (st1[2*kk - mu + ki] - st1[2*kk - mu + ki - 1]);
223
224 if ((q1 = newCurve (kn + incr, kk, st1, scoef1,
225 newkind, pc1->idim, 2)) == SISL_NULL)
226 goto err101;
227 q1->cuopen = SISL_CRV_OPEN;
228
229 /* Pick part (one period)------------------ */
230 s1712 (q1, apar, apar + delta,
231 rcnew1, &kstat);
232 if (kstat < 0)
233 goto err153;
234 freeCurve (q1);
235 if (*rcnew1)
236 (*rcnew1)->cuopen = SISL_CRV_CLOSED;
237
238 /* Finished, exit */
239 *jstat = 2;
240 goto out;
241
242 }
243
244 /* End of periodic curve treatment, UJK jan 92------------- */
245
246 /* Check that the intersection point is an interior point. */
247 /* Changed by UJK and later ALA*/
248 /*if (apar < *(pc1->et) || apar > *(pc1->et+kn+kk-1)) goto err158; */
249
250 if ((apar < pc1->et[kk - 1] && DNEQUAL(apar, pc1->et[kk - 1]))||
251 (apar > pc1->et[kn] && DNEQUAL(apar, pc1->et[kn])))
252 goto err158;
253
254 /* Allocate space for the kk elements which may not be zero in eache
255 line of the basic transformation matrix.*/
256
257 if (kk > 5)
258 {
259 if ((salfa = newarray (kk, double)) == SISL_NULL) goto err101;
260 }
261 else salfa = salfa_local;
262
263
264 /* Find the number of the knots which is smaller or like
265 the intersection point, and how many knots we have to insert.*/
266
267 s1 = pc1->et;
268 kv = kk; /* The maximum number of knots we may have to insert. */
269
270 if ((apar > s1[0] && DNEQUAL(apar, s1[0])) &&
271 (apar < s1[kn+kk-1] && DNEQUAL(apar, s1[kn+kk-1])))
272 {
273 /* Using binear search*/
274 kj1=0;
275 kj2=kk+kn-1;
276 knum = (kj1+kj2)/2;
277 while (knum != kj1)
278 {
279 if ((s1[knum] < apar) && DNEQUAL (s1[knum], apar))
280 kj1=knum;
281 else
282 kj2=knum;
283 knum = (kj1+kj2)/2;
284 }
285 knum++; /* The smaller knots. */
286
287 while (DEQUAL (s1[knum], apar))
288 /* The knots thats like the intersection point. */
289 {
290 apar = s1[knum];
291 knum++;
292 kv--;
293 }
294 }
295 else if (DEQUAL(apar,s1[0]))
296 {
297 apar = s1[0];
298 knum = 0;
299 while (s1[knum] == apar)
300 /* The knots thats like the intersection point. */
301 knum++;
302 }
303 else if (DEQUAL(apar,s1[kn+kk-1]))
304 {
305 apar = s1[kn+kk-1];
306 knum = kn+kk-1;
307 while (s1[knum-1] == apar)
308 /* The knots thats like the intersection point. */
309 knum--;
310 }
311
312 /* Find the number of vertices in the two new curves. */
313
314 kn1 = knum + kv - kk;
315 kn2 = kn + kk - knum;
316
317
318
319 /* Allocating the new arrays to the two new curves. */
320
321 if (kn1 > 0)
322 {
323 if ((scoef1 = newarray (kn1 * kdim, double)) == SISL_NULL)
324 goto err101;
325 if ((st1 = newarray (kn1 + kk, double)) == SISL_NULL)
326 goto err101;
327 }
328 if (kn2 > 0)
329 {
330 if ((scoef2 = newarray (kn2 * kdim, double)) == SISL_NULL)
331 goto err101;
332 if ((st2 = newarray (kn2 + kk, double)) == SISL_NULL)
333 goto err101;
334 }
335
336
337 /* Copying the knotvectors, all but the intersection point from
338 the old curve to the new curves */
339
340 memcopy (st1, pc1->et, kn1, double);
341 memcopy (st2 + kk, pc1->et + knum, kn2, double);
342
343
344 /* Updating the knotvectors by inserting a k-touple knot in
345 the intersection point at each curve.*/
346
347 for (s2 = st1 + kn1, s3 = st2, s4 = s3 + kk; s3 < s4; s2++, s3++)
348 *s2 = *s3 = apar;
349
350
351 /* Copying the coefisientvectors to the new curves.*/
352
353 memcopy (scoef1, scoef, kdim * kn1, double);
354 memcopy (scoef2, scoef + kdim * (knum - kk), kdim * kn2, double);
355
356
357 /* Updating the coefisientvectors to the new curves.*/
358
359 /* Updating the first curve. */
360 knum -= kk - 1;
361 for (ki=max(0, knum), kv1=max(0,-knum), s1=scoef1+ki*kdim; ki < kn1; ki++)
362 {
363 /* Initialising:
364 knum = knum-kk+1, Index of the first vertice to change.
365 ki = knum, Index of the vertices we are going to
366 change. Starting with knum, but if
367 knum is negativ we start at zero.
368 kv1 = 0, Number if new knots between index ki
369 and ki+kk. We are starting one below
370 becase we are counting up before using
371 it. If knum is negativ we are not
372 starting at zero but at -knum.
373 s1=scoef1+ki*kdim,SISLPointer at the first vertice to
374 change. */
375
376
377 /* Using the Oslo-algorithm to make a transformation-vector
378 from the old vertices to one new vertice. */
379
380 kmy = ki;
381 s1700 (kmy, kk, kn, ++kv1, &kpl, &kfi, &kla, pc1->et, apar, salfa, &kstat);
382 if (kstat)
383 goto err153;
384
385
386 /* Compute the kdim vertices with the same "index". */
387
388 for (kj = 0; kj < kdim; kj++, s1++)
389 for (*s1 = 0, kj1 = kfi, kj2 = kfi + kpl; kj1 <= kla; kj1++, kj2++)
390 *s1 += salfa[kj2] * scoef[kj1 * kdim + kj];
391 }
392
393 /* And the second curve. */
394
395 for (ki1 = min (kn1 + kv - 1, kn + kv), s1 = scoef2; ki < ki1; ki++)
396 {
397 /* Initialising:
398 ki1 = kn1+kv-1, the index of the vertice next to the
399 last vertice we have to change.
400 If we do not have so many vertices,
401 we have to use the index next to the
402 last vertice we have, kn+kv.
403 s1=scoef2 Pointer at the first vertice to
404 change. */
405
406
407 /* Using the Oslo-algorithm to make a transformation-vector
408 from the old vertices to one new vertice. */
409
410 s1700 (kmy, kk, kn, kv1--, &kpl, &kfi, &kla, pc1->et, apar, salfa, &kstat);
411 if (kstat)
412 goto err153;
413
414
415 /* Compute the kdim vertices with the same "index". */
416
417 for (kj = 0; kj < kdim; kj++, s1++)
418 for (*s1 = 0, kj1 = kfi, kj2 = kfi + kpl; kj1 <= kla; kj1++, kj2++)
419 *s1 += salfa[kj2] * scoef[kj1 * kdim + kj];
420 }
421
422
423 /* Allocating new curve-objects.*/
424
425 if (kn1 > 0)
426 q1 = newCurve (kn1, kk, st1, scoef1, newkind, pc1->idim, 2);
427
428 if (kn2 > 0)
429 q2 = newCurve (kn2, kk, st2, scoef2, newkind, pc1->idim, 2);
430
431 if (q1 == SISL_NULL && q2 == SISL_NULL) goto err101;
432
433 /* Updating output. */
434
435 *rcnew1 = q1;
436 *rcnew2 = q2;
437 if (q1 == SISL_NULL || q2 == SISL_NULL)
438 *jstat = 5; /* The curve is subdivided in an endpoint. */
439 else
440 *jstat = 0;
441 goto out;
442
443
444 /* Error. Error in low level routine. */
445
446 err153:
447 *jstat = kstat;
448 s6err ("s1710", *jstat, kpos);
449 goto outfree;
450
451 /* Error. Error in input */
452 errinp:
453 *jstat = -154;
454 s6err ("s1710", *jstat, kpos);
455 goto outfree;
456
457
458 /* Error. No curve to subdivide. */
459
460 err150:
461 *jstat = -150;
462 s6err ("s1710", *jstat, kpos);
463 goto out;
464
465
466 /* Error. The parameter value is outside the curve. */
467
468 err158:
469 *jstat = -158;
470 s6err ("s1710", *jstat, kpos);
471 goto out;
472
473
474 /* Error. Allocation error, not enough memory. */
475
476 err101:
477 *jstat = -101;
478 s6err ("s1710", *jstat, kpos);
479 goto outfree;
480
481
482 outfree:
483 if (q1)
484 freeCurve (q1);
485
486 if (q2)
487 freeCurve (q2);
488
489
490 /* Free local used memory. */
491
492 out:
493 if (!q1)
494 {
495 if (st1)
496 freearray (st1);
497 if (scoef1)
498 freearray (scoef1);
499 }
500
501 if (!q2)
502 {
503 if (st2)
504 freearray (st2);
505 if (scoef2)
506 freearray (scoef2);
507 }
508
509 if (kk > 5 && salfa)
510 freearray (salfa);
511 return;
512 }
513