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: s1309.c,v 1.3 2001-03-19 15:58:43 afr Exp $
45 *
46 */
47 #define S1309
48
49 #include "sislP.h"
50
51 #if defined(SISLNEEDPROTOTYPES)
52 double
s1309(double epnt[],double edir[],double eimpli[],int ideg,int * jstat)53 s1309(double epnt[],double edir[],double eimpli[],int ideg,int *jstat)
54 #else
55 double s1309(epnt,edir,eimpli,ideg,jstat)
56 double epnt[];
57 double edir[];
58 double eimpli[];
59 int ideg;
60 int *jstat;
61 #endif
62 /*
63 *********************************************************************
64 *
65 * PURPOSE : To find the shortest distance between a point and the
66 * projection of the point along a direction vector on to
67 * an implicit represented surface. The distance will
68 * have a sign, positive if the projection is along the
69 * direction of the direction vector, else it is negative.
70 * For torus surfaces the projection direction is not
71 * used: The closest point is calculated and the distance to
72 * this returned.
73 * If ideg==1003,1004,1005, the the difference between PI/2 and the
74 * angle between the two vectors which whose scalar product
75 * wants to be as close as possible to zero, is found.
76 *
77 * INPUT : epnt - The point to be projected
78 * edir - Projection direction. Not used if ideg==1003,1004,1005
79 * eimpli - Description of implicit surface
80 * ideg - Degree of implicit surface
81 * ideg=1: Plane
82 * ideg=2; Quadric surface
83 * ideg=1001: Torus surface
84 * ideg=1003: Silhouette line parallel projection
85 * ideg=1004: Silhouette line perspective projection
86 * ideg=1005: Silhouette line circular projection
87 *
88 *
89 * OUTPUT : s1309 - The distance, positive if along positive direction
90 * of edir, negative if along negative direction of
91 * edir.
92 * If ideg==1003,1004,1005, the the difference between PI/2 and the
93 * angle between the two vectors which whose scalar product
94 * wants to be as close as possible to zero, is found (in radians).
95 * jstat - status messages
96 * = 0 : ok, iteration converged
97 * < 0 : error
98 *
99 *
100 * METHOD :
101 *
102 * USE : The function can be used for preparing the calculation
103 * of the projected point by adding the product of the distance
104 * and the normalized versions of edir to epnt.
105 * REFERENCES :
106 *
107 *-
108 * CALLS :
109 *
110 * WRITTEN BY : Tor Dokken, SI, Oslo, Norway, 4-July-1988
111 * Revised by : Mike Floater, SI, 1991-01
112 * Improved the routine for parallel silhouettes (ideg=1003) and
113 * added perspective and circular silhouettes (ideg=1004,ideg=1005)
114 * Changed by : Per OEyvind Hvidsten, SINTEF, 11-94.
115 * Added initialization of tcurdst in declaration.
116 *
117 *********************************************************************
118 */
119 {
120 double sdir[3]; /* Normilized direction vector */
121 double tb1=0.0,ta11=0.0,ta12=0.0; /* Dummy variables */
122 double tsum,t1,t2,tdum1;/* Dummy variables */
123 double tcurdst=0.0; /* The distance */
124 double sq[4]; /* Array used for temporary results */
125 int ksize; /* Number of doubles for storage of derivatives
126 and normal vector */
127 int ksizem3; /* ksize - 3 */
128 int kstat; /* Local status variable */
129 int kdim=3; /* Dimesnon of 3-D space */
130 int ki,kj,kl,kp; /* Control variables in loop */
131 int kpos=1; /* Position of error */
132
133
134
135 /* If ideg=1,2 or 1001 then only derivatives up to second order
136 are calculated, then 18 doubles for derivatives and 3 for the
137 normal vector are to be used for calculation of points in the
138 spline surface. For ideg=1003,1004,1005 we have a silhouette curve and
139 derivatives up to the third are to be calculated,
140 thus 30 +3 a total of 33 doubles are to be calculated */
141
142 if (ideg==1003 || ideg==1004 || ideg==1005)
143 {
144 ksize = 33;
145 }
146 else
147 {
148 ksize = 21;
149 }
150
151 ksizem3 = ksize -3;
152 (void)s6norm(edir,kdim,sdir,&kstat);
153 if (kstat < 0) goto error;
154
155 if (ideg==1)
156 {
157 /* Put parametric representation of projection line into
158 implicit equation of plane */
159 tb1 = s6scpr(eimpli,epnt,kdim);
160 ta11 = s6scpr(eimpli,sdir,kdim);
161 if( ta11 == (double)0.0) goto war02;
162 tcurdst = -(eimpli[3]+tb1)/ta11;
163 }
164 else if (ideg==2)
165 {
166
167 /* Find distance from the new point to the implicit surface, by
168 intersecting the straight line throught the point and with
169 sdir as normalized direction vector, normalized version of 3
170 first coordinates of sq the implicit surface and P0 = (P,1).
171
172 This problem can be written:
173
174 T T 2 T
175 P0 A P0 + 2 t P0 A sdir + t (sdir,0) A (sdir,0) = 0
176
177 T
178 We have to calulate calculate tb1=P0 A P and qs=P0 A, thus:
179 T 2 T
180 tb1 + 2 t sq sdir + t (sdir,0) A (sdir,0) = 0
181
182 Thus the first step is to calculate qs = A (P,1)
183 */
184 for (ki=0;ki<4;ki++)
185 {
186 tsum = eimpli[12+ki];
187 for (kj=0,kl=ki ; kj<3 ; kj++,kl+=4)
188 {
189 tsum +=(eimpli[kl]*epnt[kj]);
190 }
191 sq[ki] = tsum;
192 }
193
194 tb1 = s6scpr(epnt,sq,kdim) + sq[3];
195
196 ta11 = (double)2.0*s6scpr(sq,sdir,kdim);
197
198 ta12 = (double)0.0;
199 for (ki=0,kl=0;ki<3;ki++,kl+=4)
200 {
201 kp = kl;
202 for (kj=0;kj<3;kj++,kp++)
203 {
204 ta12 += sdir[ki]*eimpli[kp]*sdir[kj];
205 }
206 }
207
208 /* Now our equation system is:
209 2
210 ta12 t + ta11 t + tb1 = 0
211
212 we want the root with the smallest absolute value:
213 2
214 t = (-ta11 +/- sqrt(ta11 -4ta12 tb1))/(2ta12)
215 */
216 if (DNEQUAL(ta12,(double)0.0))
217 {
218 tdum1 = ta11*ta11 - (double)4.0*ta12*tb1;
219 if (tdum1 < DZERO) goto war02;
220 tdum1 = sqrt(tdum1);
221 t1 = (-ta11 + tdum1)/((double)2.0*ta12);
222 t2 = (-ta11 - tdum1)/((double)2.0*ta12);
223 if (fabs(t1)<fabs(t2))
224 {
225 tcurdst = t1;
226 }
227 else
228 {
229 tcurdst = t2;
230 }
231 }
232 else if(DNEQUAL(ta11,(double)0.0))
233 {
234 tcurdst = tb1/ta11;
235 }
236 else
237 {
238 /* Unsolvable system */
239 goto war02;
240 }
241 }
242 else if (ideg==1001)
243 {
244 /* Torus surface */
245
246 double *scentr; /* The center of the torus */
247 double *snorm; /* The normal of the torus symmetry plane */
248 double tbigr; /* The big radius of the torus */
249 double tsmalr; /* The small radius of the torus */
250 double sdum1[3]; /* Temporary storage for point */
251 double sdum2[3]; /* Temporary storage for point */
252 double tproj; /* Projection of vector onto snorm */
253
254
255 scentr = eimpli;
256 snorm = eimpli+3;
257 tbigr = *(eimpli+6);
258 tsmalr = *(eimpli+7);
259
260 /* Find projection of vector from torus center on to torus axis */
261 s6diff(epnt,scentr,kdim,sdum1);
262 tproj = s6scpr(sdum1,snorm,kdim);
263
264 /* Project vector from torus center to epnt onto torus plane */
265 for (ki=0;ki<kdim;ki++)
266 sdum2[ki] = sdum1[ki] - tproj*snorm[ki];
267 (void)s6norm(sdum2,kdim,sdum2,&kstat);
268 if (kstat<0) goto error;
269
270 /* Find vector from torus circle to epnt */
271 for (ki=0;ki<kdim;ki++)
272 sdum1[ki] = sdum1[ki] - tbigr*sdum2[ki];
273
274 /* Find length of this vector and compare with tsmalr */
275
276 tcurdst = fabs(s6length(sdum1,kdim,&kstat)-tsmalr);
277 if (kstat<0) goto error;
278 }
279 else if (ideg==1003)
280 {
281 /* Silhouette line/curve */
282
283 double sdum1[3]; /* Temporary storage for point */
284
285 (void)s6norm(epnt+ksizem3,kdim,sdum1,&kstat);
286 if (kstat<0) goto error;
287
288 /* eimpli[0,1,2] is assumed to, be normalized */
289
290 t1 = s6scpr(sdum1,eimpli,kdim);
291
292 /* t1 now contains the cosine of the angle between the view direction
293 and the normal vector. This is Equal to sin of PI/2 minus the angle
294 between the vectors, thus the actual angle can be calulated */
295
296 tcurdst = asin(t1);
297 tcurdst = fabs(tcurdst);
298 }
299
300 else if (ideg==1004)
301 {
302 /* Perspective silhouette line/curve */
303
304 double sdum1[3],sdum2[3]; /* Temporary storage for point */
305
306 s6diff(epnt,eimpli,kdim,sdum1);
307 (void)s6norm(sdum1,kdim,sdum2,&kstat);
308 /* OK if sdum1 is zero -- tcurdst will be zero as well */
309 (void)s6norm(epnt+ksizem3,kdim,sdum1,&kstat);
310
311
312 t1 = s6scpr(sdum1,sdum2,kdim);
313
314 /* t1 now contains the cosine of the angle between the direction of the
315 point epnt relative to the eyepoint E (in eimpli)
316 and the normal vector. This is Equal to sin of PI/2 minus the angle
317 between the vectors, thus the actual angle can be calulated */
318
319 tcurdst = asin(t1);
320 tcurdst = fabs(tcurdst);
321 }
322
323 else if (ideg==1005)
324 {
325 /* Circular silhouette line/curve */
326
327 double sdum1[3],sdum2[3]; /* Temporary storage for point */
328 double *bvec=eimpli+3;
329
330 s6diff(epnt,eimpli,kdim,sdum1);
331 s6crss(epnt+ksizem3,sdum1,sdum2);
332 (void)s6norm(sdum2,kdim,sdum1,&kstat);
333 /* OK if sdum2 is zero -- tcurdst will be zero as well */
334
335 /* bvec = eimpli[3,4,5] is assumed to, be normalized */
336
337 t1 = s6scpr(sdum1,bvec,kdim);
338
339 /* t1 now contains the cosine of the angle between
340 (the cross product of the normal vector and the direction of the
341 point epnt relative to the point Q (in eimpli))
342 and the direction vector B. This is Equal to sin of PI/2 minus the angle
343 between the vectors, thus the actual angle can be calulated */
344
345 tcurdst = asin(t1);
346 tcurdst = fabs(tcurdst);
347 }
348
349 *jstat = 0;
350 goto out;
351
352 /* Projection not possible */
353 war02: *jstat = 2;
354 goto out;
355
356 /* Error in lower level routine. */
357
358 error : *jstat = kstat;
359 s6err("s1309",*jstat,kpos);
360 goto out;
361
362 out:
363 return(tcurdst);
364 }
365