sys/src/libgeometry/transform.c

   1 /*
   2  * The following routines transform points and planes from one space
   3  * to another.  Points and planes are represented by their
   4  * homogeneous coordinates, stored in variables of type Point3.
   5  */
   6 #include <u.h>
   7 #include <libc.h>
   8 #include <draw.h>
   9 #include <geometry.h>
  10 /*
  11  * Transform point p.
  12  */
  13 Point3 xformpoint(Point3 p, Space *to, Space *from){
  14         Point3 q, r;
  15         register double *m;
  16         if(from){
  17                 m=&from->t[0][0];
  18                 q.x=*m++*p.x; q.x+=*m++*p.y; q.x+=*m++*p.z; q.x+=*m++*p.w;
  19                 q.y=*m++*p.x; q.y+=*m++*p.y; q.y+=*m++*p.z; q.y+=*m++*p.w;
  20                 q.z=*m++*p.x; q.z+=*m++*p.y; q.z+=*m++*p.z; q.z+=*m++*p.w;
  21                 q.w=*m++*p.x; q.w+=*m++*p.y; q.w+=*m++*p.z; q.w+=*m  *p.w;
  22         }
  23         else
  24                 q=p;
  25         if(to){
  26                 m=&to->tinv[0][0];
  27                 r.x=*m++*q.x; r.x+=*m++*q.y; r.x+=*m++*q.z; r.x+=*m++*q.w;
  28                 r.y=*m++*q.x; r.y+=*m++*q.y; r.y+=*m++*q.z; r.y+=*m++*q.w;
  29                 r.z=*m++*q.x; r.z+=*m++*q.y; r.z+=*m++*q.z; r.z+=*m++*q.w;
  30                 r.w=*m++*q.x; r.w+=*m++*q.y; r.w+=*m++*q.z; r.w+=*m  *q.w;
  31         }
  32         else
  33                 r=q;
  34         return r;
  35 }
  36 /*
  37  * Transform point p with perspective division.
  38  */
  39 Point3 xformpointd(Point3 p, Space *to, Space *from){
  40         p=xformpoint(p, to, from);
  41         if(p.w!=0){
  42                 p.x/=p.w;
  43                 p.y/=p.w;
  44                 p.z/=p.w;
  45                 p.w=1;
  46         }
  47         return p;
  48 }
  49 /*
  50  * Transform plane p -- same as xformpoint, except multiply on the
  51  * other side by the inverse matrix.
  52  */
  53 Point3 xformplane(Point3 p, Space *to, Space *from){
  54         Point3 q, r;
  55         register double *m;
  56         if(from){
  57                 m=&from->tinv[0][0];
  58                 q.x =*m++*p.x; q.y =*m++*p.x; q.z =*m++*p.x; q.w =*m++*p.x;
  59                 q.x+=*m++*p.y; q.y+=*m++*p.y; q.z+=*m++*p.y; q.w+=*m++*p.y;
  60                 q.x+=*m++*p.z; q.y+=*m++*p.z; q.z+=*m++*p.z; q.w+=*m++*p.z;
  61                 q.x+=*m++*p.w; q.y+=*m++*p.w; q.z+=*m++*p.w; q.w+=*m  *p.w;
  62         }
  63         else
  64                 q=p;
  65         if(to){
  66                 m=&to->t[0][0];
  67                 r.x =*m++*q.x; r.y =*m++*q.x; r.z =*m++*q.x; r.w =*m++*q.x;
  68                 r.x+=*m++*q.y; r.y+=*m++*q.y; r.z+=*m++*q.y; r.w+=*m++*q.y;
  69                 r.x+=*m++*q.z; r.y+=*m++*q.z; r.z+=*m++*q.z; r.w+=*m++*q.z;
  70                 r.x+=*m++*q.w; r.y+=*m++*q.w; r.z+=*m++*q.w; r.w+=*m  *q.w;
  71         }
  72         else
  73                 r=q;
  74         return r;
  75 }