1 d1e9002f 2005-01-04 devnull /*% cc -gpc %
2 d1e9002f 2005-01-04 devnull * These transformation routines maintain stacks of transformations
3 fa325e9b 2020-01-10 cross * and their inverses.
4 d1e9002f 2005-01-04 devnull * t=pushmat(t) push matrix stack
5 d1e9002f 2005-01-04 devnull * t=popmat(t) pop matrix stack
6 d1e9002f 2005-01-04 devnull * rot(t, a, axis) multiply stack top by rotation
7 d1e9002f 2005-01-04 devnull * qrot(t, q) multiply stack top by rotation, q is unit quaternion
8 d1e9002f 2005-01-04 devnull * scale(t, x, y, z) multiply stack top by scale
9 d1e9002f 2005-01-04 devnull * move(t, x, y, z) multiply stack top by translation
10 d1e9002f 2005-01-04 devnull * xform(t, m) multiply stack top by m
11 d1e9002f 2005-01-04 devnull * ixform(t, m, inv) multiply stack top by m. inv is the inverse of m.
12 d1e9002f 2005-01-04 devnull * look(t, e, l, u) multiply stack top by viewing transformation
13 d1e9002f 2005-01-04 devnull * persp(t, fov, n, f) multiply stack top by perspective transformation
14 d1e9002f 2005-01-04 devnull * viewport(t, r, aspect)
15 d1e9002f 2005-01-04 devnull * multiply stack top by window->viewport transformation.
17 d1e9002f 2005-01-04 devnull #include <u.h>
18 d1e9002f 2005-01-04 devnull #include <libc.h>
19 d1e9002f 2005-01-04 devnull #include <draw.h>
20 d1e9002f 2005-01-04 devnull #include <geometry.h>
21 d1e9002f 2005-01-04 devnull Space *pushmat(Space *t){
22 d1e9002f 2005-01-04 devnull Space *v;
23 d1e9002f 2005-01-04 devnull v=malloc(sizeof(Space));
24 d1e9002f 2005-01-04 devnull if(t==0){
25 d1e9002f 2005-01-04 devnull ident(v->t);
26 d1e9002f 2005-01-04 devnull ident(v->tinv);
30 d1e9002f 2005-01-04 devnull v->next=t;
31 d1e9002f 2005-01-04 devnull return v;
33 d1e9002f 2005-01-04 devnull Space *popmat(Space *t){
34 d1e9002f 2005-01-04 devnull Space *v;
35 d1e9002f 2005-01-04 devnull if(t==0) return 0;
36 d1e9002f 2005-01-04 devnull v=t->next;
38 d1e9002f 2005-01-04 devnull return v;
40 d1e9002f 2005-01-04 devnull void rot(Space *t, double theta, int axis){
41 d1e9002f 2005-01-04 devnull double s=sin(radians(theta)), c=cos(radians(theta));
42 d1e9002f 2005-01-04 devnull Matrix m, inv;
43 d1e9002f 2005-01-04 devnull int i=(axis+1)%3, j=(axis+2)%3;
44 d1e9002f 2005-01-04 devnull ident(m);
45 d1e9002f 2005-01-04 devnull m[i][i] = c;
46 d1e9002f 2005-01-04 devnull m[i][j] = -s;
47 d1e9002f 2005-01-04 devnull m[j][i] = s;
48 d1e9002f 2005-01-04 devnull m[j][j] = c;
49 d1e9002f 2005-01-04 devnull ident(inv);
50 d1e9002f 2005-01-04 devnull inv[i][i] = c;
51 d1e9002f 2005-01-04 devnull inv[i][j] = s;
52 d1e9002f 2005-01-04 devnull inv[j][i] = -s;
53 d1e9002f 2005-01-04 devnull inv[j][j] = c;
54 d1e9002f 2005-01-04 devnull ixform(t, m, inv);
56 d1e9002f 2005-01-04 devnull void qrot(Space *t, Quaternion q){
57 d1e9002f 2005-01-04 devnull Matrix m, inv;
58 d1e9002f 2005-01-04 devnull int i, j;
59 d1e9002f 2005-01-04 devnull qtom(m, q);
60 d1e9002f 2005-01-04 devnull for(i=0;i!=4;i++) for(j=0;j!=4;j++) inv[i][j]=m[j][i];
61 d1e9002f 2005-01-04 devnull ixform(t, m, inv);
63 d1e9002f 2005-01-04 devnull void scale(Space *t, double x, double y, double z){
64 d1e9002f 2005-01-04 devnull Matrix m, inv;
65 d1e9002f 2005-01-04 devnull ident(m);
66 d1e9002f 2005-01-04 devnull m[0][0]=x;
67 d1e9002f 2005-01-04 devnull m[1][1]=y;
68 d1e9002f 2005-01-04 devnull m[2][2]=z;
69 d1e9002f 2005-01-04 devnull ident(inv);
70 d1e9002f 2005-01-04 devnull inv[0][0]=1/x;
71 d1e9002f 2005-01-04 devnull inv[1][1]=1/y;
72 d1e9002f 2005-01-04 devnull inv[2][2]=1/z;
73 d1e9002f 2005-01-04 devnull ixform(t, m, inv);
75 d1e9002f 2005-01-04 devnull void move(Space *t, double x, double y, double z){
76 d1e9002f 2005-01-04 devnull Matrix m, inv;
77 d1e9002f 2005-01-04 devnull ident(m);
78 d1e9002f 2005-01-04 devnull m[0][3]=x;
79 d1e9002f 2005-01-04 devnull m[1][3]=y;
80 d1e9002f 2005-01-04 devnull m[2][3]=z;
81 d1e9002f 2005-01-04 devnull ident(inv);
82 d1e9002f 2005-01-04 devnull inv[0][3]=-x;
83 d1e9002f 2005-01-04 devnull inv[1][3]=-y;
84 d1e9002f 2005-01-04 devnull inv[2][3]=-z;
85 d1e9002f 2005-01-04 devnull ixform(t, m, inv);
87 d1e9002f 2005-01-04 devnull void xform(Space *t, Matrix m){
88 d1e9002f 2005-01-04 devnull Matrix inv;
89 d1e9002f 2005-01-04 devnull if(invertmat(m, inv)==0) return;
90 d1e9002f 2005-01-04 devnull ixform(t, m, inv);
92 d1e9002f 2005-01-04 devnull void ixform(Space *t, Matrix m, Matrix inv){
93 d1e9002f 2005-01-04 devnull matmul(t->t, m);
94 d1e9002f 2005-01-04 devnull matmulr(t->tinv, inv);
97 d1e9002f 2005-01-04 devnull * multiply the top of the matrix stack by a view-pointing transformation
98 d1e9002f 2005-01-04 devnull * with the eyepoint at e, looking at point l, with u at the top of the screen.
99 d1e9002f 2005-01-04 devnull * The coordinate system is deemed to be right-handed.
100 d1e9002f 2005-01-04 devnull * The generated transformation transforms this view into a view from
101 d1e9002f 2005-01-04 devnull * the origin, looking in the positive y direction, with the z axis pointing up,
102 d1e9002f 2005-01-04 devnull * and x to the right.
104 d1e9002f 2005-01-04 devnull void look(Space *t, Point3 e, Point3 l, Point3 u){
105 d1e9002f 2005-01-04 devnull Matrix m, inv;
106 d1e9002f 2005-01-04 devnull Point3 r;
107 d1e9002f 2005-01-04 devnull l=unit3(sub3(l, e));
108 d1e9002f 2005-01-04 devnull u=unit3(vrem3(sub3(u, e), l));
109 d1e9002f 2005-01-04 devnull r=cross3(l, u);
110 d1e9002f 2005-01-04 devnull /* make the matrix to transform from (rlu) space to (xyz) space */
111 d1e9002f 2005-01-04 devnull ident(m);
112 d1e9002f 2005-01-04 devnull m[0][0]=r.x; m[0][1]=r.y; m[0][2]=r.z;
113 d1e9002f 2005-01-04 devnull m[1][0]=l.x; m[1][1]=l.y; m[1][2]=l.z;
114 d1e9002f 2005-01-04 devnull m[2][0]=u.x; m[2][1]=u.y; m[2][2]=u.z;
115 d1e9002f 2005-01-04 devnull ident(inv);
116 d1e9002f 2005-01-04 devnull inv[0][0]=r.x; inv[0][1]=l.x; inv[0][2]=u.x;
117 d1e9002f 2005-01-04 devnull inv[1][0]=r.y; inv[1][1]=l.y; inv[1][2]=u.y;
118 d1e9002f 2005-01-04 devnull inv[2][0]=r.z; inv[2][1]=l.z; inv[2][2]=u.z;
119 d1e9002f 2005-01-04 devnull ixform(t, m, inv);
120 d1e9002f 2005-01-04 devnull move(t, -e.x, -e.y, -e.z);
123 d1e9002f 2005-01-04 devnull * generate a transformation that maps the frustum with apex at the origin,
124 d1e9002f 2005-01-04 devnull * apex angle=fov and clipping planes y=n and y=f into the double-unit cube.
125 d1e9002f 2005-01-04 devnull * plane y=n maps to y'=-1, y=f maps to y'=1
127 d1e9002f 2005-01-04 devnull int persp(Space *t, double fov, double n, double f){
128 d1e9002f 2005-01-04 devnull Matrix m;
129 d1e9002f 2005-01-04 devnull double z;
130 d1e9002f 2005-01-04 devnull if(n<=0 || f<=n || fov<=0 || 180<=fov) /* really need f!=n && sin(v)!=0 */
131 d1e9002f 2005-01-04 devnull return -1;
132 d1e9002f 2005-01-04 devnull z=1/tan(radians(fov)/2);
133 d1e9002f 2005-01-04 devnull m[0][0]=z; m[0][1]=0; m[0][2]=0; m[0][3]=0;
134 d1e9002f 2005-01-04 devnull m[1][0]=0; m[1][1]=(f+n)/(f-n); m[1][2]=0; m[1][3]=f*(1-m[1][1]);
135 d1e9002f 2005-01-04 devnull m[2][0]=0; m[2][1]=0; m[2][2]=z; m[2][3]=0;
136 d1e9002f 2005-01-04 devnull m[3][0]=0; m[3][1]=1; m[3][2]=0; m[3][3]=0;
137 d1e9002f 2005-01-04 devnull xform(t, m);
138 d1e9002f 2005-01-04 devnull return 0;
141 d1e9002f 2005-01-04 devnull * Map the unit-cube window into the given screen viewport.
142 d1e9002f 2005-01-04 devnull * r has min at the top left, max just outside the lower right. Aspect is the
143 d1e9002f 2005-01-04 devnull * aspect ratio (dx/dy) of the viewport's pixels (not of the whole viewport!)
144 d1e9002f 2005-01-04 devnull * The whole window is transformed to fit centered inside the viewport with equal
145 d1e9002f 2005-01-04 devnull * slop on either top and bottom or left and right, depending on the viewport's
146 d1e9002f 2005-01-04 devnull * aspect ratio.
147 d1e9002f 2005-01-04 devnull * The window is viewed down the y axis, with x to the left and z up. The viewport
148 d1e9002f 2005-01-04 devnull * has x increasing to the right and y increasing down. The window's y coordinates
149 d1e9002f 2005-01-04 devnull * are mapped, unchanged, into the viewport's z coordinates.
151 d1e9002f 2005-01-04 devnull void viewport(Space *t, Rectangle r, double aspect){
152 d1e9002f 2005-01-04 devnull Matrix m;
153 d1e9002f 2005-01-04 devnull double xc, yc, wid, hgt, scale;
154 d1e9002f 2005-01-04 devnull xc=.5*(r.min.x+r.max.x);
155 d1e9002f 2005-01-04 devnull yc=.5*(r.min.y+r.max.y);
156 d1e9002f 2005-01-04 devnull wid=(r.max.x-r.min.x)*aspect;
157 d1e9002f 2005-01-04 devnull hgt=r.max.y-r.min.y;
158 d1e9002f 2005-01-04 devnull scale=.5*(wid<hgt?wid:hgt);
159 d1e9002f 2005-01-04 devnull ident(m);
160 d1e9002f 2005-01-04 devnull m[0][0]=scale;
161 d1e9002f 2005-01-04 devnull m[0][3]=xc;
162 d1e9002f 2005-01-04 devnull m[1][1]=0;
163 d1e9002f 2005-01-04 devnull m[1][2]=-scale;
164 d1e9002f 2005-01-04 devnull m[1][3]=yc;
165 d1e9002f 2005-01-04 devnull m[2][1]=1;
166 d1e9002f 2005-01-04 devnull m[2][2]=0;
167 d1e9002f 2005-01-04 devnull /* should get inverse by hand */
168 d1e9002f 2005-01-04 devnull xform(t, m);
