用OpenGL實現射線揀取物件程式設計(轉)

用OpenGL實現射線揀取物件程式設計(轉)[@more@]

　　關於用射線原理來揀取物件網上已經有完整的理論，另外DirectX也提供了一個Pick例子來演示，在這裡我將這些資料和理論來稍微的總結，並給出OpenGL下的完整實現。

　　相關的理論大體來自一篇英文資料和一篇總結性的中文資料，分別是：

　　

　　

　　前一篇完整講述用DirectX實現射線揀取物體的原理和實現。後一篇講述的是二維螢幕空間到三維世界空間的轉換原理。前一篇的名字是“Direct3D中實現圖元的滑鼠拾取”，它的講述很好也很透徹。後一篇講述射線形成的原理，並且有原始碼例子。

　　下面就OpenGL進行實現。

　　第一步：

　　實現螢幕座標到三維世界空間座標的轉化，在這一步Opengl要比DirectX簡單的多，利用函式 gluUnProject直接可以得到螢幕座標相應的三維空間座標，示例如下：

　　gluUnProject((GLdouble)xpos,(GLdouble)ypos,1.0,mvmatrix,projmatrix,viewport,&wx,&wy,&wz);　xpos 和ypos 是以螢幕左下角為原點的螢幕座標，1.0代表返回zbuffer為1.0處（遠剪下面交點）的世界座標，mvmatrix 為視矩陣，透過GetDoublev(GL_MODELVIEW_MATRIX,mvmatrix)得到，projmatrix為投影矩陣，透過glGetDoublev(GL_PROJECTION_MATRIX,projmatrix)得到，viewport為視口，透過glGetIntegerv(GL_VIEWPORT,viewport)得到，剩下的wx、wy、wz 就是我們要得到的世界座標，得到這樣兩個世界座標，射線就確定了，或者也可以用原點（視點）來代替其中一個點，因為這條射線是從視點出發的。

　　第二步：

　　用射線和要檢測的三角形求交點，用到的原理和公式如下。

　　原理一：三角形內的任意一點都可以用變數u、v和其三個頂點座標來確定，其中0　　原理二：射線上的任意一點可以用射線的方向向量（格式化後的）乘以其模（該向量長度）來表示，記為:vPoint =originPoint+dir * len

　　如果和三角形相交則必定同時滿足上面的兩個條件所以有：

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

　　相當方程組: (len ,v ,u 為變數，其它為常量)

　　（-Dir.x）*len +(V2.x-V1.x)*u + (V3.x ? V1.x )*v = originPoint.x －V1.x

　　（-Dir.y）*len +(V2.y-V1.y)*u + (V3.y ? V1.y )*v = originPoint.y －V1.y

　　（-Dir.z）*len +(V2.z-V1.z)*u + (V3.z ? V1.z )*v = originPoint.z －V1.z

　　或：

　　len

　　【－Dir，V2－V1，V3－V1】{　u　 } = originPoint ? V1

　　v

　　這是一個線性方程組，根據克拉姆法則，【－Dir，V2－V1，V3－V1】不為零。

　　所以滿足條件：00, ,0　　【－Dir，V2－V1，V3－V1】寫成矩陣形式為：

　　| -Dir.x , V2.x-V1.x , V3.x ? V1.x |

　　| -Dir.y , V2.y-V1.y , V3.y ? V1.x |

　　| -Dir.z , V2.z-V1.z , V3.z ? V1.z |

　　偽碼實現（原理在DirectX Pick例子中有原始碼實現）：

　　// 三角形兩個邊的向量

　　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 )// 接近零視為0

　　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;

　　// 這裡這個演算法是微軟給出的，從幾何角度分析其含義十分難懂，真正的方法是根據線性方程租求解，巧的是文中的方法恰好和線性方程組整理出來的東西相符合，這大概就是幾何和代數相通的原理。

　　原始碼（VC6.0 + OPENGL + WINDOWS2000，除錯透過）：

　　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);

　　}