sys/src/libgeometry/tstack.c

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