File: common\3d.c
1 /*
2 A word about the 3D library. Even though this library supports
3 three dimensions, the matrices are 4x4 for the following reason.
4 With normal 3 dimensional vectors, translation is an ADDITION,
5 and rotation is a MULTIPLICATION. A vector {x,y,z} is represented
6 as a 4-tuple {x,y,z,1}. It is then possible to define a 4x4
7 matrix such that multiplying the vector by the matrix translates
8 the vector. This allows combinations of translation and rotation
9 to be obtained in a single matrix by multiplying a translation
10 matrix and a rotation matrix together. Note that in the code,
11 vectors have three components; since the fourth component is
12 always 1, that value is not included in the vector variable to
13 save space, but the routines make use of the fourth component
14 (see vec_mult()). Similarly, the fourth column of EVERY matrix is
15 always
16 0
17 0
18 0
19 1
20 but currently the C version of a matrix includes this even though
21 it could be left out of the data structure and assumed in the
22 routines. Vectors are ROW vectors, and are always multiplied with
23 matrices FROM THE LEFT (e.g. vector*matrix). Also note the order
24 of indices of a matrix is matrix[row][column], and in usual C
25 fashion, numbering starts with 0.
26
27 TRANSLATION MATRIX = 1 0 0 0
28 0 1 0 0
29 0 0 1 0
30 Tx Ty Tz 1
31
32 SCALE MATRIX = Sx 0 0 0
33 0 Sy 0 0
34 0 0 Sz 0
35 0 0 0 1
36
37 Rotation about x axis i degrees:
38 ROTX(i) = 1 0 0 0
39 0 cosi sini 0
40 0 -sini cosi 0
41 0 0 0 1
42
43 Rotation about y axis i degrees:
44 ROTY(i) = cosi 0 -sini 0
45 0 1 0 0
46 sini 0 cosi 0
47 0 0 0 1
48
49 Rotation about z axis i degrees:
50 ROTZ(i) = cosi sini 0 0
51 -sini cosi 0 0
52 0 0 1 0
53 0 0 0 1
54
55 -- Tim Wegner April 22, 1989
56 */
57
58
59 #include <string.h>
60 /* see Fractint.c for a description of the "include" hierarchy */
61 #include "port.h"
62 #include "prototyp.h"
63
64 /* initialize a matrix and set to identity matrix
65 (all 0's, 1's on diagonal) */
66 void identity(MATRIX m)
67 {
68 int i,j;
69 for(i=0;i<CMAX;i++)
70 for(j=0;j<RMAX;j++)
71 if(i==j)
72 m[j][i] = 1.0;
73 else
74 m[j][i] = 0.0;
75 }
76
77 /* Multiply two matrices */
78 void mat_mul(MATRIX mat1, MATRIX mat2, MATRIX mat3)
79 {
80 /* result stored in MATRIX new to avoid problems
81 in case parameter mat3 == mat2 or mat 1 */
82 MATRIX new;
83 int i,j;
84 for(i=0;i<4;i++)
85 for(j=0;j<4;j++)
86 new[j][i] = mat1[j][0]*mat2[0][i]+
87 mat1[j][1]*mat2[1][i]+
88 mat1[j][2]*mat2[2][i]+
89 mat1[j][3]*mat2[3][i];
90 memcpy(mat3,new,sizeof(new));
91 }
92
93 /* multiply a matrix by a scalar */
94 void scale (double sx, double sy, double sz, MATRIX m)
95 {
96 MATRIX scale;
97 identity(scale);
98 scale[0][0] = sx;
99 scale[1][1] = sy;
100 scale[2][2] = sz;
101 mat_mul(m,scale,m);
102 }
103
104 /* rotate about X axis */
105 void xrot (double theta, MATRIX m)
106 {
107 MATRIX rot;
108 double sintheta,costheta;
109 sintheta = sin(theta);
110 costheta = cos(theta);
111 identity(rot);
112 rot[1][1] = costheta;
113 rot[1][2] = -sintheta;
114 rot[2][1] = sintheta;
115 rot[2][2] = costheta;
116 mat_mul(m,rot,m);
117 }
118
119 /* rotate about Y axis */
120 void yrot (double theta, MATRIX m)
121 {
122 MATRIX rot;
123 double sintheta,costheta;
124 sintheta = sin(theta);
125 costheta = cos(theta);
126 identity(rot);
127 rot[0][0] = costheta;
128 rot[0][2] = sintheta;
129 rot[2][0] = -sintheta;
130 rot[2][2] = costheta;
131 mat_mul(m,rot,m);
132 }
133
134 /* rotate about Z axis */
135 void zrot (double theta, MATRIX m)
136 {
137 MATRIX rot;
138 double sintheta,costheta;
139 sintheta = sin(theta);
140 costheta = cos(theta);
141 identity(rot);
142 rot[0][0] = costheta;
143 rot[0][1] = -sintheta;
144 rot[1][0] = sintheta;
145 rot[1][1] = costheta;
146 mat_mul(m,rot,m);
147 }
148
149 /* translate */
150 void trans (double tx, double ty, double tz, MATRIX m)
151 {
152 MATRIX trans;
153 identity(trans);
154 trans[3][0] = tx;
155 trans[3][1] = ty;
156 trans[3][2] = tz;
157 mat_mul(m,trans,m);
158 }
159
160 /* cross product - useful because cross is perpendicular to v and w */
161 int cross_product (VECTOR v, VECTOR w, VECTOR cross)
162 {
163 VECTOR tmp;
164 tmp[0] = v[1]*w[2] - w[1]*v[2];
165 tmp[1] = w[0]*v[2] - v[0]*w[2];
166 tmp[2] = v[0]*w[1] - w[0]*v[1];
167 cross[0] = tmp[0];
168 cross[1] = tmp[1];
169 cross[2] = tmp[2];
170 return(0);
171 }
172
173 /* cross product integer arguments (not fudged) */
174 /*** pb, unused
175 int icross_product (IVECTOR v, IVECTOR w, IVECTOR cross)
176 {
177 IVECTOR tmp;
178 tmp[0] = v[1]*w[2] - w[1]*v[2];
179 tmp[1] = w[0]*v[2] - v[0]*w[2];
180 tmp[2] = v[0]*w[1] - w[0]*v[1];
181 cross[0] = tmp[0];
182 cross[1] = tmp[1];
183 cross[2] = tmp[2];
184 return(0);
185 }
186 ***/
187
188 /* normalize a vector to length 1 */
189 int
190 normalize_vector(VECTOR v)
191 {
192 double vlength;
193 vlength = dot_product(v,v);
194
195 /* bailout if zero vlength */
196 if(vlength < FLT_MIN || vlength > FLT_MAX)
197 return(-1);
198 vlength = sqrt(vlength);
199 if(vlength < FLT_MIN)
200 return(-1);
201
202 v[0] /= vlength;
203 v[1] /= vlength;
204 v[2] /= vlength;
205 return(0);
206 }
207
208 /* multiply source vector s by matrix m, result in target t */
209 /* used to apply transformations to a vector */
210 int vmult(VECTOR s, MATRIX m, VECTOR t)
211 {
212 VECTOR tmp;
213 int i,j;
214 for(j=0;j<CMAX-1;j++)
215 {
216 tmp[j] = 0.0;
217 for(i=0;i<RMAX-1;i++)
218 tmp[j] += s[i]*m[i][j];
219 /* vector is really four dimensional with last component always 1 */
220 tmp[j] += m[3][j];
221 }
222 /* set target = tmp. Necessary to use tmp in case source = target */
223 memcpy(t,tmp,sizeof(tmp));
224 return(0);
225 }
226
227 /* multiply vector s by matrix m, result in s */
228 /* use with a function pointer in line3d.c */
229 /* must coordinate calling conventions with */
230 /* mult_vec in general.asm */
231 void mult_vec(VECTOR s)
232 {
233 VECTOR tmp;
234 int i,j;
235 for(j=0;j<CMAX-1;j++)
236 {
237 tmp[j] = 0.0;
238 for(i=0;i<RMAX-1;i++)
239 tmp[j] += s[i]*m[i][j];
240 /* vector is really four dimensional with last component always 1 */
241 tmp[j] += m[3][j];
242 }
243 /* set target = tmp. Necessary to use tmp in case source = target */
244 memcpy(s,tmp,sizeof(tmp));
245 }
246
247 /* perspective projection of vector v with respect to viewpont vector view */
248 int
249 perspective(VECTOR v)
250 {
251 double denom;
252 denom = view[2] - v[2];
253
254 if(denom >= 0.0)
255 {
256 v[0] = bad_value; /* clipping will catch these values */
257 v[1] = bad_value; /* so they won't plot values BEHIND viewer */
258 v[2] = bad_value;
259 return(-1);
260 }
261 v[0] = (v[0]*view[2] - view[0]*v[2])/denom;
262 v[1] = (v[1]*view[2] - view[1]*v[2])/denom;
263
264 /* calculation of z if needed later */
265 /* v[2] = v[2]/denom;*/
266 return(0);
267 }
268
269 /* long version of vmult and perspective combined for speed */
270 int
271 longvmultpersp(LVECTOR s, LMATRIX m, LVECTOR t0, LVECTOR t, LVECTOR lview,
272 int bitshift)
273 {
274 /* s: source vector */
275 /* m: transformation matrix */
276 /* t0: after transformation, before persp */
277 /* t: target vector */
278 /* lview: perspective viewer coordinates */
279 /* bitshift: fixed point conversion bitshift */
280 LVECTOR tmp;
281 int i,j, k;
282 overflow = 0;
283 k = CMAX-1; /* shorten the math if non-perspective and non-illum */
284 if (lview[2] == 0 && t0[0] == 0) k--;
285
286 for(j=0;j<k;j++)
287 {
288 tmp[j] = 0;
289 for(i=0;i<RMAX-1;i++)
290 tmp[j] += multiply(s[i],m[i][j],bitshift);
291 /* vector is really four dimensional with last component always 1 */
292 tmp[j] += m[3][j];
293 }
294 if(t0[0]) /* first component of t0 used as flag */
295 {
296 /* faster than for loop, if less general */
297 t0[0] = tmp[0];
298 t0[1] = tmp[1];
299 t0[2] = tmp[2];
300 }
301 if (lview[2] != 0) /* perspective 3D */
302 {
303
304 LVECTOR tmpview;
305 long denom;
306
307 denom = lview[2] - tmp[2];
308 if (denom >= 0) /* bail out if point is "behind" us */
309 {
310 t[0] = bad_value;
311 t[0] = t[0]<<bitshift;
312 t[1] = t[0];
313 t[2] = t[0];
314 return(-1);
315 }
316
317 /* doing math in this order helps prevent overflow */
318 tmpview[0] = divide(lview[0],denom,bitshift);
319 tmpview[1] = divide(lview[1],denom,bitshift);
320 tmpview[2] = divide(lview[2],denom,bitshift);
321
322 tmp[0] = multiply(tmp[0], tmpview[2], bitshift) -
323 multiply(tmpview[0], tmp[2], bitshift);
324
325 tmp[1] = multiply(tmp[1], tmpview[2], bitshift) -
326 multiply(tmpview[1], tmp[2], bitshift);
327
328 /* z coordinate if needed */
329 /* tmp[2] = divide(lview[2],denom); */
330 }
331
332 /* set target = tmp. Necessary to use tmp in case source = target */
333 /* faster than for loop, if less general */
334 t[0] = tmp[0];
335 t[1] = tmp[1];
336 t[2] = tmp[2];
337 return(overflow);
338 }
339
340 /* Long version of perspective. Because of use of fixed point math, there
341 is danger of overflow and underflow */
342 int
343 longpersp(LVECTOR lv, LVECTOR lview, int bitshift)
344 {
345 LVECTOR tmpview;
346 long denom;
347 overflow = 0;
348 denom = lview[2] - lv[2];
349 if (denom >= 0) /* bail out if point is "behind" us */
350 {
351 lv[0] = bad_value;
352 lv[0] = lv[0]<<bitshift;
353 lv[1] = lv[0];
354 lv[2] = lv[0];
355 return(-1);
356 }
357
358 /* doing math in this order helps prevent overflow */
359 tmpview[0] = divide(lview[0],denom,bitshift);
360 tmpview[1] = divide(lview[1],denom,bitshift);
361 tmpview[2] = divide(lview[2],denom,bitshift);
362
363 lv[0] = multiply(lv[0], tmpview[2], bitshift) -
364 multiply(tmpview[0], lv[2], bitshift);
365
366 lv[1] = multiply(lv[1], tmpview[2], bitshift) -
367 multiply(tmpview[1], lv[2], bitshift);
368
369 /* z coordinate if needed */
370 /* lv[2] = divide(lview[2],denom); */
371 return(overflow);
372 }
373
374 int longvmult(LVECTOR s,LMATRIX m,LVECTOR t,int bitshift)
375 {
376 LVECTOR tmp;
377 int i,j, k;
378 overflow = 0;
379 k = CMAX-1;
380
381 for(j=0;j<k;j++)
382 {
383 tmp[j] = 0;
384 for(i=0;i<RMAX-1;i++)
385 tmp[j] += multiply(s[i],m[i][j],bitshift);
386 /* vector is really four dimensional with last component always 1 */
387 tmp[j] += m[3][j];
388 }
389
390 /* set target = tmp. Necessary to use tmp in case source = target */
391 /* faster than for loop, if less general */
392 t[0] = tmp[0];
393 t[1] = tmp[1];
394 t[2] = tmp[2];
395 return(overflow);
396 }
397