OpenGL: 中平截头体六个裁剪平面方程的求解

OpenGL中平截头体六个裁剪平面方程的求解

• 空间中的坐标乘以视图模型矩阵得到视觉坐标, 再乘以投影矩阵得到裁剪坐标. 现在我求出裁剪矩阵M, 其等于投影矩阵左乘模型视图矩阵

 1     GLfloat   clip[16];
 2     GLfloat   proj[16];
 3     GLfloat   modl[16];
 4 
 5     glGetFloatv( GL_PROJECTION_MATRIX, proj );
 6     glGetFloatv( GL_MODELVIEW_MATRIX, modl );
 7      //  而后 proj 左乘 modl 得到裁剪矩阵
 8      clip[ 0] = modl[ 0] * proj[ 0] + modl[ 1] * proj[ 4] + modl[ 2] * proj[ 8] + modl[ 3] * proj[12];
 9     clip[ 1] = modl[ 0] * proj[ 1] + modl[ 1] * proj[ 5] + modl[ 2] * proj[ 9] + modl[ 3] * proj[13];
10     clip[ 2] = modl[ 0] * proj[ 2] + modl[ 1] * proj[ 6] + modl[ 2] * proj[10] + modl[ 3] * proj[14];
11     clip[ 3] = modl[ 0] * proj[ 3] + modl[ 1] * proj[ 7] + modl[ 2] * proj[11] + modl[ 3] * proj[15];
12 
13     clip[ 4] = modl[ 4] * proj[ 0] + modl[ 5] * proj[ 4] + modl[ 6] * proj[ 8] + modl[ 7] * proj[12];
14     clip[ 5] = modl[ 4] * proj[ 1] + modl[ 5] * proj[ 5] + modl[ 6] * proj[ 9] + modl[ 7] * proj[13];
15     clip[ 6] = modl[ 4] * proj[ 2] + modl[ 5] * proj[ 6] + modl[ 6] * proj[10] + modl[ 7] * proj[14];
16     clip[ 7] = modl[ 4] * proj[ 3] + modl[ 5] * proj[ 7] + modl[ 6] * proj[11] + modl[ 7] * proj[15];
17 
18     clip[ 8] = modl[ 8] * proj[ 0] + modl[ 9] * proj[ 4] + modl[10] * proj[ 8] + modl[11] * proj[12];
19     clip[ 9] = modl[ 8] * proj[ 1] + modl[ 9] * proj[ 5] + modl[10] * proj[ 9] + modl[11] * proj[13];
20     clip[10] = modl[ 8] * proj[ 2] + modl[ 9] * proj[ 6] + modl[10] * proj[10] + modl[11] * proj[14];
21     clip[11] = modl[ 8] * proj[ 3] + modl[ 9] * proj[ 7] + modl[10] * proj[11] + modl[11] * proj[15];
22 
23     clip[12] = modl[12] * proj[ 0] + modl[13] * proj[ 4] + modl[14] * proj[ 8] + modl[15] * proj[12];
24     clip[13] = modl[12] * proj[ 1] + modl[13] * proj[ 5] + modl[14] * proj[ 9] + modl[15] * proj[13];
25     clip[14] = modl[12] * proj[ 2] + modl[13] * proj[ 6] + modl[14] * proj[10] + modl[15] * proj[14];
26     clip[15] = modl[12] * proj[ 3] + modl[13] * proj[ 7] + modl[14] * proj[11] + modl[15] * proj[15];

• 假设空间有点P(x, y, z, 1), 乘以裁剪矩阵M得到点(x', y', z', w'). 齐次化为P' (x'/w', y'/w', z'/w', 1).假设点P'位于平截头体的近截面处, 则z'/w'=1.0. 如果位于平截头体的左平面处, 则x'/w'=-1.0, 以此类推.
• 假设裁剪矩阵M为上面代码的clip数组所表示. 则x' = clip[0] * x + clip[4] * y + clip[8] * z + clip[12]; w' = clip[3] * x + clip[7] * y + clip[11] * z + clip[15];
• 如果P'位于平截头体的右平面处, 则 x'/w' = 1.0, 即 x' - w' = 0.0. 所以得 (clip[0]-clip[3]) * x + (clip[4]-clip[7]) * y + (clip[8]-clip[11]) * z + (clip[12]-clip[15]) = 0; 其余平截头体的平面方程以此类推, 具体代码如下:

 1     GLfloat m_Frustum[6][4];
 2      // 提取右面的平面方程系数
 3      m_Frustum[0][0] = clip[ 3] - clip[ 0];
 4     m_Frustum[0][1] = clip[ 7] - clip[ 4];
 5     m_Frustum[0][2] = clip[11] - clip[ 8];
 6     m_Frustum[0][3] = clip[15] - clip[12];
 7     t = GLfloat(sqrt( m_Frustum[0][0] * m_Frustum[0][0] + m_Frustum[0][1] * m_Frustum[0][1] + m_Frustum[0][2] * m_Frustum[0][2] ));
 8     m_Frustum[0][0] /= t;
 9     m_Frustum[0][1] /= t;
