![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to 3d.c CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [Research Unix] / researchv9 / jerq / src / lib / jj|
Return to 3d.c CVS log | Up to [Research Unix] / researchv9 / jerq / src / lib / jj |
1.1 root 1: /*% 3cc -c %
2: * Fixed-point 3d coordinate transform routines for machines with
3: * 16 bit shorts and 32 bit longs.
4: */
5: #include <jerq.h>
6: #include <3d.h>
7: #include <font.h>
8: static Hcoordl cur;
9: static Point cen;
10: static int wid;
11: /*
12: * angles are scaled so that PI=>32768
13: * Since unsigned shorts are 16 bits, 2*PI overflows and yields
14: * 0, which seems appropriate.
15: */
16: angle cvtangle[361]={
17: 0, 182, 364, 546, 728, 910, 1092, 1274,
18: 1456, 1638, 1820, 2002, 2184, 2366, 2548, 2730,
19: 2912, 3094, 3276, 3458, 3640, 3822, 4004, 4187,
20: 4369, 4551, 4733, 4915, 5097, 5279, 5461, 5643,
21: 5825, 6007, 6189, 6371, 6553, 6735, 6917, 7099,
22: 7281, 7463, 7645, 7827, 8009, 8192, 8374, 8556,
23: 8738, 8920, 9102, 9284, 9466, 9648, 9830, 10012,
24: 10194, 10376, 10558, 10740, 10922, 11104, 11286, 11468,
25: 11650, 11832, 12014, 12196, 12379, 12561, 12743, 12925,
26: 13107, 13289, 13471, 13653, 13835, 14017, 14199, 14381,
27: 14563, 14745, 14927, 15109, 15291, 15473, 15655, 15837,
28: 16019, 16201, 16384, 16566, 16748, 16930, 17112, 17294,
29: 17476, 17658, 17840, 18022, 18204, 18386, 18568, 18750,
30: 18932, 19114, 19296, 19478, 19660, 19842, 20024, 20206,
31: 20388, 20571, 20753, 20935, 21117, 21299, 21481, 21663,
32: 21845, 22027, 22209, 22391, 22573, 22755, 22937, 23119,
33: 23301, 23483, 23665, 23847, 24029, 24211, 24393, 24576,
34: 24758, 24940, 25122, 25304, 25486, 25668, 25850, 26032,
35: 26214, 26396, 26578, 26760, 26942, 27124, 27306, 27488,
36: 27670, 27852, 28034, 28216, 28398, 28580, 28763, 28945,
37: 29127, 29309, 29491, 29673, 29855, 30037, 30219, 30401,
38: 30583, 30765, 30947, 31129, 31311, 31493, 31675, 31857,
39: 32039, 32221, 32403, 32585, 32768, 32950, 33132, 33314,
40: 33496, 33678, 33860, 34042, 34224, 34406, 34588, 34770,
41: 34952, 35134, 35316, 35498, 35680, 35862, 36044, 36226,
42: 36408, 36590, 36772, 36955, 37137, 37319, 37501, 37683,
43: 37865, 38047, 38229, 38411, 38593, 38775, 38957, 39139,
44: 39321, 39503, 39685, 39867, 40049, 40231, 40413, 40595,
45: 40777, 40960, 41142, 41324, 41506, 41688, 41870, 42052,
46: 42234, 42416, 42598, 42780, 42962, 43144, 43326, 43508,
47: 43690, 43872, 44054, 44236, 44418, 44600, 44782, 44964,
48: 45147, 45329, 45511, 45693, 45875, 46057, 46239, 46421,
49: 46603, 46785, 46967, 47149, 47331, 47513, 47695, 47877,
50: 48059, 48241, 48423, 48605, 48787, 48969, 49152, 49334,
51: 49516, 49698, 49880, 50062, 50244, 50426, 50608, 50790,
52: 50972, 51154, 51336, 51518, 51700, 51882, 52064, 52246,
53: 52428, 52610, 52792, 52974, 53156, 53339, 53521, 53703,
54: 53885, 54067, 54249, 54431, 54613, 54795, 54977, 55159,
55: 55341, 55523, 55705, 55887, 56069, 56251, 56433, 56615,
56: 56797, 56979, 57161, 57344, 57526, 57708, 57890, 58072,
57: 58254, 58436, 58618, 58800, 58982, 59164, 59346, 59528,
58: 59710, 59892, 60074, 60256, 60438, 60620, 60802, 60984,
59: 61166, 61348, 61531, 61713, 61895, 62077, 62259, 62441,
60: 62623, 62805, 62987, 63169, 63351, 63533, 63715, 63897,
61: 64079, 64261, 64443, 64625, 64807, 64989, 65171, 65353,
62: 65536,
63: };
64: static fract sintab[128+2]={
65: 0, 201, 402, 603, 803, 1004, 1205, 1405,
66: 1605, 1805, 2005, 2204, 2404, 2602, 2801, 2998,
67: 3196, 3393, 3589, 3785, 3980, 4175, 4369, 4563,
68: 4756, 4948, 5139, 5329, 5519, 5708, 5896, 6083,
69: 6269, 6455, 6639, 6822, 7005, 7186, 7366, 7545,
70: 7723, 7900, 8075, 8249, 8423, 8594, 8765, 8934,
71: 9102, 9268, 9434, 9597, 9759, 9920, 10079, 10237,
72: 10393, 10548, 10701, 10853, 11002, 11150, 11297, 11442,
73: 11585, 11726, 11866, 12003, 12139, 12273, 12406, 12536,
74: 12665, 12791, 12916, 13038, 13159, 13278, 13395, 13510,
75: 13622, 13733, 13842, 13948, 14053, 14155, 14255, 14353,
76: 14449, 14543, 14634, 14723, 14810, 14895, 14978, 15058,
77: 15136, 15212, 15286, 15357, 15426, 15492, 15557, 15618,
78: 15678, 15735, 15790, 15842, 15892, 15940, 15985, 16028,
79: 16069, 16107, 16142, 16175, 16206, 16234, 16260, 16284,
80: 16305, 16323, 16339, 16353, 16364, 16372, 16379, 16382,
81: 16384, 16382,
82: };
83: fract isin(x)
84: register angle x;
85: {
86: register fract a, b;
87: register i, f=0;
88: if(x>=PI){
89: x-=PI;
90: f=1;
91: }
92: if(x>PI/2)
93: x=PI-x;
94: i=x>>7;
95: a=sintab[i];
96: b=sintab[i+1]-a;
97: i=x&0177;
98: a+=(b*i+64)/128;
99: return(f?-a:a);
100: }
101: fract icos(x)
102: register angle x;
103: {
104: return(isin(x+PI/2));
105: }
106: long isqrt(a)
107: register long a;
108: {
109: register long x, y;
110: if(a<=0)return(0);
111: for(y=a,x=1;y!=0;y>>=2,x<<=1);
112: while((y=(a/x+x)>>1)<x)x=y;
113: return(x);
114: }
115: #define NSTACK 32
116: Hcoord hcoord(x, y, z, w)
117: fract x, y, z, w;
118: {
119: Hcoord r;
120: r.x=x;
121: r.y=y;
122: r.z=z;
123: r.w=w;
124: return(r);
125: }
126: matrix stack3[NSTACK]={
127: ONE,0,0,0,
128: 0,ONE,0,0,
129: 0,0,ONE,0,
130: 0,0,0,ONE,
131: };
132: matrix *sp3=stack3; /* matrix stack pointer */
133: /*
134: * Initialize the matrix stack.
135: * v is the viewport into which to map the image. (d/s) is the distance to the
136: * screen, whose half-width is 1. (equivalently, d is the distance to a screen of
137: * half-size s.)
138: */
139: Bitmap *viewport;
140: init3(v, d, s)
141: Bitmap *v;
142: {
143: sp3=stack3;
144: ident(*sp3);
145: (*sp3)[0][0]=(*sp3)[2][2]=d;
146: (*sp3)[1][1]=s;
147: (*sp3)[3][3]=s;
148: viewport=v;
149: cen=div(add(v->rect.origin, v->rect.corner), 2);
150: wid=min(v->rect.corner.x-v->rect.origin.x,
151: v->rect.corner.y-v->rect.origin.y)/2;
152: }
153: /*
154: * Push a copy of the top of the matrix stack onto the stack
155: */
156: push3(){
157: register i, j;
158: ;if(sp3==stack3+NSTACK-1)
159: return(-1);
160: for(i=0;i!=4;i++) for(j=0;j!=4;j++)
161: sp3[1][i][j]=sp3[0][i][j];
162: sp3++;
163: return(0);
164: }
165: /*
166: * pop the top matrix off the stack
167: */
168: pop3(){
169: if(sp3==stack3)
170: return(-1);
171: --sp3;
172: return(0);
173: }
174: /*
175: * Multiply the top matrix on the stack by rotation of angle theta about
176: * the given axis.
177: */
178: rot3(theta, axis)
179: angle theta;
180: {
181: rotsc3(isin(theta), icos(theta), axis);
182: }
183: /*
184: * rotate, given the sin and cosine of the rotation angle.
185: */
186: rotsc3(s, c, axis)
187: short s, c;
188: {
189: matrix m;
190: register i=(axis+1)%3, j=(axis+2)%3;
191: ident(m);
192: m[i][i] = c;
193: m[i][j] = -s;
194: m[j][i] = s;
195: m[j][j] = c;
196: xform3(m);
197: }
198: /*
199: * Scale by a factor of [xyz]/w about axis [xyz]
200: */
201: scale3(p)
202: Hcoord p;
203: {
204: matrix m;
205: ident(m);
206: m[0][0]=p.x;
207: m[1][1]=p.y;
208: m[2][2]=p.z;
209: m[3][3]=p.w;
210: xform3(m);
211: }
212: /*
213: * translate by (dx/dw, dy/dw, dz/dw)
214: */
215: tran3(p)
216: Hcoord p;
217: {
218: matrix m;
219: ident(m);
220: m[0][0]=p.w;
221: m[1][1]=p.w;
222: m[2][2]=p.w;
223: m[3][3]=p.w;
224: m[0][3]=p.x;
225: m[1][3]=p.y;
226: m[2][3]=p.z;
227: xform3(m);
228: }
229: /*
230: * Store an identity matrix in m
231: */
232: ident(m)
233: matrix m;
234: {
235: register fract *s = &m[0][0];
236: *s++=ONE;*s++=0;*s++=0;*s++=0;
237: *s++=0;*s++=ONE;*s++=0;*s++=0;
238: *s++=0;*s++=0;*s++=ONE;*s++=0;
239: *s++=0;*s++=0;*s++=0;*s++=ONE;
240: }
241: /*
242: * Multiply the top of the matrix stack by m
243: */
244: xform3(m)
245: matrix m;
246: {
247: register i, j, k;
248: long sum, max=0;
249: long tmp[4][4];
250: char buf[128];
251:
252: for(i=0;i!=4;i++) for(j=0;j!=4;j++){
253: sum=0;
254: for(k=0;k!=4;k++)
255: sum+=(long)(*sp3)[i][k]*m[k][j];
256: tmp[i][j]=sum;
257: if(sum<0)
258: sum = -sum;
259: if(sum>max)
260: max=sum;
261: }
262: if((k = max/(ONE-1)) == 0) k = 1;
263: for(i=0;i!=4;i++) for(j=0;j!=4;j++)
264: (*sp3)[i][j]=tmp[i][j]/k;
265: }
266: long dot(a, b)
267: Hcoord a, b;
268: {
269: return((long)a.x*b.x+(long)a.y*b.y+(long)a.z*b.z);
270: }
271: Hcoord unitize(a)
272: Hcoord a;
273: {
274: short l=isqrt(dot(a, a));
275: Hcoord v;
276: if(l==0){
277: v=a;
278: v.w=0;
279: }
280: else{
281: if(a.w<0) l = -l;
282: v.x=muldiv(a.x, ONE, l);
283: v.y=muldiv(a.y, ONE, l);
284: v.z=muldiv(a.z, ONE, l);
285: v.w=ONE;
286: }
287: return(v);
288: }
289: Hcoord cross(a, b)
290: Hcoord a, b;
291: {
292: Hcoord r;
293: r.x=((long)a.y*b.z-(long)a.z*b.y)/ONE;
294: r.y=((long)a.z*b.x-(long)a.x*b.z)/ONE;
295: r.z=((long)a.x*b.y-(long)a.y*b.x)/ONE;
296: r.w=ONE;
297: return(r);
298: }
299: /*
300: * multiply the top of the maxrix stack by a view-pointing transformation
301: * with the eyepoint at e, looking in direction (lx, ly, lz),
302: * with up direction (ux, uy, uz)
303: */
304: look3(e, l, u)
305: Hcoord e, l, u;
306: {
307: matrix m;
308: fract len;
309: Hcoord r;
310:
311: /*printf("look3(e=(%d,%d,%d,%d), l=(%d,%d,%d,%d), u=(%d,%d,%d,%d))\n",
312: e.x,e.y,e.z,e.w,l.x,l.y,l.z,l.w,u.x,u.y,u.z,u.w);*/
313: l=unitize(l);
314: if(l.w==0){
315: l.y=ONE;
316: l.w=ONE;
317: }
318: /* Adjust u to be perpendicular to l */
319: u=unitize(u);
320: len=dot(u, l)/ONE;
321: u.x-=muldiv(len, l.x, ONE);
322: u.y-=muldiv(len, l.y, ONE);
323: u.z-=muldiv(len, l.z, ONE);
324: u=unitize(u);
325: if(u.w==0){
326: u.z=ONE;
327: u.w=ONE;
328: }
329: r=cross(l, u);
330: /*printf("r=(%d,%d,%d,%d), l=(%d,%d,%d,%d), u=(%d,%d,%d,%d))\n",
331: r.x,r.y,r.z,r.w,l.x,l.y,l.z,l.w,u.x,u.y,u.z,u.w);*/
332: /* make the matrix to transform from (rlu) space to (xyz) space */
333: ident(m);
334: m[0][0]=r.x; m[0][1]=r.y; m[0][2]=r.z;
335: m[1][0]=l.x; m[1][1]=l.y; m[1][2]=l.z;
336: m[2][0]=u.x; m[2][1]=u.y; m[2][2]=u.z;
337: xform3(m);
338: tran3(hcoord(e.x, e.y, e.z, -e.w));
339: }
340: Hcoordl xformp(p)
341: Hcoord p;
342: {
343: Hcoordl r;
344: r.x=(long)(*sp3)[0][0]*p.x+(long)(*sp3)[0][1]*p.y
345: +(long)(*sp3)[0][2]*p.z+(long)(*sp3)[0][3]*p.w;
346: r.y=(long)(*sp3)[1][0]*p.x+(long)(*sp3)[1][1]*p.y
347: +(long)(*sp3)[1][2]*p.z+(long)(*sp3)[1][3]*p.w;
348: r.z=(long)(*sp3)[2][0]*p.x+(long)(*sp3)[2][1]*p.y
349: +(long)(*sp3)[2][2]*p.z+(long)(*sp3)[2][3]*p.w;
350: r.w=(long)(*sp3)[3][0]*p.x+(long)(*sp3)[3][1]*p.y
351: +(long)(*sp3)[3][2]*p.z+(long)(*sp3)[3][3]*p.w;
352: return(r);
353: }
354: move3(p)
355: Hcoord p;
356: {
357: cur=xformp(p);
358: /*printf("move3(%d,%d,%d,%d)\n", p.x, p.y, p.z, p.w);*/
359: }
360: line3(p)
361: Hcoord p;
362: {
363: Hcoordl new;
364: /*printf("line3(%d,%d,%d,%d)\n", p.x, p.y, p.z, p.w);*/
365: new=xformp(p);
366: clipanddraw(cur, new);
367: cur=new;
368: }
369: /*
370: * clip points to the viewing volume, then do perspective division
371: * and mapping to screen coordinates and draw lines. We only clip against the
372: * eye plane, trusting segment to clip in x and z. There probably should
373: * be some overflow precaution taken in the multiplication by n/d.
374: */
375: int clip3scale = 0;
376:
377: static clipanddraw(a, b)
378: Hcoordl a, b;
379: {
380: Point A, B;
381: Hcoordl t;
382: int clip=0;
383: long n, d;
384: long m;
385: long as, bs;
386: #define OV (32767) /* 2**(15.5) */
387: /*
388: * we clip against the plane p=(0,1,0,-1) [this corresponds
389: * to the euclidean plane y=1. This is picked rather than y=0
390: * to avoid a divide check in the perspective divide]
391: * First we check to see if one endpoint is on each side of the plane.
392: * and if necessary switch them so that a is behind and b is in front.
393: * then we replace a by a+(b-a)*n/d, where p.(a+(b-a)*n/d)=0
394: * i.e. a.y+(b.y-a.y)*n/d-a.w-(b.w-a.w)*n/d=0
395: * (b.y-a.y-b.w+a.w)*n/d=a.w-a.y
396: * n/d=(a.w-a.y)/(b.y-a.y-b.w+a.w)
397: */
398: /*printf("clipanddraw(a=(%d,%d,%d,%d), b=(%d,%d,%d,%d)\n",
399: a.x,a.y,a.z,a.w,b.x,b.y,b.z,b.w);*/
400: if(a.w>0?a.y<a.w:a.y>a.w) /* eqv. a.y/a.w<1, but avoids divide */
401: clip++;
402: if(b.w>0?b.y<b.w:b.y>b.w){
403: if(clip) /* both points behind y=1 plane */
404: {
405: /*printf("clipped out of existence\n");*/
406: return;
407: }
408: t=a;
409: a=b;
410: b=t;
411: clip++;
412: }
413: /* check for overflow */
414: m = max(max(abs(a.x), abs(a.y)), max(abs(a.z), abs(b.x)));
415: m = max(m, max(max(abs(b.y), abs(b.z)), max(abs(a.w), abs(b.w))));
416: if(m > OV)
417: {
418: m = (m+OV-1)/OV;
419: a.x /= m;
420: a.y /= m;
421: a.z /= m;
422: a.w /= m;
423: b.x /= m;
424: b.y /= m;
425: b.z /= m;
426: b.w /= m;
427: }
428: if(clip){
429: n=a.w-a.y;
430: d=b.y-a.y-b.w+a.w;
431: if(d==0) /* never happens? */
432: return;
433: a.x+=(b.x-a.x)*n/d;
434: a.y+=(b.y-a.y)*n/d;
435: a.z+=(b.z-a.z)*n/d;
436: a.w+=(b.w-a.w)*n/d; /* don't use w below */
437: /*printf("clipped to a=(%d,%d,%d,%d), b=(%d,%d,%d,%d)\n",
438: a.x,a.y,a.z,a.w,b.x,b.y,b.z,b.w);*/
439: }
440: if(clip3scale)
441: {
442: as = a.w*clip3scale;
443: bs = b.w*clip3scale;
444: }
445: else
446: {
447: as = a.y;
448: bs = b.y;
449: }
450: A=add(Pt((short)muldiv(a.x, wid, as), (short)muldiv(-a.z, wid, as)), cen);
451: B=add(Pt((short)muldiv(b.x, wid, bs), (short)muldiv(-b.z, wid, bs)), cen);
452: segment(viewport, A, B, F_OR);
453: /*printf("---->drew points (%d,%d) - (%d,%d)\n", A.x, A.y, B.x, B.y);*/
454: }
This archive runs on limited infrastructure. Preserving old code on modern bandwidth. Automated agents are requested to crawl responsibly.