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