10     m_Frustum[0][2] /= t;
11     m_Frustum[0][3] /= t;
12 
13      // 提取左面的平面方程系数
14      m_Frustum[1][0] = clip[ 3] + clip[ 0];
15     m_Frustum[1][1] = clip[ 7] + clip[ 4];
16     m_Frustum[1][2] = clip[11] + clip[ 8];
17     m_Frustum[1][3] = clip[15] + clip[12];
18     t = GLfloat(sqrt( m_Frustum[1][0] * m_Frustum[1][0] + m_Frustum[1][1] * m_Frustum[1][1] + m_Frustum[1][2] * m_Frustum[1][2] ));
19     m_Frustum[1][0] /= t;
20     m_Frustum[1][1] /= t;
21     m_Frustum[1][2] /= t;
22     m_Frustum[1][3] /= t;
23 
24      // 提取下面的平面方程系数
25      m_Frustum[2][0] = clip[ 3] + clip[ 1];
26     m_Frustum[2][1] = clip[ 7] + clip[ 5];
27     m_Frustum[2][2] = clip[11] + clip[ 9];
28     m_Frustum[2][3] = clip[15] + clip[13];
29     t = GLfloat(sqrt( m_Frustum[2][0] * m_Frustum[2][0] + m_Frustum[2][1] * m_Frustum[2][1] + m_Frustum[2][2] * m_Frustum[2][2] ));
30     m_Frustum[2][0] /= t;
31     m_Frustum[2][1] /= t;
32     m_Frustum[2][2] /= t;
33     m_Frustum[2][3] /= t;
34 
35      // 提取上面的平面方程系数
36      m_Frustum[3][0] = clip[ 3] - clip[ 1];
37     m_Frustum[3][1] = clip[ 7] - clip[ 5];
38     m_Frustum[3][2] = clip[11] - clip[ 9];
39     m_Frustum[3][3] = clip[15] - clip[13];
40     t = GLfloat(sqrt( m_Frustum[3][0] * m_Frustum[3][0] + m_Frustum[3][1] * m_Frustum[3][1] + m_Frustum[3][2] * m_Frustum[3][2] ));
41     m_Frustum[3][0] /= t;
42     m_Frustum[3][1] /= t;
43     m_Frustum[3][2] /= t;
44     m_Frustum[3][3] /= t;
45 
46      // 提取远面的平面方程系数
47      m_Frustum[4][0] = clip[ 3] - clip[ 2];
48     m_Frustum[4][1] = clip[ 7] - clip[ 6];
49     m_Frustum[4][2] = clip[11] - clip[10];
50     m_Frustum[4][3] = clip[15] - clip[14];
51     t = GLfloat(sqrt( m_Frustum[4][0] * m_Frustum[4][0] + m_Frustum[4][1] * m_Frustum[4][1] + m_Frustum[4][2] * m_Frustum[4][2] ));
52     m_Frustum[4][0] /= t;
53     m_Frustum[4][1] /= t;
54     m_Frustum[4][2] /= t;
55     m_Frustum[4][3] /= t;
56 
57      // 提取近面的平面方程系数
58      m_Frustum[5][0] = clip[ 3] + clip[ 2];
59     m_Frustum[5][1] = clip[ 7] + clip[ 6];
60     m_Frustum[5][2] = clip[11] + clip[10];
61     m_Frustum[5][3] = clip[15] + clip[14];
62     t = GLfloat(sqrt( m_Frustum[5][0] * m_Frustum[5][0] + m_Frustum[5][1] * m_Frustum[5][1] + m_Frustum[5][2] * m_Frustum[5][2] ));
63     m_Frustum[5][0] /= t;
64     m_Frustum[5][1] /= t;
65     m_Frustum[5][2] /= t;
66     m_Frustum[5][3] /= t;

• 上面代码求右平面方程时使用的是 w'-x' = 0;, 而不是 x'-w' = 0的方程, 这是由于满足方程 w' - x' < 0 的所有点位于右平面左方, 满足 w' + x' < 0 的所有点位于左平面左方. 这样求解一个点是否在平截头体之内, 可以只需要使用一个小于号就可以了, 代码如下:

1      //  假设有点 P, p.x, p.y, p.z 为其坐标
2       for(i = 0; i < 6; i++)
3     {
4          if(m_Frustum[i][0] * p.x + m_Frustum[i][1] * p.y + m_Frustum[i][2] * p.z + m_Frustum[i][3] <= 0)
5         {
6              return(FALSE);
7         }
8     }

http://www.cppblog.com/summericeyl/archive/2011/09/30/157215.html