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: s1513.c,v 1.2 2001-03-19 15:58:50 afr Exp $
45 *
46 */
47
48
49 #define S1513
50
51 #include "sislP.h"
52
53 #if defined(SISLNEEDPROTOTYPES)
54 void
s1513(SISLSurf * psurf,double qpoint[],double bvec[],int idim,SISLSurf ** rsurf,int * jstat)55 s1513 (SISLSurf * psurf, double qpoint[], double bvec[], int idim,
56 SISLSurf ** rsurf, int *jstat)
57 #else
58 void
59 s1513 (psurf, qpoint, bvec, idim, rsurf, jstat)
60 SISLSurf *psurf;
61 double qpoint[];
62 double bvec[];
63 int idim;
64 SISLSurf **rsurf;
65 int *jstat;
66 #endif
67 /*
68 *********************************************************************
69 *
70 * PURPOSE : To make the function whose zeroes describes the circular silhouette
71 * lines of a surface. The circular silhouette lines are defined
72 * by a point qpoint and an axis bvec.
73 * f(u,v) = N(u,v) x (P(u,v) - Q) . B = 0
74 *
75 *
76 *
77 * INPUT : psurf - Pointer to input surface
78 * qpoint - Point on axis vector.
79 * bvec - Spin axis vector.
80 B
81 *
82 *
83 *
84 * OUTPUT : rsurf - The resulting surface
85 * jstat - status messages
86 * > 0 : warning
87 * = 0 : ok
88 * < 0 : error
89 *
90 *
91 * METHOD : We first make the appropriate knot vector, then we calulate
92 * parameter values for the interpolation, then the positions and normal
93 * vectors of the surface are calculated at these parameter
94 * values. The function f(U,V)
95 * is calculated at the parameter values and these f values are interpolated
96 * to put f, which describes the silhouette lines, into B-spline form.
97 *
98 * For nonrational surfaces the function whose zeroes form the
99 * silhouette is
100 *
101 * f = N(U,V) x (P(U,V) - Q) . B = (P xP ) x (P - Q) . B
102 * U V
103 *
104 * which, if P is degree m x n, is degree 3m-1 x 3n-1
105 * or order 3(m+1)-3 x 3(n+1)-3.
106 * In the rational case, setting P(U,V) = T(U,V) / w(U,V), we find
107 *
108 * 4
109 * P xP =(wT -w T) x (wT -w T) / w
110 * U V U U V V
111 *
112 * 3
113 * =(wT xT +w T xT+w TxT ) / w
114 * U V U V V U
115 *
116 * so we find the zeroes of
117 *
118 * f = (wT xT +w T xT+w TxT ) x (T - wQ) . B
119 * U V U V V U
120 *
121 * which, if P is degree m x n, is degree 4m-1 x 4n-1
122 * or order 4(m+1)-4 x 4(n+1)-4.
123 *
124 *
125 * REFERENCES :
126 *
127 *-
128 * CALLS :
129 *
130 * WRITTEN BY : Mike Floater, SI, 1990-11
131 * REVISED BY : Mike Floater, SI, 1991-04 for rational surfaces
132 *
133 *********************************************************************
134 */
135 {
136 double cross[3]; /* Cross product of norm and diff */
137 double diff[3]; /* Difference between surface point and eye point */
138 int kn1; /* Number of vertices of psurf in first par.dir */
139 int kk1; /* Order in psurf in first par.dir */
140 int kn2; /* Number of vertices of psurf in second par.dir */
141 int kk2; /* Order in psurf in second par.dir */
142 int kjkk1; /* Order of interpolated basis in first par.dir */
143 int kjkn1; /* Number of vertices in interpolated basis first.dr*/
144 int kjkk2; /* Order of interpolated basis in first par dir */
145 int kjkn2; /* Number of vertices in interpolated basis secnd.dr*/
146 int kdim; /* Number of dimesions in psurf */
147 int kstat; /* Local status variable */
148 int kpos = 0; /* Position indicator for errors */
149 int kzero = 0; /* Value 0 */
150 int kone = 1; /* Value 1 */
151 int cuopen; /* Open/Closed flag */
152 int ki, kj; /* Loop control variable */
153 int kp; /* Index of points put into conic equation */
154 int klfs = 0; /* Pointer into knot vector */
155 int klft = 0; /* Pointer into knot vector */
156
157 double *st1 = SISL_NULL; /* First knot vector is psurf */
158 double *st2 = SISL_NULL; /* Second knot vector is psurf */
159 double *sval1 = SISL_NULL; /* Array of values of surface put into torus eq. */
160 double *sval2 = SISL_NULL;
161 double *sval3 = SISL_NULL;
162
163 double *sgt1 = SISL_NULL; /* Knot vector in first parameter direction of
164 surface put into torus equation */
165 double *sgt2 = SISL_NULL; /* Knot vector in second parameter direction of
166 surface put into torus equation */
167 double sder[12]; /* Point on the surface */
168 double spar[2]; /* Current parameter pair */
169 double snorm[3]; /* Normal vector */
170 int i; /* a loop variable */
171 int ikind; /* kind of surface psurf is */
172 double tutv[3]; /* T_u x T_v in rational case */
173 double tvt[3]; /* T_v x T in rational case */
174 double ttu[3]; /* T x T_u in rational case */
175 SISLSurf *tempsurf = SISL_NULL; /* only used for rational surfaces */
176 double *par1=SISL_NULL;
177 double *par2=SISL_NULL;
178 int *der1=SISL_NULL;
179 int *der2=SISL_NULL;
180
181
182 if (idim != psurf->idim) goto err104;
183
184 /* Make local pointers */
185
186 kn1 = psurf->in1;
187 kk1 = psurf->ik1;
188 kn2 = psurf->in2;
189 kk2 = psurf->ik2;
190 kdim = psurf->idim;
191 st1 = psurf->et1;
192 st2 = psurf->et2;
193 ikind = psurf->ikind;
194
195 /* Test input */
196
197 if (kdim != 3)
198 goto err104;
199
200 if (ikind == 2 || ikind == 4)
201 {
202 /* A tricky way to evaluate the derivatives of the HOMOGENEOUS form
203 of psurf. In other words we need the derivs of T(u,v) where
204 p(u,v) = T(u,v) / w(u,v).
205 We should really have a separate evaluator for this
206 but I didn't want to mess around with the existing evaluator
207 which does the division automatically. MF 16/4/91
208 */
209
210 tempsurf = newSurf (kn1, kn2, kk1, kk2, st1, st2, psurf->rcoef, ikind - 1, kdim + 1, 0);
211 if (tempsurf == SISL_NULL) goto err171;
212
213 tempsurf->cuopen_1 = psurf->cuopen_1;
214 tempsurf->cuopen_2 = psurf->cuopen_2;
215 kjkk1 = 4 * kk1 - 4;
216 kjkk2 = 4 * kk2 - 4;
217 }
218 else
219 {
220 tempsurf = psurf;
221 kjkk1 = 3 * kk1 - 3;
222 kjkk2 = 3 * kk2 - 3;
223 }
224
225 /* Make description of knot array for interpolation in first parameter
226 direction */
227
228 s1381 (st1, kn1, kk1, &sgt1, &kjkn1, kjkk1, &kstat);
229 if (kstat < 0)
230 goto error;
231
232 /* Make parameter values and derivative indicators */
233
234 s1890(sgt1,kjkk1,kjkn1,&par1,&der1,&kstat);
235 if (kstat < 0)
236 goto error;
237
238
239 /* Make description of knot array for interpolation in second parameter
240 direction
241 */
242
243 s1381 (st2, kn2, kk2, &sgt2, &kjkn2, kjkk2, &kstat);
244 if (kstat < 0)
245 goto error;
246
247 /* Make parameter values and derivative indicators */
248
249 s1890(sgt2,kjkk2,kjkn2,&par2,&der2,&kstat);
250 if (kstat < 0)
251 goto error;
252
253 /* Allocate array for values of surface put into torus equation */
254
255 sval1 = newarray (kjkn1 * kjkn2, DOUBLE);
256 if (sval1 == SISL_NULL)
257 goto err101;
258
259 /* Calculate values to be interpolated */
260
261 /* Index of point to be stored */
262
263 kp = 0;
264
265 for (kj = 0; kj < kjkn2; kj++)
266 {
267
268 spar[1] = par2[kj];
269
270 for (ki = 0; ki < kjkn1; ki++)
271 {
272 /* Calculate values on 3-D surface */
273
274 spar[0] = par1[ki];
275
276 s1421 (tempsurf, 1, spar, &klfs, &klft, sder, snorm, &kstat);
277 if (kstat < 0)
278 goto error;
279
280 if (ikind == 2 || ikind == 4)
281 {
282 /* Calculate (wT xT +w T xT+w TxT )
283 U V U V V U
284 instead of normal to surface in rational case. */
285
286 s6crss (sder + 4, sder + 8, tutv);
287 s6crss (sder + 8, sder, tvt);
288 s6crss (sder, sder + 4, ttu);
289 for (i = 0; i < 3; i++)
290 {
291 snorm[i] = sder[3] * tutv[i] + sder[7] * tvt[i] + sder[11] * ttu[i];
292 diff[i] = sder[i] - sder[3] * qpoint[i];
293 }
294 }
295 else
296 {
297 s6diff (sder, qpoint, kdim, diff);
298 }
299
300 /* Make function f(u,v) = N(u,v) x ( P(u,v) - Q ) . B */
301
302 s6crss (snorm, diff, cross);
303 sval1[kp++] = s6scpr (cross, bvec, kdim);
304 }
305 }
306
307 cuopen = TRUE;
308
309 /* Interpolate in second parameter direction, the first parameter direction
310 is treated as a point of dimension kjkn1 */
311
312 s1891(par2,sval1,kjkn1,kjkn2,kone,der2,cuopen,sgt2,&sval2,&kjkn2,kjkk2,kzero,kzero,&kstat);
313 if (kstat < 0)
314 goto error;
315
316
317 /* Interpolate in first parameter direction, perform kjkn2 interpolations
318 of one dimensional data */
319
320 s1891(par1,sval2,kone,kjkn1,kjkn2,der1,cuopen,sgt1,&sval3,&kjkn1,kjkk1,kzero,kzero,&kstat);
321 if (kstat < 0)
322 goto error;
323
324 if (kstat < 0)
325 goto error;
326
327 *rsurf = SISL_NULL;
328 *rsurf = newSurf (kjkn1, kjkn2, kjkk1, kjkk2, sgt1, sgt2, sval3, 1, 1, 1);
329
330 if (*rsurf == SISL_NULL)
331 goto err171;
332
333 (*rsurf)->cuopen_1 = psurf->cuopen_1;
334 (*rsurf)->cuopen_2 = psurf->cuopen_2;
335
336 /* Ok ! */
337
338
339 *jstat = 0;
340 goto out;
341
342 /* Error in lower level function */
343
344 error:
345 *jstat = kstat;
346 s6err ("s1513", *jstat, kpos);
347 goto out;
348
349 /* Error in space allocation */
350 err101:
351 *jstat = -101;
352 s6err ("s1513", *jstat, kpos);
353 goto out;
354
355 /* Error in dimension */
356 err104:
357 *jstat = -104;
358 s6err ("s1513", *jstat, kpos);
359 goto out;
360
361 /* Could not create surfe */
362
363 err171:
364 *jstat = -171;
365 s6err ("s1513", *jstat, kpos);
366 goto out;
367
368 out:
369
370 /* Release allocated arrays */
371
372 if (sgt1 != SISL_NULL)
373 freearray (sgt1);
374 if (sgt2 != SISL_NULL)
375 freearray (sgt2);
376 if (sval1 != SISL_NULL)
377 freearray (sval1);
378 if (sval2 != SISL_NULL)
379 freearray (sval2);
380 if (sval3 != SISL_NULL)
381 freearray (sval3);
382 if (par1 != SISL_NULL)
383 freearray(par1);
384 if (par2 != SISL_NULL)
385 freearray(par2);
386 if (der1 != SISL_NULL)
387 freearray(der1);
388 if (der2 != SISL_NULL)
389 freearray(der2);
390 if ((ikind == 2 || ikind == 4) && (tempsurf != SISL_NULL))
391 freeSurf (tempsurf);
392
393 return;
394 }
395