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