home *** CD-ROM | disk | FTP | other *** search
/ Magazyn Amiga 5 / MA_Cover_5.iso / ppc / mesa / src-glu / project.c < prev    next >
Encoding:
C/C++ Source or Header  |  1998-01-31  |  10.8 KB  |  408 lines
  1. /* $Id: project.c,v 1.5 1997/07/24 01:28:44 brianp Exp $ */
  2.  
  3. /*
  4.  * Mesa 3-D graphics library
  5.  * Version:  2.4
  6.  * Copyright (C) 1995-1997  Brian Paul
  7.  *
  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.
  12.  *
  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.
  17.  *
  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.
  21.  */
  22.  
  23.  
  24. /*
  25.  * $Log: project.c,v $
  26.  * Revision 1.5  1997/07/24 01:28:44  brianp
  27.  * changed precompiled header symbol from PCH to PC_HEADER
  28.  *
  29.  * Revision 1.4  1997/05/28 02:29:38  brianp
  30.  * added support for precompiled headers (PCH), inserted APIENTRY keyword
  31.  *
  32.  * Revision 1.3  1997/04/11 23:22:42  brianp
  33.  * added divide by zero checks to gluProject() and gluUnproject()
  34.  *
  35.  * Revision 1.2  1997/01/29 19:05:29  brianp
  36.  * faster invert_matrix() function from Stephane Rehel
  37.  *
  38.  * Revision 1.1  1996/09/27 01:19:39  brianp
  39.  * Initial revision
  40.  *
  41.  */
  42.  
  43.  
  44. #ifdef PC_HEADER
  45. #include "all.h"
  46. #else
  47. #include <stdio.h>
  48. #include <string.h>
  49. #include <math.h>
  50. #include "gluP.h"
  51. #endif
  52.  
  53.  
  54. /*
  55.  * This code was contributed by Marc Buffat (buffat@mecaflu.ec-lyon.fr).
  56.  * Thanks Marc!!!
  57.  */
  58.  
  59.  
  60.  
  61. /* implementation de gluProject et gluUnproject */
  62. /* M. Buffat 17/2/95 */
  63.  
  64.  
  65.  
  66. /*
  67.  * Transform a point (column vector) by a 4x4 matrix.  I.e.  out = m * in
  68.  * Input:  m - the 4x4 matrix
  69.  *         in - the 4x1 vector
  70.  * Output:  out - the resulting 4x1 vector.
  71.  */
  72. static void transform_point( GLdouble out[4], const GLdouble m[16],
  73.                  const GLdouble in[4] )
  74. {
  75. #define M(row,col)  m[col*4+row]
  76.    out[0] = M(0,0) * in[0] + M(0,1) * in[1] + M(0,2) * in[2] + M(0,3) * in[3];
  77.    out[1] = M(1,0) * in[0] + M(1,1) * in[1] + M(1,2) * in[2] + M(1,3) * in[3];
  78.    out[2] = M(2,0) * in[0] + M(2,1) * in[1] + M(2,2) * in[2] + M(2,3) * in[3];
  79.    out[3] = M(3,0) * in[0] + M(3,1) * in[1] + M(3,2) * in[2] + M(3,3) * in[3];
  80. #undef M
  81. }
  82.  
  83.  
  84.  
  85.  
  86. /*
  87.  * Perform a 4x4 matrix multiplication  (product = a x b).
  88.  * Input:  a, b - matrices to multiply
  89.  * Output:  product - product of a and b
  90.  */
  91. static void matmul( GLdouble *product, const GLdouble *a, const GLdouble *b )
  92. {
  93.    /* This matmul was contributed by Thomas Malik */
  94.    GLdouble temp[16];
  95.    GLint i;
  96.  
  97. #define A(row,col)  a[(col<<2)+row]
  98. #define B(row,col)  b[(col<<2)+row]
  99. #define T(row,col)  temp[(col<<2)+row]
  100.  
  101.    /* i-te Zeile */
  102.    for (i = 0; i < 4; i++)
  103.      {
  104.     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);
  105.     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);
  106.     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);
  107.     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);
  108.      }
  109.  
  110. #undef A
  111. #undef B
  112. #undef T
  113.    MEMCPY( product, temp, 16*sizeof(GLdouble) );
  114. }
  115.  
  116.  
  117. static GLdouble Identity[16] = {
  118.    1.0, 0.0, 0.0, 0.0,
  119.    0.0, 1.0, 0.0, 0.0,
  120.    0.0, 0.0, 1.0, 0.0,
  121.    0.0, 0.0, 0.0, 1.0
  122. };
  123.  
  124.  
  125.  
  126. /*
  127.  * Compute the inverse of a 4x4 matrix.  Contributed by scotter@lafn.org
  128.  */
  129. static void invert_matrix_general( const GLdouble *m, GLdouble *out )
  130. {
  131.  
  132. /* NB. OpenGL Matrices are COLUMN major. */
  133. #define MAT(m,r,c) (m)[(c)*4+(r)]
  134.  
  135. /* Here's some shorthand converting standard (row,column) to index. */
  136. #define m11 MAT(m,0,0)
  137. #define m12 MAT(m,0,1)
  138. #define m13 MAT(m,0,2)
  139. #define m14 MAT(m,0,3)
  140. #define m21 MAT(m,1,0)
  141. #define m22 MAT(m,1,1)
  142. #define m23 MAT(m,1,2)
  143. #define m24 MAT(m,1,3)
  144. #define m31 MAT(m,2,0)
  145. #define m32 MAT(m,2,1)
  146. #define m33 MAT(m,2,2)
  147. #define m34 MAT(m,2,3)
  148. #define m41 MAT(m,3,0)
  149. #define m42 MAT(m,3,1)
  150. #define m43 MAT(m,3,2)
  151. #define m44 MAT(m,3,3)
  152.  
  153.    GLdouble det;
  154.    GLdouble d12, d13, d23, d24, d34, d41;
  155.    GLdouble tmp[16]; /* Allow out == in. */
  156.  
  157.    /* Inverse = adjoint / det. (See linear algebra texts.)*/
  158.  
  159.    /* pre-compute 2x2 dets for last two rows when computing */
  160.    /* cofactors of first two rows. */
  161.    d12 = (m31*m42-m41*m32);
  162.    d13 = (m31*m43-m41*m33);
  163.    d23 = (m32*m43-m42*m33);
  164.    d24 = (m32*m44-m42*m34);
  165.    d34 = (m33*m44-m43*m34);
  166.    d41 = (m34*m41-m44*m31);
  167.  
  168.    tmp[0] =  (m22 * d34 - m23 * d24 + m24 * d23);
  169.    tmp[1] = -(m21 * d34 + m23 * d41 + m24 * d13);
  170.    tmp[2] =  (m21 * d24 + m22 * d41 + m24 * d12);
  171.    tmp[3] = -(m21 * d23 - m22 * d13 + m23 * d12);
  172.  
  173.    /* Compute determinant as early as possible using these cofactors. */
  174.    det = m11 * tmp[0] + m12 * tmp[1] + m13 * tmp[2] + m14 * tmp[3];
  175.  
  176.    /* Run singularity test. */
  177.    if (det == 0.0) {
  178.       /* printf("invert_matrix: Warning: Singular matrix.\n"); */
  179.       MEMCPY( out, Identity, 16*sizeof(GLdouble) );
  180.    }
  181.    else {
  182.       GLdouble invDet = 1.0 / det;
  183.       /* Compute rest of inverse. */
  184.       tmp[0] *= invDet;
  185.       tmp[1] *= invDet;
  186.       tmp[2] *= invDet;
  187.       tmp[3] *= invDet;
  188.  
  189.       tmp[4] = -(m12 * d34 - m13 * d24 + m14 * d23) * invDet;
  190.       tmp[5] =  (m11 * d34 + m13 * d41 + m14 * d13) * invDet;
  191.       tmp[6] = -(m11 * d24 + m12 * d41 + m14 * d12) * invDet;
  192.       tmp[7] =  (m11 * d23 - m12 * d13 + m13 * d12) * invDet;
  193.  
  194.       /* Pre-compute 2x2 dets for first two rows when computing */
  195.       /* cofactors of last two rows. */
  196.       d12 = m11*m22-m21*m12;
  197.       d13 = m11*m23-m21*m13;
  198.       d23 = m12*m23-m22*m13;
  199.       d24 = m12*m24-m22*m14;
  200.       d34 = m13*m24-m23*m14;
  201.       d41 = m14*m21-m24*m11;
  202.  
  203.       tmp[8] =  (m42 * d34 - m43 * d24 + m44 * d23) * invDet;
  204.       tmp[9] = -(m41 * d34 + m43 * d41 + m44 * d13) * invDet;
  205.       tmp[10] =  (m41 * d24 + m42 * d41 + m44 * d12) * invDet;
  206.       tmp[11] = -(m41 * d23 - m42 * d13 + m43 * d12) * invDet;
  207.       tmp[12] = -(m32 * d34 - m33 * d24 + m34 * d23) * invDet;
  208.       tmp[13] =  (m31 * d34 + m33 * d41 + m34 * d13) * invDet;
  209.       tmp[14] = -(m31 * d24 + m32 * d41 + m34 * d12) * invDet;
  210.       tmp[15] =  (m31 * d23 - m32 * d13 + m33 * d12) * invDet;
  211.  
  212.       MEMCPY(out, tmp, 16*sizeof(GLdouble));
  213.    }
  214.  
  215. #undef m11
  216. #undef m12
  217. #undef m13
  218. #undef m14
  219. #undef m21
  220. #undef m22
  221. #undef m23
  222. #undef m24
  223. #undef m31
  224. #undef m32
  225. #undef m33
  226. #undef m34
  227. #undef m41
  228. #undef m42
  229. #undef m43
  230. #undef m44
  231. #undef MAT
  232. }
  233.  
  234.  
  235. /*
  236.  * Invert matrix m.  This algorithm contributed by Stephane Rehel
  237.  * <rehel@worldnet.fr>
  238.  */
  239. static void invert_matrix( const GLdouble *m, GLdouble *out )
  240. {
  241. /* NB. OpenGL Matrices are COLUMN major. */
  242. #define MAT(m,r,c) (m)[(c)*4+(r)]
  243.  
  244. /* Here's some shorthand converting standard (row,column) to index. */
  245. #define m11 MAT(m,0,0)
  246. #define m12 MAT(m,0,1)
  247. #define m13 MAT(m,0,2)
  248. #define m14 MAT(m,0,3)
  249. #define m21 MAT(m,1,0)
  250. #define m22 MAT(m,1,1)
  251. #define m23 MAT(m,1,2)
  252. #define m24 MAT(m,1,3)
  253. #define m31 MAT(m,2,0)
  254. #define m32 MAT(m,2,1)
  255. #define m33 MAT(m,2,2)
  256. #define m34 MAT(m,2,3)
  257. #define m41 MAT(m,3,0)
  258. #define m42 MAT(m,3,1)
  259. #define m43 MAT(m,3,2)
  260. #define m44 MAT(m,3,3)
  261.  
  262.    register GLdouble det;
  263.    GLdouble tmp[16]; /* Allow out == in. */
  264.  
  265.    if( m41 != 0. || m42 != 0. || m43 != 0. || m44 != 1. ) {
  266.       invert_matrix_general(m, out);
  267.       return;
  268.    }
  269.  
  270.    /* Inverse = adjoint / det. (See linear algebra texts.)*/
  271.  
  272.    tmp[0]= m22 * m33 - m23 * m32;
  273.    tmp[1]= m23 * m31 - m21 * m33;
  274.    tmp[2]= m21 * m32 - m22 * m31;
  275.  
  276.    /* Compute determinant as early as possible using these cofactors. */
  277.    det= m11 * tmp[0] + m12 * tmp[1] + m13 * tmp[2];
  278.  
  279.    /* Run singularity test. */
  280.    if (det == 0.0) {
  281.       /* printf("invert_matrix: Warning: Singular matrix.\n"); */
  282.       MEMCPY( out, Identity, 16*sizeof(GLdouble) );
  283.    }
  284.    else {
  285.       GLdouble d12, d13, d23, d24, d34, d41;
  286.       register GLdouble im11, im12, im13, im14;
  287.  
  288.       det= 1. / det;
  289.  
  290.       /* Compute rest of inverse. */
  291.       tmp[0] *= det;
  292.       tmp[1] *= det;
  293.       tmp[2] *= det;
  294.       tmp[3]  = 0.;
  295.  
  296.       im11= m11 * det;
  297.       im12= m12 * det;
  298.       im13= m13 * det;
  299.       im14= m14 * det;
  300.       tmp[4] = im13 * m32 - im12 * m33;
  301.       tmp[5] = im11 * m33 - im13 * m31;
  302.       tmp[6] = im12 * m31 - im11 * m32;
  303.       tmp[7] = 0.;
  304.  
  305.       /* Pre-compute 2x2 dets for first two rows when computing */
  306.       /* cofactors of last two rows. */
  307.       d12 = im11*m22 - m21*im12;
  308.       d13 = im11*m23 - m21*im13;
  309.       d23 = im12*m23 - m22*im13;
  310.       d24 = im12*m24 - m22*im14;
  311.       d34 = im13*m24 - m23*im14;
  312.       d41 = im14*m21 - m24*im11;
  313.  
  314.       tmp[8] =  d23;
  315.       tmp[9] = -d13;
  316.       tmp[10] = d12;
  317.       tmp[11] = 0.;
  318.  
  319.       tmp[12] = -(m32 * d34 - m33 * d24 + m34 * d23);
  320.       tmp[13] =  (m31 * d34 + m33 * d41 + m34 * d13);
  321.       tmp[14] = -(m31 * d24 + m32 * d41 + m34 * d12);
  322.       tmp[15] =  1.;
  323.  
  324.       MEMCPY(out, tmp, 16*sizeof(GLdouble));
  325.   }
  326.  
  327. #undef m11
  328. #undef m12
  329. #undef m13
  330. #undef m14
  331. #undef m21
  332. #undef m22
  333. #undef m23
  334. #undef m24
  335. #undef m31
  336. #undef m32
  337. #undef m33
  338. #undef m34
  339. #undef m41
  340. #undef m42
  341. #undef m43
  342. #undef m44
  343. #undef MAT
  344. }
  345.  
  346.  
  347.  
  348. /* projection du point (objx,objy,obz) sur l'ecran (winx,winy,winz) */
  349. GLint APIENTRY gluProject(GLdouble objx,GLdouble objy,GLdouble objz,
  350.                           const GLdouble model[16],const GLdouble proj[16],
  351.                           const GLint viewport[4],
  352.                           GLdouble *winx,GLdouble *winy,GLdouble *winz)
  353. {
  354.     /* matrice de transformation */
  355.     GLdouble in[4],out[4];
  356.  
  357.     /* initilise la matrice et le vecteur a transformer */
  358.     in[0]=objx; in[1]=objy; in[2]=objz; in[3]=1.0;
  359.     transform_point(out,model,in);
  360.     transform_point(in,proj,out);
  361.  
  362.     /* d'ou le resultat normalise entre -1 et 1*/
  363.     if (in[3]==0.0)
  364.        return GL_FALSE;
  365.  
  366.     in[0]/=in[3]; in[1]/=in[3]; in[2]/=in[3];
  367.  
  368.     /* en coordonnees ecran */
  369.     *winx = viewport[0]+(1+in[0])*viewport[2]/2;
  370.     *winy = viewport[1]+(1+in[1])*viewport[3]/2;
  371.     /* entre 0 et 1 suivant z */
  372.     *winz = (1+in[2])/2;
  373.     return GL_TRUE;
  374. }
  375.  
  376.  
  377.  
  378. /* transformation du point ecran (winx,winy,winz) en point objet */
  379. GLint APIENTRY gluUnProject(GLdouble winx,GLdouble winy,GLdouble winz,
  380.                             const GLdouble model[16],const GLdouble proj[16],
  381.                             const GLint viewport[4],
  382.                             GLdouble *objx,GLdouble *objy,GLdouble *objz)
  383. {
  384.     /* matrice de transformation */
  385.     GLdouble m[16], A[16];
  386.     GLdouble in[4],out[4];
  387.  
  388.     /* transformation coordonnees normalisees entre -1 et 1 */
  389.     in[0]=(winx-viewport[0])*2/viewport[2] - 1.0;
  390.     in[1]=(winy-viewport[1])*2/viewport[3] - 1.0;
  391.     in[2]=2*winz - 1.0;
  392.     in[3]=1.0;
  393.  
  394.     /* calcul transformation inverse */
  395.     matmul(A,proj,model);
  396.     invert_matrix(A,m);
  397.  
  398.     /* d'ou les coordonnees objets */
  399.     transform_point(out,m,in);
  400.     if (out[3]==0.0)
  401.        return GL_FALSE;
  402.     *objx=out[0]/out[3];
  403.     *objy=out[1]/out[3];
  404.     *objz=out[2]/out[3];
  405.     return GL_TRUE;
  406. }
  407.  
  408.