1 /*
2  * transform.c
3  *
4  * Copyright (C) 1989, 1991, Craig E. Kolb
5  * All rights reserved.
6  *
7  * This software may be freely copied, modified, and redistributed
8  * provided that this copyright notice is preserved on all copies.
9  *
10  * You may not distribute this software, in whole or in part, as part of
11  * any commercial product without the express consent of the authors.
12  *
13  * There is no warranty or other guarantee of fitness of this software
14  * for any purpose.  It is provided solely "as is".
15  *
16  * $Id: transform.c,v 4.0.1.1 91/09/29 15:37:06 cek Exp Locker: cek $
17  *
18  * $Log:	transform.c,v $
19  * Revision 4.0.1.1  91/09/29  15:37:06  cek
20  * patch1: CoordSysTransform did not detect up-Z == -1.
21  *
22  * Revision 4.0  91/07/17  14:32:25  kolb
23  * Initial version.
24  *
25  */
26 #include "common.h"
27 
28 /*
29  * Matrices are indexed row-first; that is:
30  * matrix[ROW][COLUMN]
31  */
32 /*
33  * Allocate new structure that holds both object-to-world and
34  * world-to-object space transformation structures.  It probably
35  * should hold pointers to these structures.
36  */
37 Trans *
TransCreate(tr,meth)38 TransCreate(tr, meth)
39 TransRef tr;
40 TransMethods *meth;
41 {
42 	Trans *res;
43 
44 	res = (Trans *)share_malloc(sizeof(Trans));
45 	res->tr = tr;
46 	res->methods = meth;
47 	res->animated = FALSE;
48 	res->assoc = (ExprAssoc *)NULL;
49 	res->prev = res->next = (Trans *)NULL;
50 	MatrixInit(&res->trans);
51 	MatrixInit(&res->itrans);
52 	return res;
53 }
54 
55 void
TransFree(trans)56 TransFree(trans)
57 Trans *trans;
58 {
59 	if (trans->tr)
60 		free((voidstar)trans->tr);
61 	free((voidstar)trans);
62 }
63 
64 void
TransAssoc(trans,ptr,expr)65 TransAssoc(trans, ptr, expr)
66 Trans *trans;
67 Float *ptr;
68 Expr *expr;
69 {
70 	ExprAssoc *assoc;
71 
72 	if (expr->timevary) {
73 		/*
74 		 * Gotta store the sucker.
75 		 */
76 		trans->assoc = AssocCreate(ptr, expr, trans->assoc);
77 		trans->animated = TRUE;
78 	} else {
79 		*ptr = expr->value;
80 	}
81 	fflush(stderr);
82 }
83 
84 /*
85  * Allocate new transformation 'matrix'.
86  */
87 RSMatrix *
MatrixCreate()88 MatrixCreate()
89 {
90 	RSMatrix *res;
91 
92 	res = (RSMatrix *)share_malloc(sizeof(RSMatrix));
93 	MatrixInit(res);
94 	return res;
95 }
96 
97 /*
98  * Multiply m1 and m2, copy result into "res".
99  */
100 void
MatrixMult(t1,t2,res)101 MatrixMult(t1, t2, res)
102 RSMatrix *t1, *t2, *res;
103 {
104 	register int i;
105 	RSMatrix tmp;
106 
107 	for (i = 0; i < 3; i++) {
108 		tmp.matrix[i][0] = t1->matrix[i][0] * t2->matrix[0][0] +
109 			  	   t1->matrix[i][1] * t2->matrix[1][0] +
110 			  	   t1->matrix[i][2] * t2->matrix[2][0];
111 		tmp.matrix[i][1] = t1->matrix[i][0] * t2->matrix[0][1] +
112 			  	   t1->matrix[i][1] * t2->matrix[1][1] +
113 			  	   t1->matrix[i][2] * t2->matrix[2][1];
114 		tmp.matrix[i][2] = t1->matrix[i][0] * t2->matrix[0][2] +
115 			  	   t1->matrix[i][1] * t2->matrix[1][2] +
116 		  		   t1->matrix[i][2] * t2->matrix[2][2];
117 	}
118 
119 	tmp.translate.x = t1->translate.x * t2->matrix[0][0] +
120 			  t1->translate.y * t2->matrix[1][0] +
121 			  t1->translate.z * t2->matrix[2][0] + t2->translate.x;
122 	tmp.translate.y = t1->translate.x * t2->matrix[0][1] +
123 			  t1->translate.y * t2->matrix[1][1] +
124 			  t1->translate.z * t2->matrix[2][1] + t2->translate.y;
125 	tmp.translate.z = t1->translate.x * t2->matrix[0][2] +
126 			  t1->translate.y * t2->matrix[1][2] +
127 			  t1->translate.z * t2->matrix[2][2] + t2->translate.z;
128 	MatrixCopy(&tmp, res);
129 }
130 
131 /*
132  * Return transformation information to map the "coordinate system"
133  * with the given origin, "up" vector, radius, and up axis lengths to
134  * one in which the "up" vector is the Z axis and the x/y/up axes
135  * have unit length.  This is useful for transforming a general
136  * form of a primitive into a canonical, Z-axis aligned, unit size
137  * primitive, facilitating intersection testing.
138  * Assumes that "up" is normalized.
139  */
140 void
CoordSysTransform(origin,up,r,len,trans)141 CoordSysTransform(origin, up, r, len, trans)
142 Vector *origin, *up;
143 Float r, len;
144 Trans *trans;
145 {
146 	RSMatrix tmp;
147 	Vector atmp;
148 
149 	ScaleMatrix(r, r, len, &trans->trans);
150 	if (1. - fabs(up->z) < EPSILON) {
151 		atmp.x = 1.;
152 		atmp.y = atmp.z = 0.;
153 	} else {
154 		atmp.x = up->y;
155 		atmp.y = -up->x;
156 		atmp.z= 0.;
157 	}
158 	/*
159 	 * Might want to make sure that |up->z| is < 1.
160 	 */
161 	RotationMatrix(atmp.x, atmp.y, atmp.z, -acos(up->z), &tmp);
162 	MatrixMult(&trans->trans, &tmp, &trans->trans);
163 	TranslationMatrix(origin->x, origin->y, origin->z, &tmp);
164 	MatrixMult(&trans->trans, &tmp, &trans->trans);
165 	MatrixInvert(&trans->trans, &trans->itrans);
166 }
167 
168 void
TransCopy(from,into)169 TransCopy(from, into)
170 Trans *into, *from;
171 {
172 	MatrixCopy(&from->trans, &into->trans);
173 	MatrixCopy(&from->itrans, &into->itrans);
174 }
175 
176 void
TransInvert(from,into)177 TransInvert(from, into)
178 Trans *into, *from;
179 {
180 	RSMatrix ttmp;
181 	/*
182 	 * In case into == from...
183 	 */
184 	if (from == into) {
185 		ttmp = from->trans;
186 		into->trans = from->itrans;
187 		into->itrans = ttmp;
188 	} else {
189 		into->trans = from->itrans;
190 		into->itrans = from->trans;
191 	}
192 }
193 
194 /*
195  * Copy a given transformation structure.
196  */
197 void
MatrixCopy(from,into)198 MatrixCopy(from, into)
199 RSMatrix *into, *from;
200 {
201 	into->matrix[0][0] = from->matrix[0][0];
202 	into->matrix[0][1] = from->matrix[0][1];
203 	into->matrix[0][2] = from->matrix[0][2];
204 	into->matrix[1][0] = from->matrix[1][0];
205 	into->matrix[1][1] = from->matrix[1][1];
206 	into->matrix[1][2] = from->matrix[1][2];
207 	into->matrix[2][0] = from->matrix[2][0];
208 	into->matrix[2][1] = from->matrix[2][1];
209 	into->matrix[2][2] = from->matrix[2][2];
210 	into->translate = from->translate;
211 }
212 
213 void
TransInit(trans)214 TransInit(trans)
215 Trans *trans;
216 {
217 	MatrixInit(&trans->trans);
218 	MatrixInit(&trans->itrans);
219 }
220 
221 void
TransCompose(t1,t2,res)222 TransCompose(t1, t2, res)
223 Trans *t1, *t2, *res;
224 {
225 	MatrixMult(&t1->trans, &t2->trans, &res->trans);
226 	MatrixMult(&t2->itrans, &t1->itrans, &res->itrans);
227 }
228 
229 /*
230  * Initialize transformation structure.
231  */
232 void
MatrixInit(trans)233 MatrixInit(trans)
234 RSMatrix *trans;
235 {
236 	trans->matrix[0][0] = trans->matrix[1][1] = trans->matrix[2][2] = 1.;
237 	trans->matrix[0][1] = trans->matrix[0][2] = trans->matrix[1][0] =
238 	trans->matrix[1][2] = trans->matrix[2][0] = trans->matrix[2][1] = 0.;
239 	trans->translate.x = trans->translate.y = trans->translate.z = 0.;
240 }
241 
242 /*
243  * Calculate inverse of the given transformation structure.
244  */
245 void
MatrixInvert(trans,inverse)246 MatrixInvert(trans, inverse)
247 RSMatrix *inverse, *trans;
248 {
249 	RSMatrix ttmp;
250 	int i;
251 	Float d;
252 	extern int yylineno;
253 
254 	ttmp.matrix[0][0] = trans->matrix[1][1]*trans->matrix[2][2] -
255 			    trans->matrix[1][2]*trans->matrix[2][1];
256 	ttmp.matrix[1][0] = trans->matrix[1][0]*trans->matrix[2][2] -
257 			    trans->matrix[1][2]*trans->matrix[2][0];
258 	ttmp.matrix[2][0] = trans->matrix[1][0]*trans->matrix[2][1] -
259 			    trans->matrix[1][1]*trans->matrix[2][0];
260 
261 	ttmp.matrix[0][1] = trans->matrix[0][1]*trans->matrix[2][2] -
262 			    trans->matrix[0][2]*trans->matrix[2][1];
263 	ttmp.matrix[1][1] = trans->matrix[0][0]*trans->matrix[2][2] -
264 			    trans->matrix[0][2]*trans->matrix[2][0];
265 	ttmp.matrix[2][1] = trans->matrix[0][0]*trans->matrix[2][1] -
266 			    trans->matrix[0][1]*trans->matrix[2][0];
267 
268 	ttmp.matrix[0][2] = trans->matrix[0][1]*trans->matrix[1][2] -
269 			    trans->matrix[0][2]*trans->matrix[1][1];
270 	ttmp.matrix[1][2] = trans->matrix[0][0]*trans->matrix[1][2] -
271 			    trans->matrix[0][2]*trans->matrix[1][0];
272 	ttmp.matrix[2][2] = trans->matrix[0][0]*trans->matrix[1][1] -
273 			    trans->matrix[0][1]*trans->matrix[1][0];
274 
275 	d = trans->matrix[0][0]*ttmp.matrix[0][0] -
276 	    trans->matrix[0][1]*ttmp.matrix[1][0] +
277 	    trans->matrix[0][2]*ttmp.matrix[2][0];
278 
279 	if (fabs(d) < EPSILON*EPSILON)
280 		RLerror(RL_PANIC, "Singular matrix.\n",yylineno);
281 
282 	ttmp.matrix[0][0] /= d;
283 	ttmp.matrix[0][2] /= d;
284 	ttmp.matrix[1][1] /= d;
285 	ttmp.matrix[2][0] /= d;
286 	ttmp.matrix[2][2] /= d;
287 
288 	d = -d;
289 
290 	ttmp.matrix[0][1] /= d;
291 	ttmp.matrix[1][0] /= d;
292 	ttmp.matrix[1][2] /= d;
293 	ttmp.matrix[2][1] /= d;
294 
295 	ttmp.translate.x = -(ttmp.matrix[0][0]*trans->translate.x +
296 			     ttmp.matrix[1][0]*trans->translate.y +
297 			     ttmp.matrix[2][0]*trans->translate.z);
298 	ttmp.translate.y = -(ttmp.matrix[0][1]*trans->translate.x +
299 			     ttmp.matrix[1][1]*trans->translate.y +
300 			     ttmp.matrix[2][1]*trans->translate.z);
301 	ttmp.translate.z = -(ttmp.matrix[0][2]*trans->translate.x +
302 			     ttmp.matrix[1][2]*trans->translate.y +
303 			     ttmp.matrix[2][2]*trans->translate.z);
304 
305 	MatrixCopy(&ttmp, inverse);
306 }
307 
308 /*
309  * Apply a transformation to a point (translation affects the point).
310  */
311 void
PointTransform(vec,trans)312 PointTransform(vec, trans)
313 Vector *vec;
314 RSMatrix *trans;
315 {
316 	Vector tmp;
317 
318 	tmp.x = vec->x * trans->matrix[0][0] + vec->y * trans->matrix[1][0] +
319 			vec->z * trans->matrix[2][0] + trans->translate.x;
320 	tmp.y = vec->x * trans->matrix[0][1] + vec->y * trans->matrix[1][1] +
321 			vec->z * trans->matrix[2][1] + trans->translate.y;
322 	tmp.z = vec->x * trans->matrix[0][2] + vec->y * trans->matrix[1][2] +
323 			vec->z * trans->matrix[2][2] + trans->translate.z;
324 	*vec = tmp;
325 }
326 
327 /*
328  * 'c1x' is the X (0th) component of the first column, and so on.
329  */
330 void
ArbitraryMatrix(c1x,c2x,c3x,c1y,c2y,c3y,c1z,c2z,c3z,tx,ty,tz,trans)331 ArbitraryMatrix(c1x, c2x, c3x, c1y, c2y, c3y, c1z, c2z, c3z, tx, ty, tz, trans)
332 Float c1x, c1y, c1z, c2x, c2y, c2z, c3x, c3y, c3z, tx, ty, tz;
333 RSMatrix *trans;
334 {
335 	trans->matrix[0][0] = c1x;
336 	trans->matrix[1][0] = c1y;
337 	trans->matrix[2][0] = c1z;
338 
339 	trans->matrix[0][1] = c2x;
340 	trans->matrix[1][1] = c2y;
341 	trans->matrix[2][1] = c2z;
342 
343 	trans->matrix[0][2] = c3x;
344 	trans->matrix[1][2] = c3y;
345 	trans->matrix[2][2] = c3z;
346 
347 	trans->translate.x = tx;
348 	trans->translate.y = ty;
349 	trans->translate.z = tz;
350 }
351 
352 /*
353  * Apply transformation to a vector (translations have no effect).
354  */
355 void
VecTransform(vec,trans)356 VecTransform(vec, trans)
357 Vector *vec;
358 RSMatrix *trans;
359 {
360 	Vector tmp;
361 
362 	tmp.x = vec->x*trans->matrix[0][0] +
363 		vec->y*trans->matrix[1][0] + vec->z*trans->matrix[2][0];
364 	tmp.y = vec->x*trans->matrix[0][1] +
365 		vec->y*trans->matrix[1][1] + vec->z*trans->matrix[2][1];
366 	tmp.z = vec->x*trans->matrix[0][2] +
367 		vec->y*trans->matrix[1][2] + vec->z*trans->matrix[2][2];
368 
369 	*vec = tmp;
370 }
371 
372 /*
373  * Transform normal -- multiply by the transpose of the given
374  * matrix (which is the inverse of the 'desired' transformation).
375  */
376 void
NormalTransform(norm,it)377 NormalTransform(norm, it)
378 Vector *norm;
379 RSMatrix *it;
380 {
381 	Vector onorm;
382 
383 	onorm = *norm;
384 
385 	norm->x = onorm.x*it->matrix[0][0] + onorm.y*it->matrix[0][1] +
386 				onorm.z*it->matrix[0][2];
387 	norm->y = onorm.x*it->matrix[1][0] + onorm.y*it->matrix[1][1] +
388 				onorm.z*it->matrix[1][2];
389 	norm->z = onorm.x*it->matrix[2][0] + onorm.y*it->matrix[2][1] +
390 				onorm.z*it->matrix[2][2];
391 	(void)VecNormalize(norm);
392 }
393 
394 /*
395  * Transform "ray" by transforming the origin point and direction vector.
396  */
397 Float
RayTransform(ray,trans)398 RayTransform(ray, trans)
399 Ray *ray;
400 RSMatrix *trans;
401 {
402 	PointTransform(&ray->pos, trans);
403 	VecTransform(&ray->dir, trans);
404 	return VecNormalize(&ray->dir);
405 }
406 
407 void
TransPropagate(trans)408 TransPropagate(trans)
409 Trans *trans;
410 {
411 	(*trans->methods->propagate)(trans->tr, &trans->trans, &trans->itrans);
412 }
413 
414 void
TransResolveAssoc(trans)415 TransResolveAssoc(trans)
416 Trans *trans;
417 {
418 	Trans *curtrans;
419 	ExprAssoc *curassoc;
420 
421 	for (curtrans = trans; curtrans; curtrans = curtrans->next) {
422 		for (curassoc = curtrans->assoc; curassoc; curassoc = curassoc->next) {
423 			*curassoc->lhs = ExprEval(curassoc->expr);
424 		}
425 		if (curtrans->assoc)
426 			TransPropagate(curtrans);
427 	}
428 }
429 
430 void
TransComposeList(list,result)431 TransComposeList(list, result)
432 Trans *list, *result;
433 {
434 	TransCopy(list, result);
435 	for (list = list->next; list; list = list->next)
436 		TransCompose(list, result, result);
437 }
438