new, faster version of __gluInvertMatrixd(), bug 6748
[mesa.git] / src / glu / sgi / libutil / project.c
1 /*
2 ** License Applicability. Except to the extent portions of this file are
3 ** made subject to an alternative license as permitted in the SGI Free
4 ** Software License B, Version 1.1 (the "License"), the contents of this
5 ** file are subject only to the provisions of the License. You may not use
6 ** this file except in compliance with the License. You may obtain a copy
7 ** of the License at Silicon Graphics, Inc., attn: Legal Services, 1600
8 ** Amphitheatre Parkway, Mountain View, CA 94043-1351, or at:
9 **
10 ** http://oss.sgi.com/projects/FreeB
11 **
12 ** Note that, as provided in the License, the Software is distributed on an
13 ** "AS IS" basis, with ALL EXPRESS AND IMPLIED WARRANTIES AND CONDITIONS
14 ** DISCLAIMED, INCLUDING, WITHOUT LIMITATION, ANY IMPLIED WARRANTIES AND
15 ** CONDITIONS OF MERCHANTABILITY, SATISFACTORY QUALITY, FITNESS FOR A
16 ** PARTICULAR PURPOSE, AND NON-INFRINGEMENT.
17 **
18 ** Original Code. The Original Code is: OpenGL Sample Implementation,
19 ** Version 1.2.1, released January 26, 2000, developed by Silicon Graphics,
20 ** Inc. The Original Code is Copyright (c) 1991-2000 Silicon Graphics, Inc.
21 ** Copyright in any portions created by third parties is as indicated
22 ** elsewhere herein. All Rights Reserved.
23 **
24 ** Additional Notice Provisions: The application programming interfaces
25 ** established by SGI in conjunction with the Original Code are The
26 ** OpenGL(R) Graphics System: A Specification (Version 1.2.1), released
27 ** April 1, 1999; The OpenGL(R) Graphics System Utility Library (Version
28 ** 1.3), released November 4, 1998; and OpenGL(R) Graphics with the X
29 ** Window System(R) (Version 1.3), released October 19, 1998. This software
30 ** was created using the OpenGL(R) version 1.2.1 Sample Implementation
31 ** published by SGI, but has not been independently verified as being
32 ** compliant with the OpenGL(R) version 1.2.1 Specification.
33 **
34 ** $Date: 2006/05/01 16:01:17 $ $Revision: 1.6 $
35 ** $Header: /home/krh/git/sync/mesa-cvs-repo/Mesa/src/glu/sgi/libutil/project.c,v 1.6 2006/05/01 16:01:17 brianp Exp $
36 */
37
38 #include "gluos.h"
39 #include <math.h>
40 #include <GL/gl.h>
41 #include <GL/glu.h>
42 #include "gluint.h"
43
44 /*
45 ** Make m an identity matrix
46 */
47 static void __gluMakeIdentityd(GLdouble m[16])
48 {
49 m[0+4*0] = 1; m[0+4*1] = 0; m[0+4*2] = 0; m[0+4*3] = 0;
50 m[1+4*0] = 0; m[1+4*1] = 1; m[1+4*2] = 0; m[1+4*3] = 0;
51 m[2+4*0] = 0; m[2+4*1] = 0; m[2+4*2] = 1; m[2+4*3] = 0;
52 m[3+4*0] = 0; m[3+4*1] = 0; m[3+4*2] = 0; m[3+4*3] = 1;
53 }
54
55 static void __gluMakeIdentityf(GLfloat m[16])
56 {
57 m[0+4*0] = 1; m[0+4*1] = 0; m[0+4*2] = 0; m[0+4*3] = 0;
58 m[1+4*0] = 0; m[1+4*1] = 1; m[1+4*2] = 0; m[1+4*3] = 0;
59 m[2+4*0] = 0; m[2+4*1] = 0; m[2+4*2] = 1; m[2+4*3] = 0;
60 m[3+4*0] = 0; m[3+4*1] = 0; m[3+4*2] = 0; m[3+4*3] = 1;
61 }
62
63 void GLAPIENTRY
64 gluOrtho2D(GLdouble left, GLdouble right, GLdouble bottom, GLdouble top)
65 {
66 glOrtho(left, right, bottom, top, -1, 1);
67 }
68
69 #define __glPi 3.14159265358979323846
70
71 void GLAPIENTRY
72 gluPerspective(GLdouble fovy, GLdouble aspect, GLdouble zNear, GLdouble zFar)
73 {
74 GLdouble m[4][4];
75 double sine, cotangent, deltaZ;
76 double radians = fovy / 2 * __glPi / 180;
77
78 deltaZ = zFar - zNear;
79 sine = sin(radians);
80 if ((deltaZ == 0) || (sine == 0) || (aspect == 0)) {
81 return;
82 }
83 cotangent = COS(radians) / sine;
84
85 __gluMakeIdentityd(&m[0][0]);
86 m[0][0] = cotangent / aspect;
87 m[1][1] = cotangent;
88 m[2][2] = -(zFar + zNear) / deltaZ;
89 m[2][3] = -1;
90 m[3][2] = -2 * zNear * zFar / deltaZ;
91 m[3][3] = 0;
92 glMultMatrixd(&m[0][0]);
93 }
94
95 static void normalize(float v[3])
96 {
97 float r;
98
99 r = sqrt( v[0]*v[0] + v[1]*v[1] + v[2]*v[2] );
100 if (r == 0.0) return;
101
102 v[0] /= r;
103 v[1] /= r;
104 v[2] /= r;
105 }
106
107 static void cross(float v1[3], float v2[3], float result[3])
108 {
109 result[0] = v1[1]*v2[2] - v1[2]*v2[1];
110 result[1] = v1[2]*v2[0] - v1[0]*v2[2];
111 result[2] = v1[0]*v2[1] - v1[1]*v2[0];
112 }
113
114 void GLAPIENTRY
115 gluLookAt(GLdouble eyex, GLdouble eyey, GLdouble eyez, GLdouble centerx,
116 GLdouble centery, GLdouble centerz, GLdouble upx, GLdouble upy,
117 GLdouble upz)
118 {
119 float forward[3], side[3], up[3];
120 GLfloat m[4][4];
121
122 forward[0] = centerx - eyex;
123 forward[1] = centery - eyey;
124 forward[2] = centerz - eyez;
125
126 up[0] = upx;
127 up[1] = upy;
128 up[2] = upz;
129
130 normalize(forward);
131
132 /* Side = forward x up */
133 cross(forward, up, side);
134 normalize(side);
135
136 /* Recompute up as: up = side x forward */
137 cross(side, forward, up);
138
139 __gluMakeIdentityf(&m[0][0]);
140 m[0][0] = side[0];
141 m[1][0] = side[1];
142 m[2][0] = side[2];
143
144 m[0][1] = up[0];
145 m[1][1] = up[1];
146 m[2][1] = up[2];
147
148 m[0][2] = -forward[0];
149 m[1][2] = -forward[1];
150 m[2][2] = -forward[2];
151
152 glMultMatrixf(&m[0][0]);
153 glTranslated(-eyex, -eyey, -eyez);
154 }
155
156 static void __gluMultMatrixVecd(const GLdouble matrix[16], const GLdouble in[4],
157 GLdouble out[4])
158 {
159 int i;
160
161 for (i=0; i<4; i++) {
162 out[i] =
163 in[0] * matrix[0*4+i] +
164 in[1] * matrix[1*4+i] +
165 in[2] * matrix[2*4+i] +
166 in[3] * matrix[3*4+i];
167 }
168 }
169
170 /*
171 ** inverse = invert(src)
172 ** New, faster implementation by Shan Hao Bo, April 2006.
173 */
174 static int __gluInvertMatrixd(const GLdouble src[16], GLdouble inverse[16])
175 {
176 int i, j, k;
177 double t;
178 GLdouble temp[4][4];
179
180 for (i=0; i<4; i++) {
181 for (j=0; j<4; j++) {
182 temp[i][j] = src[i*4+j];
183 }
184 }
185 __gluMakeIdentityd(inverse);
186
187 for (i = 0; i < 4; i++) {
188 if (temp[i][i] == 0.0f) {
189 /*
190 ** Look for non-zero element in column
191 */
192 for (j = i + 1; j < 4; j++) {
193 if (temp[j][i] != 0.0f) {
194 break;
195 }
196 }
197
198 if (j != 4) {
199 /*
200 ** Swap rows.
201 */
202 for (k = 0; k < 4; k++) {
203 t = temp[i][k];
204 temp[i][k] = temp[j][k];
205 temp[j][k] = t;
206
207 t = inverse[i*4+k];
208 inverse[i*4+k] = inverse[j*4+k];
209 inverse[j*4+k] = t;
210 }
211 }
212 else {
213 /*
214 ** No non-zero pivot. The matrix is singular,
215 which shouldn't
216 ** happen. This means the user gave us a bad
217 matrix.
218 */
219 return GL_FALSE;
220 }
221 }
222
223 t = 1.0f / temp[i][i];
224 for (k = 0; k < 4; k++) {
225 temp[i][k] *= t;
226 inverse[i*4+k] *= t;
227 }
228 for (j = 0; j < 4; j++) {
229 if (j != i) {
230 t = temp[j][i];
231 for (k = 0; k < 4; k++) {
232 temp[j][k] -= temp[i][k]*t;
233 inverse[j*4+k] -= inverse[i*4+k]*t;
234 }
235 }
236 }
237 }
238 return GL_TRUE;
239 }
240
241 static void __gluMultMatricesd(const GLdouble a[16], const GLdouble b[16],
242 GLdouble r[16])
243 {
244 int i, j;
245
246 for (i = 0; i < 4; i++) {
247 for (j = 0; j < 4; j++) {
248 r[i*4+j] =
249 a[i*4+0]*b[0*4+j] +
250 a[i*4+1]*b[1*4+j] +
251 a[i*4+2]*b[2*4+j] +
252 a[i*4+3]*b[3*4+j];
253 }
254 }
255 }
256
257 GLint GLAPIENTRY
258 gluProject(GLdouble objx, GLdouble objy, GLdouble objz,
259 const GLdouble modelMatrix[16],
260 const GLdouble projMatrix[16],
261 const GLint viewport[4],
262 GLdouble *winx, GLdouble *winy, GLdouble *winz)
263 {
264 double in[4];
265 double out[4];
266
267 in[0]=objx;
268 in[1]=objy;
269 in[2]=objz;
270 in[3]=1.0;
271 __gluMultMatrixVecd(modelMatrix, in, out);
272 __gluMultMatrixVecd(projMatrix, out, in);
273 if (in[3] == 0.0) return(GL_FALSE);
274 in[0] /= in[3];
275 in[1] /= in[3];
276 in[2] /= in[3];
277 /* Map x, y and z to range 0-1 */
278 in[0] = in[0] * 0.5 + 0.5;
279 in[1] = in[1] * 0.5 + 0.5;
280 in[2] = in[2] * 0.5 + 0.5;
281
282 /* Map x,y to viewport */
283 in[0] = in[0] * viewport[2] + viewport[0];
284 in[1] = in[1] * viewport[3] + viewport[1];
285
286 *winx=in[0];
287 *winy=in[1];
288 *winz=in[2];
289 return(GL_TRUE);
290 }
291
292 GLint GLAPIENTRY
293 gluUnProject(GLdouble winx, GLdouble winy, GLdouble winz,
294 const GLdouble modelMatrix[16],
295 const GLdouble projMatrix[16],
296 const GLint viewport[4],
297 GLdouble *objx, GLdouble *objy, GLdouble *objz)
298 {
299 double finalMatrix[16];
300 double in[4];
301 double out[4];
302
303 __gluMultMatricesd(modelMatrix, projMatrix, finalMatrix);
304 if (!__gluInvertMatrixd(finalMatrix, finalMatrix)) return(GL_FALSE);
305
306 in[0]=winx;
307 in[1]=winy;
308 in[2]=winz;
309 in[3]=1.0;
310
311 /* Map x and y from window coordinates */
312 in[0] = (in[0] - viewport[0]) / viewport[2];
313 in[1] = (in[1] - viewport[1]) / viewport[3];
314
315 /* Map to range -1 to 1 */
316 in[0] = in[0] * 2 - 1;
317 in[1] = in[1] * 2 - 1;
318 in[2] = in[2] * 2 - 1;
319
320 __gluMultMatrixVecd(finalMatrix, in, out);
321 if (out[3] == 0.0) return(GL_FALSE);
322 out[0] /= out[3];
323 out[1] /= out[3];
324 out[2] /= out[3];
325 *objx = out[0];
326 *objy = out[1];
327 *objz = out[2];
328 return(GL_TRUE);
329 }
330
331 GLint GLAPIENTRY
332 gluUnProject4(GLdouble winx, GLdouble winy, GLdouble winz, GLdouble clipw,
333 const GLdouble modelMatrix[16],
334 const GLdouble projMatrix[16],
335 const GLint viewport[4],
336 GLclampd nearVal, GLclampd farVal,
337 GLdouble *objx, GLdouble *objy, GLdouble *objz,
338 GLdouble *objw)
339 {
340 double finalMatrix[16];
341 double in[4];
342 double out[4];
343
344 __gluMultMatricesd(modelMatrix, projMatrix, finalMatrix);
345 if (!__gluInvertMatrixd(finalMatrix, finalMatrix)) return(GL_FALSE);
346
347 in[0]=winx;
348 in[1]=winy;
349 in[2]=winz;
350 in[3]=clipw;
351
352 /* Map x and y from window coordinates */
353 in[0] = (in[0] - viewport[0]) / viewport[2];
354 in[1] = (in[1] - viewport[1]) / viewport[3];
355 in[2] = (in[2] - nearVal) / (farVal - nearVal);
356
357 /* Map to range -1 to 1 */
358 in[0] = in[0] * 2 - 1;
359 in[1] = in[1] * 2 - 1;
360 in[2] = in[2] * 2 - 1;
361
362 __gluMultMatrixVecd(finalMatrix, in, out);
363 if (out[3] == 0.0) return(GL_FALSE);
364 *objx = out[0];
365 *objy = out[1];
366 *objz = out[2];
367 *objw = out[3];
368 return(GL_TRUE);
369 }
370
371 void GLAPIENTRY
372 gluPickMatrix(GLdouble x, GLdouble y, GLdouble deltax, GLdouble deltay,
373 GLint viewport[4])
374 {
375 if (deltax <= 0 || deltay <= 0) {
376 return;
377 }
378
379 /* Translate and scale the picked region to the entire window */
380 glTranslatef((viewport[2] - 2 * (x - viewport[0])) / deltax,
381 (viewport[3] - 2 * (y - viewport[1])) / deltay, 0);
382 glScalef(viewport[2] / deltax, viewport[3] / deltay, 1.0);
383 }