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