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: s1612.c,v 1.2 2001-03-19 15:58:51 afr Exp $
45 *
46 */
47
48
49 #define S1612
50
51 #include "sislP.h"
52
53 #if defined(SISLNEEDPROTOTYPES)
54 void
s1612(SISLCurve * pc,double aepsge,double ** gpoint,int * jnbpnt,int * jleng,int * jstat)55 s1612(SISLCurve *pc,double aepsge,double **gpoint,int *jnbpnt,int *jleng,int *jstat)
56 #else
57 void s1612(pc,aepsge,gpoint,jnbpnt,jleng,jstat)
58 SISLCurve *pc;
59 double aepsge;
60 double **gpoint;
61 int *jnbpnt;
62 int *jleng;
63 int *jstat;
64 #endif
65 /*
66 *********************************************************************
67 *
68 * PURPOSE : To calculate a set of points on a B-spline curve.
69 * The straight lines between the points will not deviate
70 * more than aepsge from the B-spline curve at any point.
71 *
72 * INPUT : pc - The input B-spline curve.
73 * aepsge - Geometry resolution, maximum distance allowed between
74 * the curve and the straight lines to be calculated.
75 *
76 * INPUT/OUTPUT :
77 * jnbpnt - No. of calculated points until now
78 * jleng - no. of allocated doubles in gpoint until now
79 * gpoint - Calculated points
80 *
81 * OUTPUT : jstat - status messages
82 * > 0 : warning
83 * = 0 : ok
84 * < 0 : error
85 *
86 * METHOD : This routine recursively calls itself.
87 * First the curve is split at kk-1 internal multiple knots.
88 * Each curve segment is then treated:
89 * The distance between the control polygon and a straight
90 * line from the first to the last vertex is found. If the
91 * distance is grater than aepsge, the curve is split again
92 * as follows:
93 * If any internal knots exist the curve is split in the
94 * middlemost knot; if not, the curve is split in the middle of
95 * the parameter interval.
96 * If the distance is less than aepsge the last vertex is
97 * saved in array.
98 *
99 * ASSUMTIONS :kk multiple knots in start and end of parameterintervall.
100 *
101 * USE :
102 *
103 * REFERENCES :
104 *
105 *-
106 * CALLS : s1840,s1710,s1612(recursive),s1235,s6err
107 *
108 *
109 * WRITTEN BY : Qyvind Hjelle, SI, Oslo, Norway. 31. Jan 1989
110 *
111 *********************************************************************
112 */
113 {
114 int kstat; /* Status variable */
115 int kn; /* The number of B-splines, i.e., the dimension of
116 the spline space associated with the knot
117 vector. */
118 int kk; /* The polynomial order of the curve. */
119 int kdim; /* The dimension of the space in which the curve
120 lies. Equivalently, the number of components
121 of each B-spline coefficient. */
122 int ki; /* Local counter */
123 int kleft; /* Pointer into knot vector array */
124 int knbreak=0; /* No. of kk-1 multiple knots */
125 int knbpnt; /* Local for jnbpnt */
126 int kleng; /* Local for jleng */
127 int kvlast; /* Position of last vertex in vertex array */
128 int kpos=0; /* Position of error */
129
130 double *spoint=SISL_NULL;/* Pointer to array of points */
131 double tdist; /* Distance */
132 double tpar; /* A parameter value of the curve */
133 double *sbreak = SISL_NULL; /* Array containing kk-1 multiple knot */
134 double *st; /* Pointer to the first element of the knot vector
135 of the curve. The knot vector has [kn+kk]
136 elements. */
137
138
139 SISLCurve *qcnew1=SISL_NULL; /* Pointer to first new curve-object */
140 SISLCurve *qcnew2=SISL_NULL; /* Pointer to second new curve-objec */
141
142
143 /* Check input */
144
145 if (aepsge <= (double)0.0) goto err120;
146
147 /* Make locals */
148
149 spoint = *gpoint;
150 knbpnt = *jnbpnt;
151 kleng = *jleng;
152
153 /* Describe curve with local parameters. */
154
155 kn = pc -> in;
156 kk = pc -> ik;
157 kdim = pc -> idim;
158 st = pc -> et;
159
160 /* Find all kk-1 multiple knots including start and endpoint */
161
162 s1235(st,kn,kk,&knbreak,&sbreak,&kstat);
163 if (kstat < 0) goto error;
164
165 /* Always split curve at kk-1 multiple knots */
166
167 if (knbreak > 2)
168 {
169 for (ki=1; ki<knbreak-1; ki++)
170 {
171 tpar = sbreak[ki];
172
173 s1710(pc,tpar,&qcnew1,&qcnew2,&kstat);
174 if(kstat < 0) goto error;
175
176 /* Recursion on first part */
177
178 if (qcnew1)
179 {
180 s1612(qcnew1,aepsge,&spoint,&knbpnt,&kleng,&kstat);
181 if (kstat < 0) goto error;
182 }
183
184 /* Recursion on second part */
185
186 if (qcnew2)
187 {
188 s1612(qcnew2,aepsge,&spoint,&knbpnt,&kleng,&kstat);
189 if (kstat < 0) goto error;
190 }
191 }
192 }
193 else
194 /* If no internal kk-1 multiple knots */
195 {
196 /* Find distance between control polygon and a straight line beetween
197 first and last vertex*/
198
199 s1840(pc,&tdist,&kstat);
200 if (kstat < 0) goto error;
201
202 if (tdist < aepsge)
203 {
204 /* Save last vertex */
205
206 kvlast = (kn-1) * kdim;
207
208 knbpnt += 1;
209
210 /* Allocate place for 100 more points in gpoint if not enought place */
211
212 if (kleng < (knbpnt+1)*kdim )
213 {
214 kleng += 100*kdim;
215 spoint = increasearray(spoint,kleng,DOUBLE);
216 if (!spoint) goto err101;
217 }
218
219 /* Save point in array */
220
221 memcopy (&spoint[(knbpnt-1)*kdim],&pc->ecoef[kvlast],kdim,DOUBLE);
222 }
223 else
224 {
225 /* Split curve in two parts. If any internal knots, split in the
226 midlemost */
227
228 tpar = (st[0] + st[kn+kk-1]) / (double)2.0;
229
230 if (kn > kk)
231 {
232 /* Localize the parameter value */
233
234 kleft = 0;
235 s1219(st,kk,kn,&kleft,tpar,&kstat);
236 if (kstat < 0) goto error;
237
238 /* Find the knot nearest to the midle of the parameter interval */
239
240 if (fabs(tpar-st[kleft]) < fabs(st[kleft+1]-tpar))
241 {
242 tpar = st[kleft];
243 }
244 else
245 {
246 tpar = st[kleft+1];
247 }
248 }
249
250 s1710(pc,tpar,&qcnew1,&qcnew2,&kstat);
251 if(kstat < 0) goto error;
252
253 /* Recursion on first part */
254
255 if (qcnew1)
256 {
257 s1612(qcnew1,aepsge,&spoint,&knbpnt,&kleng,&kstat);
258 if (kstat < 0) goto error;
259 }
260
261 /* Recursion on second part */
262
263 if (qcnew2)
264 {
265 s1612(qcnew2,aepsge,&spoint,&knbpnt,&kleng,&kstat);
266 if (kstat < 0) goto error;
267 }
268 }
269 }
270
271 *gpoint = spoint;
272 *jnbpnt = knbpnt;
273 *jleng = kleng;
274 *jstat = 0;
275 goto out;
276
277 /* Error in memory allocation */
278
279 err101:
280 *jstat = -101;
281 s6err("s1612",*jstat,kpos);
282 goto out;
283
284 /* Error in input. Relative tollerance <=0 */
285
286 err120:
287 *jstat = -120;
288 s6err("s1612",*jstat,kpos);
289 goto out;
290
291
292 /* Error in lower level function */
293
294 error:
295 *jstat = kstat;
296 s6err("s1612",*jstat,kpos);
297 goto out;
298
299
300 out:
301
302 /* Free space occupied by local arrays and objects. */
303
304 if (sbreak) freearray(sbreak);
305 if (qcnew1) freeCurve(qcnew1);
306 if (qcnew2) freeCurve(qcnew2);
307
308 return;
309 }
310
311
312