1 /* $Id: project.c,v 1.2 1999/09/14 00:10:31 brianp Exp $ */
4 * Mesa 3-D graphics library
6 * Copyright (C) 1995-1999 Brian Paul
8 * This library is free software; you can redistribute it and/or
9 * modify it under the terms of the GNU Library General Public
10 * License as published by the Free Software Foundation; either
11 * version 2 of the License, or (at your option) any later version.
13 * This library is distributed in the hope that it will be useful,
14 * but WITHOUT ANY WARRANTY; without even the implied warranty of
15 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
16 * Library General Public License for more details.
18 * You should have received a copy of the GNU Library General Public
19 * License along with this library; if not, write to the Free
20 * Software Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
26 * Revision 1.2 1999/09/14 00:10:31 brianp
27 * added gluUnProject4()
29 * Revision 1.1.1.1 1999/08/19 00:55:42 jtg
32 * Revision 1.7 1999/01/03 03:23:15 brianp
33 * now using GLAPIENTRY and GLCALLBACK keywords (Ted Jump)
35 * Revision 1.6 1998/07/08 01:43:43 brianp
36 * new version of invert_matrix() (also in src/matrix.c)
38 * Revision 1.5 1997/07/24 01:28:44 brianp
39 * changed precompiled header symbol from PCH to PC_HEADER
41 * Revision 1.4 1997/05/28 02:29:38 brianp
42 * added support for precompiled headers (PCH), inserted APIENTRY keyword
44 * Revision 1.3 1997/04/11 23:22:42 brianp
45 * added divide by zero checks to gluProject() and gluUnproject()
47 * Revision 1.2 1997/01/29 19:05:29 brianp
48 * faster invert_matrix() function from Stephane Rehel
50 * Revision 1.1 1996/09/27 01:19:39 brianp
67 * This code was contributed by Marc Buffat (buffat@mecaflu.ec-lyon.fr).
73 /* implementation de gluProject et gluUnproject */
74 /* M. Buffat 17/2/95 */
79 * Transform a point (column vector) by a 4x4 matrix. I.e. out = m * in
80 * Input: m - the 4x4 matrix
82 * Output: out - the resulting 4x1 vector.
84 static void transform_point( GLdouble out
[4], const GLdouble m
[16],
85 const GLdouble in
[4] )
87 #define M(row,col) m[col*4+row]
88 out
[0] = M(0,0) * in
[0] + M(0,1) * in
[1] + M(0,2) * in
[2] + M(0,3) * in
[3];
89 out
[1] = M(1,0) * in
[0] + M(1,1) * in
[1] + M(1,2) * in
[2] + M(1,3) * in
[3];
90 out
[2] = M(2,0) * in
[0] + M(2,1) * in
[1] + M(2,2) * in
[2] + M(2,3) * in
[3];
91 out
[3] = M(3,0) * in
[0] + M(3,1) * in
[1] + M(3,2) * in
[2] + M(3,3) * in
[3];
99 * Perform a 4x4 matrix multiplication (product = a x b).
100 * Input: a, b - matrices to multiply
101 * Output: product - product of a and b
103 static void matmul( GLdouble
*product
, const GLdouble
*a
, const GLdouble
*b
)
105 /* This matmul was contributed by Thomas Malik */
109 #define A(row,col) a[(col<<2)+row]
110 #define B(row,col) b[(col<<2)+row]
111 #define T(row,col) temp[(col<<2)+row]
114 for (i
= 0; i
< 4; i
++)
116 T(i
, 0) = A(i
, 0) * B(0, 0) + A(i
, 1) * B(1, 0) + A(i
, 2) * B(2, 0) + A(i
, 3) * B(3, 0);
117 T(i
, 1) = A(i
, 0) * B(0, 1) + A(i
, 1) * B(1, 1) + A(i
, 2) * B(2, 1) + A(i
, 3) * B(3, 1);
118 T(i
, 2) = A(i
, 0) * B(0, 2) + A(i
, 1) * B(1, 2) + A(i
, 2) * B(2, 2) + A(i
, 3) * B(3, 2);
119 T(i
, 3) = A(i
, 0) * B(0, 3) + A(i
, 1) * B(1, 3) + A(i
, 2) * B(2, 3) + A(i
, 3) * B(3, 3);
125 MEMCPY( product
, temp
, 16*sizeof(GLdouble
) );
131 * Compute inverse of 4x4 transformation matrix.
132 * Code contributed by Jacques Leroy jle@star.be
133 * Return GL_TRUE for success, GL_FALSE for failure (singular matrix)
135 static GLboolean
invert_matrix( const GLdouble
*m
, GLdouble
*out
)
137 /* NB. OpenGL Matrices are COLUMN major. */
138 #define SWAP_ROWS(a, b) { GLdouble *_tmp = a; (a)=(b); (b)=_tmp; }
139 #define MAT(m,r,c) (m)[(c)*4+(r)]
142 GLdouble m0
, m1
, m2
, m3
, s
;
143 GLdouble
*r0
, *r1
, *r2
, *r3
;
145 r0
= wtmp
[0], r1
= wtmp
[1], r2
= wtmp
[2], r3
= wtmp
[3];
147 r0
[0] = MAT(m
,0,0), r0
[1] = MAT(m
,0,1),
148 r0
[2] = MAT(m
,0,2), r0
[3] = MAT(m
,0,3),
149 r0
[4] = 1.0, r0
[5] = r0
[6] = r0
[7] = 0.0,
151 r1
[0] = MAT(m
,1,0), r1
[1] = MAT(m
,1,1),
152 r1
[2] = MAT(m
,1,2), r1
[3] = MAT(m
,1,3),
153 r1
[5] = 1.0, r1
[4] = r1
[6] = r1
[7] = 0.0,
155 r2
[0] = MAT(m
,2,0), r2
[1] = MAT(m
,2,1),
156 r2
[2] = MAT(m
,2,2), r2
[3] = MAT(m
,2,3),
157 r2
[6] = 1.0, r2
[4] = r2
[5] = r2
[7] = 0.0,
159 r3
[0] = MAT(m
,3,0), r3
[1] = MAT(m
,3,1),
160 r3
[2] = MAT(m
,3,2), r3
[3] = MAT(m
,3,3),
161 r3
[7] = 1.0, r3
[4] = r3
[5] = r3
[6] = 0.0;
163 /* choose pivot - or die */
164 if (fabs(r3
[0])>fabs(r2
[0])) SWAP_ROWS(r3
, r2
);
165 if (fabs(r2
[0])>fabs(r1
[0])) SWAP_ROWS(r2
, r1
);
166 if (fabs(r1
[0])>fabs(r0
[0])) SWAP_ROWS(r1
, r0
);
167 if (0.0 == r0
[0]) return GL_FALSE
;
169 /* eliminate first variable */
170 m1
= r1
[0]/r0
[0]; m2
= r2
[0]/r0
[0]; m3
= r3
[0]/r0
[0];
171 s
= r0
[1]; r1
[1] -= m1
* s
; r2
[1] -= m2
* s
; r3
[1] -= m3
* s
;
172 s
= r0
[2]; r1
[2] -= m1
* s
; r2
[2] -= m2
* s
; r3
[2] -= m3
* s
;
173 s
= r0
[3]; r1
[3] -= m1
* s
; r2
[3] -= m2
* s
; r3
[3] -= m3
* s
;
175 if (s
!= 0.0) { r1
[4] -= m1
* s
; r2
[4] -= m2
* s
; r3
[4] -= m3
* s
; }
177 if (s
!= 0.0) { r1
[5] -= m1
* s
; r2
[5] -= m2
* s
; r3
[5] -= m3
* s
; }
179 if (s
!= 0.0) { r1
[6] -= m1
* s
; r2
[6] -= m2
* s
; r3
[6] -= m3
* s
; }
181 if (s
!= 0.0) { r1
[7] -= m1
* s
; r2
[7] -= m2
* s
; r3
[7] -= m3
* s
; }
183 /* choose pivot - or die */
184 if (fabs(r3
[1])>fabs(r2
[1])) SWAP_ROWS(r3
, r2
);
185 if (fabs(r2
[1])>fabs(r1
[1])) SWAP_ROWS(r2
, r1
);
186 if (0.0 == r1
[1]) return GL_FALSE
;
188 /* eliminate second variable */
189 m2
= r2
[1]/r1
[1]; m3
= r3
[1]/r1
[1];
190 r2
[2] -= m2
* r1
[2]; r3
[2] -= m3
* r1
[2];
191 r2
[3] -= m2
* r1
[3]; r3
[3] -= m3
* r1
[3];
192 s
= r1
[4]; if (0.0 != s
) { r2
[4] -= m2
* s
; r3
[4] -= m3
* s
; }
193 s
= r1
[5]; if (0.0 != s
) { r2
[5] -= m2
* s
; r3
[5] -= m3
* s
; }
194 s
= r1
[6]; if (0.0 != s
) { r2
[6] -= m2
* s
; r3
[6] -= m3
* s
; }
195 s
= r1
[7]; if (0.0 != s
) { r2
[7] -= m2
* s
; r3
[7] -= m3
* s
; }
197 /* choose pivot - or die */
198 if (fabs(r3
[2])>fabs(r2
[2])) SWAP_ROWS(r3
, r2
);
199 if (0.0 == r2
[2]) return GL_FALSE
;
201 /* eliminate third variable */
203 r3
[3] -= m3
* r2
[3], r3
[4] -= m3
* r2
[4],
204 r3
[5] -= m3
* r2
[5], r3
[6] -= m3
* r2
[6],
208 if (0.0 == r3
[3]) return GL_FALSE
;
210 s
= 1.0/r3
[3]; /* now back substitute row 3 */
211 r3
[4] *= s
; r3
[5] *= s
; r3
[6] *= s
; r3
[7] *= s
;
213 m2
= r2
[3]; /* now back substitute row 2 */
215 r2
[4] = s
* (r2
[4] - r3
[4] * m2
), r2
[5] = s
* (r2
[5] - r3
[5] * m2
),
216 r2
[6] = s
* (r2
[6] - r3
[6] * m2
), r2
[7] = s
* (r2
[7] - r3
[7] * m2
);
218 r1
[4] -= r3
[4] * m1
, r1
[5] -= r3
[5] * m1
,
219 r1
[6] -= r3
[6] * m1
, r1
[7] -= r3
[7] * m1
;
221 r0
[4] -= r3
[4] * m0
, r0
[5] -= r3
[5] * m0
,
222 r0
[6] -= r3
[6] * m0
, r0
[7] -= r3
[7] * m0
;
224 m1
= r1
[2]; /* now back substitute row 1 */
226 r1
[4] = s
* (r1
[4] - r2
[4] * m1
), r1
[5] = s
* (r1
[5] - r2
[5] * m1
),
227 r1
[6] = s
* (r1
[6] - r2
[6] * m1
), r1
[7] = s
* (r1
[7] - r2
[7] * m1
);
229 r0
[4] -= r2
[4] * m0
, r0
[5] -= r2
[5] * m0
,
230 r0
[6] -= r2
[6] * m0
, r0
[7] -= r2
[7] * m0
;
232 m0
= r0
[1]; /* now back substitute row 0 */
234 r0
[4] = s
* (r0
[4] - r1
[4] * m0
), r0
[5] = s
* (r0
[5] - r1
[5] * m0
),
235 r0
[6] = s
* (r0
[6] - r1
[6] * m0
), r0
[7] = s
* (r0
[7] - r1
[7] * m0
);
237 MAT(out
,0,0) = r0
[4]; MAT(out
,0,1) = r0
[5],
238 MAT(out
,0,2) = r0
[6]; MAT(out
,0,3) = r0
[7],
239 MAT(out
,1,0) = r1
[4]; MAT(out
,1,1) = r1
[5],
240 MAT(out
,1,2) = r1
[6]; MAT(out
,1,3) = r1
[7],
241 MAT(out
,2,0) = r2
[4]; MAT(out
,2,1) = r2
[5],
242 MAT(out
,2,2) = r2
[6]; MAT(out
,2,3) = r2
[7],
243 MAT(out
,3,0) = r3
[4]; MAT(out
,3,1) = r3
[5],
244 MAT(out
,3,2) = r3
[6]; MAT(out
,3,3) = r3
[7];
254 /* projection du point (objx,objy,obz) sur l'ecran (winx,winy,winz) */
255 GLint GLAPIENTRY
gluProject(GLdouble objx
,GLdouble objy
,GLdouble objz
,
256 const GLdouble model
[16],const GLdouble proj
[16],
257 const GLint viewport
[4],
258 GLdouble
*winx
,GLdouble
*winy
,GLdouble
*winz
)
260 /* matrice de transformation */
261 GLdouble in
[4],out
[4];
263 /* initilise la matrice et le vecteur a transformer */
264 in
[0]=objx
; in
[1]=objy
; in
[2]=objz
; in
[3]=1.0;
265 transform_point(out
,model
,in
);
266 transform_point(in
,proj
,out
);
268 /* d'ou le resultat normalise entre -1 et 1*/
272 in
[0]/=in
[3]; in
[1]/=in
[3]; in
[2]/=in
[3];
274 /* en coordonnees ecran */
275 *winx
= viewport
[0]+(1+in
[0])*viewport
[2]/2;
276 *winy
= viewport
[1]+(1+in
[1])*viewport
[3]/2;
277 /* entre 0 et 1 suivant z */
284 /* transformation du point ecran (winx,winy,winz) en point objet */
285 GLint GLAPIENTRY
gluUnProject(GLdouble winx
,GLdouble winy
,GLdouble winz
,
286 const GLdouble model
[16],const GLdouble proj
[16],
287 const GLint viewport
[4],
288 GLdouble
*objx
,GLdouble
*objy
,GLdouble
*objz
)
290 /* matrice de transformation */
291 GLdouble m
[16], A
[16];
292 GLdouble in
[4],out
[4];
294 /* transformation coordonnees normalisees entre -1 et 1 */
295 in
[0]=(winx
-viewport
[0])*2/viewport
[2] - 1.0;
296 in
[1]=(winy
-viewport
[1])*2/viewport
[3] - 1.0;
300 /* calcul transformation inverse */
301 matmul(A
,proj
,model
);
304 /* d'ou les coordonnees objets */
305 transform_point(out
,m
,in
);
317 * This is like gluUnProject but also takes near and far DepthRange values.
320 gluUnProject4( GLdouble winx
, GLdouble winy
, GLdouble winz
, GLdouble clipw
,
321 const GLdouble modelMatrix
[16],
322 const GLdouble projMatrix
[16],
323 const GLint viewport
[4],
324 GLclampd nearZ
, GLclampd farZ
,
325 GLdouble
*objx
, GLdouble
*objy
, GLdouble
*objz
, GLdouble
*objw
)
327 /* matrice de transformation */
328 GLdouble m
[16], A
[16];
329 GLdouble in
[4],out
[4];
330 GLdouble z
= nearZ
+ winz
* (farZ
- nearZ
);
332 /* transformation coordonnees normalisees entre -1 et 1 */
333 in
[0] = (winx
-viewport
[0])*2/viewport
[2] - 1.0;
334 in
[1] = (winy
-viewport
[1])*2/viewport
[3] - 1.0;
335 in
[2] = 2.0 * z
- 1.0;
338 /* calcul transformation inverse */
339 matmul(A
,projMatrix
,modelMatrix
);
342 /* d'ou les coordonnees objets */
343 transform_point(out
,m
,in
);