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