用OpenGL实现射线拣取对象

第一步： 
实现屏幕坐标到三维世界空间坐标的转化，利用函数 gluUnProject直接可以得到屏幕坐标相应的三维空间坐标。

[cpp]  view plain copy

1. gluUnProject((GLdouble)xpos,(GLdouble)ypos,1.0,mvmatrix,projmatrix,viewport,&wx,&wy,&wz);    
2. //xpos,ypos 是以屏幕左下角为原点的屏幕坐标.   
3. //1.0代表返回zbuffer为1.0处（远剪切面交点）的世界坐标   
4. //mvmatrix 为视矩阵，通过GetDoublev(GL_MODELVIEW_MATRIX,mvmatrix)得到   
5. //projmatrix为投影矩阵，通过glGetDoublev(GL_PROJECTION_MATRIX,projmatrix)得到   
6. //viewport为视口，通过glGetIntegerv(GL_VIEWPORT,viewport)得到   
7. //wx、wy、wz 就是我们要得到的世界坐标，得到这样两个世界坐标，射线就确定了   
8. //或者也可以用原点（视点）来代替其中一个点，因为这条射线是从视点出发的。  

 

第二步： 
用射线和要检测的三角形求交点，用到的原理和公式如下。 
原理一：三角形内的任意一点都可以用变量u、v和其三个顶点坐标来确定

[cpp]  view plain copy

1. vPoint = V1 + u*(V2-V1) + v*(V3-V1) ;   
2. //0   
3. //V1,V2,V3为三角形的三个顶点，是已知量  

 

原理二：射线上的任意一点可以用射线的方向向量（格式化后的）乘以其模（该向量长度）来表示

  
  

   
   

    
    

     
     [cpp] 
     
     view plain
     
     copy
     
     

    
    
   
   
   
   
vPoint =originPoint+dir * len;  
  
  
射线和三角形相交则必定同时满足上面的两个条件所以有：

  
  

   
   

    
    

     
     [cpp] 
     
     view plain
     
     copy
     
     

    
    
   
   
   
   
(-Dir)*len+ (V2-V1)*u + (V3-V1)*v ＝ originPoint-V1;  
  
  
 

相当方程组: (len ,v ,u 为变量，其它为常量)

  
或

  
这是一个线性方程组，根据克拉姆法则，[-Dir，V2-V1，V3-V1]不为零。 
所以满足条件：00, ,0[-Dir，V2-V1，V3-V1]写成矩阵形式为：

 

伪码实现为:

  
  

   
   

    
    

     
     [cpp] 
     
     view plain
     
     copy
     
     

    
    
   
   
   
   
// 三角形两个边的向量  
VECTOR3 edge1 = v1 - v0;  
VCTOR3  edge2 = v2 - v0;  
VCTOR3  pvec;  
VEC3Cross( &pvec, &dir, &edge2 );    // 差积  
FLOAT det = VEC3Dot( &edge1, &pvec );// 点积  
                                     // det其含义为[－Dir，V2－V1，V3－V1]矩阵展开  
VECTOR3 tvec;  
if( det > 0 )   
{  
    tvec = orig - v0;                // 从正面穿越三角形，三角形和视点相对的面为正面  
}  
else  
{  
    tvec = v0 - orig;                // 反面穿越三角形穿越三角形  
     det = -det;  
}  
  
if( det < 0.0001f )    return FALSE;  
  
// 求u的值，求线性方程组的解展开后等同于求点积展开  
 *u = VEC3Dot( &tvec, &pvec );  
if( *u < 0.0f || *u > det ) return FALSE;  
  
// 求v的值  
VECTOR3 qvec;  
VEC3Cross( &qvec, &tvec, &edge1 );  
*v = VEC3Dot( &dir, &qvec );  
 if( *v < 0.0f || *u + *v > det ) return FALSE;  
  
// 计算t,并把t,u,v放缩为合法值  
*t = D3DXVec3Dot( &edge2, &qvec );  
// 前面的t,v,u在计算时多乘了一个系数det  
  
FLOAT fInvDet = 1.0f / det;  
*t *= fInvDet;  
*u *= fInvDet;  
*v *= fInvDet;  
  
  
 

这个算法是微软给出的，从几何角度分析其含义十分难懂，真正的方法是根据线性方程租求解，巧的是文中的方法恰好和线性方程组整理出来的东西相符合，这大概就是几何和代数相通的原理。

  
  

   
   

    
    

     
     [cpp] 
     
     view plain
     
     copy
     
     

    
    
   
   
   
   
