1 |
greg |
2.2 |
#ifndef lint |
2 |
|
|
static const char RCSid[] = "$Id$"; |
3 |
|
|
#endif |
4 |
greg |
2.1 |
#include <stdio.h> |
5 |
|
|
#include <stdlib.h> |
6 |
|
|
|
7 |
|
|
#include "g3affine.h" |
8 |
|
|
#include "gbasic.h" |
9 |
|
|
#include "g3sphere.h" |
10 |
|
|
|
11 |
|
|
static g3AVec g3av_zero = {0.0,0.0,0.0,0.0}; |
12 |
|
|
|
13 |
|
|
static g3ATrans g3at_unit = {{1.0,0.0,0.0,0.0}, |
14 |
|
|
{0.0,1.0,0.0,0.0}, |
15 |
|
|
{0.0,0.0,1.0,0.0}, |
16 |
|
|
{0.0,0.0,0.0,1.0}}; |
17 |
|
|
|
18 |
|
|
static g3ATrans g3at_zero = {{0.0,0.0,0.0,0.0}, |
19 |
|
|
{0.0,0.0,0.0,0.0}, |
20 |
|
|
{0.0,0.0,0.0,0.0}, |
21 |
|
|
{0.0,0.0,0.0,0.0}}; |
22 |
|
|
|
23 |
|
|
void g3at_tounit(g3ATrans t) |
24 |
|
|
{ |
25 |
|
|
g3at_copy(t,g3at_unit); |
26 |
|
|
} |
27 |
|
|
|
28 |
|
|
static int nr_gaussj(g3Float a[4][4],g3Float b[4]) |
29 |
|
|
{ |
30 |
|
|
int *indxc,*indxr,*ipiv; |
31 |
|
|
int i,icol,irow,j,k,l,ll; |
32 |
|
|
g3Float big,dum,pivinv; |
33 |
|
|
|
34 |
|
|
indxc = (int*) calloc(4,sizeof(int)); |
35 |
|
|
indxr = (int*) calloc(4,sizeof(int)); |
36 |
|
|
ipiv = (int*) calloc(4,sizeof(int)); |
37 |
|
|
|
38 |
|
|
icol = irow = 0; |
39 |
|
|
for (j=0;j<4;j++) |
40 |
|
|
ipiv[j]=0; |
41 |
|
|
for (i=0;i<4;i++) { |
42 |
|
|
big=0.0; |
43 |
|
|
for (j=0;j<4;j++) |
44 |
|
|
if (ipiv[j] != 1) |
45 |
|
|
for (k=0;k<4;k++) { |
46 |
|
|
if (ipiv[k] == 0) { |
47 |
|
|
if (fabs(a[j][k]) >= big) { |
48 |
|
|
big = fabs(a[j][k]); |
49 |
|
|
irow = j; |
50 |
|
|
icol = k; |
51 |
|
|
} |
52 |
|
|
} else if (ipiv[k] > 1) { |
53 |
|
|
return 0; |
54 |
|
|
} |
55 |
|
|
} |
56 |
|
|
++(ipiv[icol]); |
57 |
|
|
if (irow != icol) { |
58 |
|
|
for (l=0;l<4;l++) |
59 |
|
|
GB_SWAP(a[irow][l],a[icol][l]) |
60 |
|
|
GB_SWAP(b[irow],b[icol]) |
61 |
|
|
} |
62 |
|
|
indxr[i] = irow; |
63 |
|
|
indxc[i] = icol; |
64 |
|
|
if (a[icol][icol] == 0.0) |
65 |
|
|
return 0; |
66 |
|
|
pivinv = 1.0/a[icol][icol]; |
67 |
|
|
a[icol][icol] = 1.0; |
68 |
|
|
for (l=0;l<4;l++) |
69 |
|
|
a[icol][l] *= pivinv; |
70 |
|
|
b[icol] *= pivinv; |
71 |
|
|
for (ll=0;ll<4;ll++) { |
72 |
|
|
if (ll != icol) { |
73 |
|
|
dum = a[ll][icol]; |
74 |
|
|
a[ll][icol] = 0.0; |
75 |
|
|
for (l=0;l<4;l++) |
76 |
|
|
a[ll][l] -= a[icol][l]*dum; |
77 |
|
|
b[ll] -= b[icol]*dum; |
78 |
|
|
} |
79 |
|
|
} |
80 |
|
|
} |
81 |
|
|
for (l=3;l>=0;l--) { |
82 |
|
|
if (indxr[l] != indxc[l]) |
83 |
|
|
for (k=0;k<4;k++) |
84 |
|
|
GB_SWAP(a[k][indxr[l]],a[k][indxc[l]]); |
85 |
|
|
} |
86 |
|
|
free(ipiv); |
87 |
|
|
free(indxr); |
88 |
|
|
free(indxc); |
89 |
|
|
return 1; |
90 |
|
|
} |
91 |
|
|
|
92 |
|
|
void g3at_translate(g3ATrans t,g3Vec tv) |
93 |
|
|
{ |
94 |
|
|
int i; |
95 |
|
|
g3at_copy(t,g3at_unit); |
96 |
|
|
for(i=0;i<3;i++) |
97 |
|
|
t[i][3] += tv[i]; |
98 |
|
|
} |
99 |
|
|
|
100 |
|
|
void g3at_add_trans(g3ATrans t,g3Vec tv) |
101 |
|
|
{ |
102 |
|
|
g3ATrans tt; |
103 |
|
|
|
104 |
|
|
g3at_translate(tt,tv); |
105 |
|
|
g3at_comb_to(t,tt,t); |
106 |
|
|
} |
107 |
|
|
|
108 |
|
|
void g3at_prepend_trans(g3ATrans t,g3Vec tv) |
109 |
|
|
{ |
110 |
|
|
g3ATrans tt; |
111 |
|
|
|
112 |
|
|
g3at_translate(tt,tv); |
113 |
|
|
g3at_comb(t,tt); |
114 |
|
|
} |
115 |
|
|
|
116 |
|
|
void g3at_rotate(g3ATrans t,g3Vec rnorm,g3Float ang) |
117 |
|
|
{ |
118 |
|
|
g3Vec axe; |
119 |
|
|
int aid,i; |
120 |
|
|
|
121 |
|
|
axe = g3v_create(); |
122 |
|
|
g3at_copy(t,g3at_unit); |
123 |
|
|
|
124 |
|
|
for(aid=0;aid<3;aid++) { |
125 |
|
|
g3v_zero(axe); |
126 |
|
|
axe[aid] = 1.0; |
127 |
|
|
g3v_rotate(axe,axe,rnorm,ang); |
128 |
|
|
for(i=0;i<3;i++) |
129 |
|
|
t[i][aid] = axe[i]; |
130 |
|
|
} |
131 |
|
|
g3v_free(axe); |
132 |
|
|
} |
133 |
|
|
|
134 |
|
|
void g3at_add_rot(g3ATrans t,g3Vec rnorm,g3Float ang) |
135 |
|
|
{ |
136 |
|
|
g3ATrans rt; |
137 |
|
|
|
138 |
|
|
g3at_rotate(rt,rnorm,ang); |
139 |
|
|
g3at_comb_to(t,rt,t); |
140 |
|
|
} |
141 |
|
|
|
142 |
|
|
void g3at_prepend_rot(g3ATrans t,g3Vec rnorm,g3Float ang) |
143 |
|
|
{ |
144 |
|
|
g3ATrans rt; |
145 |
|
|
|
146 |
|
|
g3at_rotate(rt,rnorm,ang); |
147 |
|
|
g3at_comb(t,rt); |
148 |
|
|
} |
149 |
|
|
|
150 |
|
|
void g3at_btrans(g3ATrans t,g3Vec xv,g3Vec yv,g3Vec zv) |
151 |
|
|
{ |
152 |
|
|
int i; |
153 |
|
|
|
154 |
|
|
g3at_copy(t,g3at_unit); |
155 |
|
|
for(i=0;i<3;i++) { |
156 |
|
|
t[i][0] = xv[i]; |
157 |
|
|
t[i][1] = yv[i]; |
158 |
|
|
t[i][2] = zv[i]; |
159 |
|
|
} |
160 |
|
|
} |
161 |
|
|
|
162 |
|
|
int g3at_inverse(g3ATrans t) |
163 |
|
|
{ |
164 |
|
|
g3Float v[] = {0,0,0,0}; |
165 |
|
|
|
166 |
|
|
if (!nr_gaussj(t,v)) { |
167 |
|
|
fprintf(stderr,"g3at_inverse: singular matrix\n"); |
168 |
|
|
return 0; |
169 |
|
|
} |
170 |
|
|
return 1; |
171 |
|
|
} |
172 |
|
|
|
173 |
|
|
void g3at_apply_h(g3ATrans t,g3AVec v) |
174 |
|
|
{ |
175 |
|
|
int i,j; |
176 |
|
|
g3AVec vt; |
177 |
|
|
|
178 |
|
|
g3av_copy(vt,v); |
179 |
|
|
g3av_copy(v,g3av_zero); |
180 |
|
|
for(i=0;i<4;i++) |
181 |
|
|
for(j=0;j<4;j++) |
182 |
|
|
v[i] += t[i][j]*vt[j]; |
183 |
|
|
g3av_tonorm(v); |
184 |
|
|
} |
185 |
|
|
|
186 |
|
|
void g3av_print(g3AVec vh,FILE* outf) |
187 |
|
|
{ |
188 |
|
|
int i; |
189 |
|
|
|
190 |
|
|
for(i=0;i<4;i++) |
191 |
|
|
fprintf(outf,"%f ",vh[i]); |
192 |
|
|
} |
193 |
|
|
|
194 |
|
|
void g3at_apply(g3ATrans t,g3Vec v) |
195 |
|
|
{ |
196 |
|
|
g3AVec res; |
197 |
|
|
|
198 |
|
|
g3av_vectohom(res,v); |
199 |
|
|
g3at_apply_h(t,res); |
200 |
|
|
g3v_set(v,res[0],res[1],res[2]); |
201 |
|
|
} |
202 |
|
|
|
203 |
|
|
void g3at_iapply(g3Vec v,g3ATrans t) |
204 |
|
|
{ |
205 |
|
|
g3AVec res; |
206 |
|
|
|
207 |
|
|
g3av_copy(res,v); |
208 |
|
|
g3at_iapply_h(res,t); |
209 |
|
|
g3v_set(v,res[0],res[1],res[2]); |
210 |
|
|
} |
211 |
|
|
|
212 |
|
|
void g3at_iapply_h(g3Vec v,g3ATrans t) |
213 |
|
|
{ |
214 |
|
|
int i,j; |
215 |
|
|
g3AVec vt; |
216 |
|
|
|
217 |
|
|
g3av_copy(vt,v); |
218 |
|
|
g3av_copy(v,g3av_zero); |
219 |
|
|
for(i=0;i<4;i++) |
220 |
|
|
for(j=0;j<4;j++) |
221 |
|
|
v[i] += t[j][i]*vt[j]; |
222 |
|
|
g3av_tonorm(v); |
223 |
|
|
} |
224 |
|
|
|
225 |
|
|
void g3at_comb(g3ATrans t,g3ATrans tf) |
226 |
|
|
{ |
227 |
|
|
g3ATrans res; |
228 |
|
|
|
229 |
|
|
g3at_comb_to(res,t,tf); |
230 |
|
|
g3at_copy(t,res); |
231 |
|
|
} |
232 |
|
|
|
233 |
|
|
void g3at_comb_to(g3ATrans res,g3ATrans t,g3ATrans tf) |
234 |
|
|
{ |
235 |
|
|
int i,j,k; |
236 |
|
|
g3ATrans temp; |
237 |
|
|
g3at_copy(temp,g3at_zero); |
238 |
|
|
for(i=0;i<4;i++) |
239 |
|
|
for(j=0;j<4;j++) |
240 |
|
|
for(k=0;k<4;k++) |
241 |
|
|
temp[i][j] += t[i][k]*tf[k][j]; |
242 |
|
|
g3at_copy(res,temp); |
243 |
|
|
} |
244 |
|
|
|
245 |
|
|
|
246 |
|
|
void g3av_vectohom(g3AVec res,g3Vec v) |
247 |
|
|
{ |
248 |
|
|
int i; |
249 |
|
|
|
250 |
|
|
for(i=0;i<3;i++) |
251 |
|
|
res[i] = v[i]; |
252 |
|
|
res[3] = 1.0; |
253 |
|
|
} |
254 |
|
|
|
255 |
|
|
void g3av_homtovec(g3Vec res,g3AVec h) |
256 |
|
|
{ |
257 |
|
|
int i; |
258 |
|
|
|
259 |
|
|
g3av_tonorm(h); |
260 |
|
|
for(i=0;i<3;i++) |
261 |
|
|
res[i] = h[i]; |
262 |
|
|
} |
263 |
|
|
|
264 |
|
|
int g3av_tonorm(g3AVec v) |
265 |
|
|
{ |
266 |
|
|
int i; |
267 |
|
|
|
268 |
|
|
if (gb_epseq(v[3],0.0)) |
269 |
|
|
return 0; |
270 |
|
|
for(i=0;i<3;i++) |
271 |
|
|
v[i] /= v[3]; |
272 |
|
|
v[3] = 1.0; |
273 |
|
|
return 1; |
274 |
|
|
} |
275 |
|
|
|
276 |
|
|
void g3av_copy(g3AVec v,g3AVec vsrc) |
277 |
|
|
{ |
278 |
|
|
int i; |
279 |
|
|
|
280 |
|
|
for(i=0;i<4;i++) |
281 |
|
|
v[i] = vsrc[i]; |
282 |
|
|
} |
283 |
|
|
|
284 |
|
|
void g3at_print(g3ATrans t,FILE* outf) |
285 |
|
|
{ |
286 |
|
|
int i,j; |
287 |
|
|
|
288 |
|
|
for(i=0;i<4;i++) { |
289 |
|
|
for(j=0;j<4;j++) |
290 |
|
|
fprintf(outf,"%f ",t[i][j]); |
291 |
|
|
fprintf(outf,"\n"); |
292 |
|
|
} |
293 |
|
|
} |
294 |
|
|
|
295 |
|
|
void g3at_copy(g3ATrans t,g3ATrans tsrc) |
296 |
|
|
{ |
297 |
|
|
int i,j; |
298 |
|
|
|
299 |
|
|
for(i=0;i<4;i++) |
300 |
|
|
for(j=0;j<4;j++) |
301 |
|
|
t[i][j] = tsrc[i][j]; |
302 |
|
|
} |
303 |
|
|
|
304 |
|
|
int g3at_ctrans(g3ATrans t,g3Vec vdir,g3Vec vup) |
305 |
|
|
{ |
306 |
|
|
g3Vec x,y,z; |
307 |
|
|
|
308 |
|
|
x = g3v_create(); |
309 |
|
|
y = g3v_create(); |
310 |
|
|
z = g3v_create(); |
311 |
|
|
g3v_copy(z,vdir); |
312 |
|
|
g3v_copy(x,vup); |
313 |
|
|
g3v_normalize(z); |
314 |
|
|
g3v_normalize(x); |
315 |
|
|
g3v_cross(y,z,x); |
316 |
|
|
if (!g3v_normalize(y)) { |
317 |
|
|
fprintf(stderr,"g3at_sph_ctrans: parallel\n"); |
318 |
|
|
g3v_free(x); |
319 |
|
|
g3v_free(y); |
320 |
|
|
g3v_free(z); |
321 |
|
|
return 0; |
322 |
|
|
} |
323 |
|
|
g3v_cross(x,y,z); |
324 |
|
|
if (!g3v_normalize(x)) { |
325 |
|
|
fprintf(stderr,"g3at_sph_ctrans: parallel\n"); |
326 |
|
|
g3v_free(x); |
327 |
|
|
g3v_free(y); |
328 |
|
|
g3v_free(z); |
329 |
|
|
return 0; |
330 |
|
|
} |
331 |
|
|
g3at_btrans(t,x,y,z); |
332 |
|
|
g3v_free(x); |
333 |
|
|
g3v_free(y); |
334 |
|
|
g3v_free(z); |
335 |
|
|
return 1; |
336 |
|
|
} |
337 |
|
|
|
338 |
|
|
int g3at_sph_ctrans(g3ATrans t,g3Vec vdir,g3Vec vup) |
339 |
|
|
{ |
340 |
|
|
g3Vec x,z; |
341 |
|
|
int res; |
342 |
|
|
|
343 |
|
|
x = g3v_create(); |
344 |
|
|
z = g3v_create(); |
345 |
|
|
g3v_copy(z,vdir); |
346 |
|
|
g3v_copy(x,vup); |
347 |
|
|
z[G3S_RAD] = x[G3S_RAD] = 1.0; |
348 |
|
|
g3s_sphtocc(x,x); |
349 |
|
|
g3s_sphtocc(z,z); |
350 |
|
|
res = g3at_ctrans(t,z,x); |
351 |
|
|
g3v_free(x); |
352 |
|
|
g3v_free(z); |
353 |
|
|
|
354 |
|
|
return res; |
355 |
|
|
} |
356 |
|
|
|
357 |
|
|
#ifdef GL_CONVERSIONS |
358 |
|
|
void g3at_get_openGL(g3ATrans t,GLdouble* glm) |
359 |
|
|
{ |
360 |
|
|
int r,c; |
361 |
|
|
|
362 |
|
|
for(r=0;r<4;r++) |
363 |
|
|
for(c=0;c<4;c++) |
364 |
|
|
glm[c+(r*4)] = t[c][r]; |
365 |
|
|
} |
366 |
|
|
|
367 |
|
|
void g3at_from_openGL(g3ATrans t,GLdouble* glm) |
368 |
|
|
{ |
369 |
|
|
int r,c; |
370 |
|
|
|
371 |
|
|
for(r=0;r<4;r++) |
372 |
|
|
for(c=0;c<4;c++) |
373 |
|
|
t[c][r] = glm[c+(r*4)]; |
374 |
|
|
} |
375 |
|
|
#endif |
376 |
|
|
|
377 |
|
|
#ifdef G3AFFINETEST |
378 |
|
|
|
379 |
|
|
|
380 |
|
|
int main(int argc,char** argv) |
381 |
|
|
{ |
382 |
|
|
g3Vec v,iv,x,y,z; |
383 |
|
|
g3ATrans t,it; |
384 |
|
|
|
385 |
|
|
g3at_copy(t,g3at_unit); |
386 |
|
|
v = g3v_create(); |
387 |
|
|
iv = g3v_create(); |
388 |
|
|
x = g3v_create(); |
389 |
|
|
y = g3v_create(); |
390 |
|
|
z = g3v_create(); |
391 |
|
|
g3v_set(z,1,0,0); |
392 |
|
|
g3v_set(x,atof(argv[1]),atof(argv[2]),atof(argv[3])); |
393 |
|
|
//g3v_set(y,atof(argv[4]),atof(argv[5]),atof(argv[6])); |
394 |
|
|
//g3v_set(z,atof(argv[7]),atof(argv[8]),atof(argv[9])); |
395 |
|
|
g3at_tounit(t); |
396 |
|
|
g3at_rotate(t,z,DEG2RAD(-90)); |
397 |
|
|
g3v_set(z,0,0,1); |
398 |
|
|
g3at_rotate(it,z,DEG2RAD(180)); |
399 |
|
|
g3at_comb_to(t,it,t); |
400 |
|
|
g3at_apply(t,x); |
401 |
|
|
g3v_print(x,stdout); |
402 |
|
|
printf("\n"); |
403 |
|
|
//g3v_print(y,stdout); |
404 |
|
|
//printf("\n"); |
405 |
|
|
exit(19); |
406 |
|
|
g3v_set(z,1,DEG2RAD(atof(argv[1])),DEG2RAD(atof(argv[2]))); |
407 |
|
|
g3v_set(x,1,DEG2RAD(atof(argv[3])),DEG2RAD(atof(argv[4]))); |
408 |
|
|
g3at_sph_ctrans(t,z,x); |
409 |
|
|
g3at_copy(it,t); |
410 |
|
|
g3at_inverse(it); |
411 |
|
|
g3at_print(t,stderr); |
412 |
|
|
fprintf(stderr,"\n"); |
413 |
|
|
g3v_set(v,1,DEG2RAD(atof(argv[5])),DEG2RAD(atof(argv[6]))); |
414 |
|
|
g3s_sphtocc(v,v); |
415 |
|
|
g3v_copy(iv,v); |
416 |
|
|
g3at_apply(it,iv); |
417 |
|
|
g3at_apply(t,v); |
418 |
|
|
g3s_cctosph(v,v); |
419 |
|
|
g3s_cctosph(iv,iv); |
420 |
|
|
g3s_sphwrap(v); |
421 |
|
|
g3s_sphwrap(iv); |
422 |
|
|
fprintf(stderr,"tt %f %f\n",RAD2DEG(v[1]),RAD2DEG(v[2])); |
423 |
|
|
fprintf(stderr,"ise %f %f\n",RAD2DEG(iv[1]),RAD2DEG(iv[2])); |
424 |
|
|
return EXIT_SUCCESS; |
425 |
|
|
} |
426 |
|
|
|
427 |
|
|
#endif |