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