bool IntersectTriangle()  
{  
    GLfloat edge1[3];  
    GLfloat edge2[3];  
  
    edge1[0]=V1[0]-V0[0];  
    edge1[1]=V1[1]-V0[1];  
    edge1[2]=V1[2]-V0[2];  
  
    edge2[0]=V2[0]-V0[0];  
    edge2[1]=V2[1]-V0[1];  
    edge2[2]=V2[2]-V0[2];  
  
    GLfloat dir[3];  
    dir[0]=g_farxyz[0]-g_nearxyz[0];  
    dir[1]=g_farxyz[1]-g_nearxyz[1];  
    dir[2]=g_farxyz[2]-g_nearxyz[2];  
  
    GLfloat w = (GLfloat)sqrt((double)pow(dir[0],2.0)+(double)pow(dir[1],2.0)+(double)pow(dir[2],2.0));  
    dir[0] /= w;  
    dir[1] /= w;  
    dir[2] /= w;  
  
    GLfloat pvec[3];  
    pvec[0]= dir[1]*edge2[2] - dir[2]*edge2[1];  
    pvec[1]= dir[2]*edge2[0] - dir[0]*edge2[2];  
    pvec[2]= dir[0]*edge2[1] - dir[1]*edge2[0];  
  
    GLfloat det ;  
    det = edge1[0]*pvec[0]+edge1[1]*pvec[1]+edge1[2]*pvec[2];  
  
    GLfloat tvec[3];  
    if( det > 0 )  
   {  
     tvec[0] = g_nearxyz[0] - V0[0];  
     tvec[1] = g_nearxyz[1] - V0[1];  
     tvec[2] = g_nearxyz[2] - V0[2];  
   }  
   else  
   {  
     tvec[0] = V0[0] - g_nearxyz[0];  
     tvec[1] = V0[1] - g_nearxyz[1];  
     tvec[2] = V0[2] - g_nearxyz[2];  
     det = -det ;  
   }  
  
   if( det < 0.0001f )    return false;  
  
   GLfloat u ;  
   u = tvec[0]*pvec[0]+ tvec[1]*pvec[1]+ tvec[2]*pvec[2];  
  
   if( u < 0.0f || u > det )   return false;  
  
   GLfloat qvec[3];  
   qvec[0]= tvec[1]*edge1[2] - tvec[2]*edge1[1];  
   qvec[1]= tvec[2]*edge1[0] - tvec[0]*edge1[2];  
   qvec[2]= tvec[0]*edge1[1] - tvec[1]*edge1[0];  
  
   GLfloat v;  
   v = dir[0]*qvec[0]+dir[1]*qvec[1]+dir[2]*qvec[2];  
   if( v < 0.0f || u + v > det )      return false;  
  
   GLfloat t = edge2[0]*qvec[0]+edge2[1]*qvec[1]+edge2[2]*qvec[2];  
   GLfloat fInvDet = 1.0f / det;  
   t *= fInvDet;  
   u *= fInvDet;  
   v *= fInvDet;  
   return true;  
}  
  
void pick(GLfloat xpos,GLfloat ypos)  
{  
    xpos,ypos;  
    GLint viewport[4];  
    GLdouble mvmatrix[16],projmatrix[16];  
    GLint realy;  
    GLdouble wx,wy,wz;  
  
    glGetIntegerv(GL_VIEWPORT,viewport);  
    glGetDoublev(GL_MODELVIEW_MATRIX,mvmatrix);  
    glGetDoublev(GL_PROJECTION_MATRIX,projmatrix);  
  
    realy = viewport[3]-(GLint)ypos -1;// 左下角为坐标原点  
     gluUnProject((GLdouble)xpos,(GLdouble)realy,0.0,mvmatrix,projmatrix,viewport,&wx,&wy,&wz);  
  
    g_nearxyz[0] = (GLfloat)wx;  
    g_nearxyz[1] = (GLfloat)wy;  
    g_nearxyz[2] = (GLfloat)wz;  
  
    gluUnProject((GLdouble)xpos,(GLdouble)realy,1.0,mvmatrix,projmatrix,viewport,&wx,&wy,&wz);  
  
    g_farxyz[0] = (GLfloat)wx;  
    g_farxyz[1] = (GLfloat)wy;  
    g_farxyz[2] = (GLfloat)wz;  
  
    g_color = 0.0;  
    if(IntersectTriangle()) g_color=1.0;  
}  
GLfloat V0[3]={1.0,0.0,-1.0 };  
GLfloat V1[3]={0.0,1.0,-1.0 };  
GLfloat V2[3]={0.0,0.0,-2.0 };  
Void DrawGLScene(GLvoid)   
{  
    glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);   
    glBegin(GL_TRIANGLES);  
    glColor3f(g_color,0.0,1.0);  
    glVertex3fv(V0);// 如果加了glTranslatef之类的变换函数，射线应该反向变化  
     glVertex3fv(V1);  
    glVertex3fv(V2);  
    glEnd();  
    SwapBuffers(hDC);   
}