ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/common/fvect.c
Revision: 2.9
Committed: Wed Mar 23 01:54:55 2005 UTC (19 years ago) by greg
Content type: text/plain
Branch: MAIN
CVS Tags: rad3R7P2, rad3R7P1, rad3R8, rad3R9
Changes since 2.8: +2 -2 lines
Log Message:
Added epsilon back to normalize() test

File Contents

# Content
1 #ifndef lint
2 static const char RCSid[] = "$Id: fvect.c,v 2.8 2003/09/16 06:30:20 greg Exp $";
3 #endif
4 /*
5 * fvect.c - routines for floating-point vector calculations
6 */
7
8 #include "copyright.h"
9
10 #include <math.h>
11 #include "fvect.h"
12
13
14 double
15 fdot( /* return the dot product of two vectors */
16 register FVECT v1,
17 register FVECT v2
18 )
19 {
20 return(DOT(v1,v2));
21 }
22
23
24 double
25 dist2( /* return square of distance between points */
26 register FVECT p1,
27 register FVECT p2
28 )
29 {
30 FVECT delta;
31
32 delta[0] = p2[0] - p1[0];
33 delta[1] = p2[1] - p1[1];
34 delta[2] = p2[2] - p1[2];
35
36 return(DOT(delta, delta));
37 }
38
39
40 double
41 dist2line( /* return square of distance to line */
42 FVECT p, /* the point */
43 FVECT ep1,
44 FVECT ep2 /* points on the line */
45 )
46 {
47 register double d, d1, d2;
48
49 d = dist2(ep1, ep2);
50 d1 = dist2(ep1, p);
51 d2 = d + d1 - dist2(ep2, p);
52
53 return(d1 - 0.25*d2*d2/d);
54 }
55
56
57 double
58 dist2lseg( /* return square of distance to line segment */
59 FVECT p, /* the point */
60 FVECT ep1,
61 FVECT ep2 /* the end points */
62 )
63 {
64 register double d, d1, d2;
65
66 d = dist2(ep1, ep2);
67 d1 = dist2(ep1, p);
68 d2 = dist2(ep2, p);
69
70 if (d2 > d1) { /* check if past endpoints */
71 if (d2 - d1 > d)
72 return(d1);
73 } else {
74 if (d1 - d2 > d)
75 return(d2);
76 }
77 d2 = d + d1 - d2;
78
79 return(d1 - 0.25*d2*d2/d); /* distance to line */
80 }
81
82
83 void
84 fcross( /* vres = v1 X v2 */
85 register FVECT vres,
86 register FVECT v1,
87 register FVECT v2
88 )
89 {
90 vres[0] = v1[1]*v2[2] - v1[2]*v2[1];
91 vres[1] = v1[2]*v2[0] - v1[0]*v2[2];
92 vres[2] = v1[0]*v2[1] - v1[1]*v2[0];
93 }
94
95
96 void
97 fvsum( /* vres = v0 + f*v1 */
98 register FVECT vres,
99 register FVECT v0,
100 register FVECT v1,
101 register double f
102 )
103 {
104 vres[0] = v0[0] + f*v1[0];
105 vres[1] = v0[1] + f*v1[1];
106 vres[2] = v0[2] + f*v1[2];
107 }
108
109
110 double
111 normalize( /* normalize a vector, return old magnitude */
112 register FVECT v
113 )
114 {
115 register double len, d;
116
117 d = DOT(v, v);
118
119 if (d <= FTINY*FTINY)
120 return(0.0);
121
122 if (d <= 1.0+FTINY && d >= 1.0-FTINY)
123 len = 0.5 + 0.5*d; /* first order approximation */
124 else
125 len = sqrt(d);
126
127 v[0] *= d = 1.0/len;
128 v[1] *= d;
129 v[2] *= d;
130
131 return(len);
132 }
133
134
135 int
136 closestapproach( /* closest approach of two rays */
137 RREAL t[2], /* returned distances along each ray */
138 FVECT rorg0, /* first origin */
139 FVECT rdir0, /* first direction (normalized) */
140 FVECT rorg1, /* second origin */
141 FVECT rdir1 /* second direction (normalized) */
142 )
143 {
144 double dotprod = DOT(rdir0, rdir1);
145 double denom = 1. - dotprod*dotprod;
146 double o1o2_d1;
147 FVECT o0o1;
148
149 if (denom <= FTINY) { /* check if lines are parallel */
150 t[0] = t[1] = 0.0;
151 return(0);
152 }
153 VSUB(o0o1, rorg0, rorg1);
154 o1o2_d1 = DOT(o0o1, rdir1);
155 t[0] = (o1o2_d1*dotprod - DOT(o0o1,rdir0)) / denom;
156 t[1] = o1o2_d1 + t[0]*dotprod;
157 return(1);
158 }
159
160
161 #if 0
162 int
163 closestapproach( /* closest approach of two rays */
164 RREAL t[2], /* returned distances along each ray */
165 FVECT rorg0, /* first origin */
166 FVECT rdir0, /* first direction (unnormalized) */
167 FVECT rorg1, /* second origin */
168 FVECT rdir1 /* second direction (unnormalized) */
169 )
170 {
171 double dotprod = DOT(rdir0, rdir1);
172 double d0n2 = DOT(rdir0, rdir0);
173 double d1n2 = DOT(rdir1, rdir1);
174 double denom = d0n2*d1n2 - dotprod*dotprod;
175 double o1o2_d1;
176 FVECT o0o1;
177
178 if (denom <= FTINY) { /* check if lines are parallel */
179 t[0] = t[1] = 0.0;
180 return(0);
181 }
182 VSUB(o0o1, rorg0, rorg1);
183 o1o2_d1 = DOT(o0o1, rdir1);
184 t[0] = (o1o2_d1*dotprod - DOT(o0o1,rdir0)*d1n2) / denom;
185 t[1] = (o1o2_d1 + t[0]*dotprod) / d1n2;
186 return(1);
187 }
188 #endif
189
190
191 void
192 spinvector( /* rotate vector around normal */
193 FVECT vres, /* returned vector */
194 FVECT vorig, /* original vector */
195 FVECT vnorm, /* normalized vector for rotation */
196 double theta /* left-hand radians */
197 )
198 {
199 double sint, cost, normprod;
200 FVECT vperp;
201 register int i;
202
203 if (theta == 0.0) {
204 if (vres != vorig)
205 VCOPY(vres, vorig);
206 return;
207 }
208 cost = cos(theta);
209 sint = sin(theta);
210 normprod = DOT(vorig, vnorm)*(1.-cost);
211 fcross(vperp, vnorm, vorig);
212 for (i = 0; i < 3; i++)
213 vres[i] = vorig[i]*cost + vnorm[i]*normprod + vperp[i]*sint;
214 }