Fast Extraction of Viewing Frustum Planes from the WorldView-Projection Matrix

Fast Extraction of Viewing Frustum Planes from the WorldView-Projection Matrix Authors (in alphabetical order): Gil Gribb ( ggribb@ravensoft.com ) Klaus Hartmann ( k_hartmann@osnabrueck.netsurf.de ) 1. Introduction We present an algorithm for extracting the six viewing frustum planes directly from the combined world, view, and projection matrices. The algorithm is fast, accurate, and general. It allows one to quickly determine the frustum planes in camera space, world space, or object space. The paper is divided into two main chapters. The first chapter "Plane Extraction in Direct3D" shows how to extract the viewing frustum planes in Direct3D, and the second chapter "Plane Extraction in OpenGL" shows how to do the same in OpenGL. In addition, there are two appendices. Appendix A reviews some of the maths behind plane equations, and Appendix B contains sample implementations of the plane extraction algorithm for both, Direct3D and OpenGL. We hope you enjoy this short paper. Gil Gribb Klaus Hartmann 06/15/2001 2. Plane Extraction in Direct3D We start off by extracting the frustum planes from the projection matrix only. That is, the world and view matrices are both identity matrices. This means that the camera is located at the origin of the world coordinate system, and we are looking along the positive z-axis. Let v = = ( ) xyzw 1 be a vertex and let ( ) M = mij be a 4 4 × projection matrix. Transforming the vertex v with the matrix M results in the transformed vertex v ' ''' ' = ( x yzw ). This can be written as: ( ) 11 21 31 41 1 12 22 32 42 2 13 23 33 43 3 14 24 34 44 4 ''' ' T T xm ym zm wm col xm ym zm wm col xyzw xm ym zm wm col xm ym zm wm col   ⋅ +⋅ +⋅ + ⋅   •     ⋅ +⋅ +⋅ + ⋅ •   =⇒ = = ⋅ +⋅ +⋅ + ⋅   •   ⋅ +⋅ +⋅ + ⋅ •     v v v' vM v v where • denotes the dot product, and col m m m m j jjjj = ( 1234 ) denotes the vector represented by the th j column of matrix M . For example, 2 v • col denotes the dot product of v with the vector represented by the 2nd column of matrix M . After this transformation, the vertex v ' is in homogeneous clipping space. In this space the viewing frustum actually is an axis-aligned box whose size is API-specific. If vertex v ' is inside this box, then the untransformed vertex v is inside the 'untransformed' viewing frustum. Using Direct3D, the vertex v ' is inside the box, if the following inequalities are all true for the components of v ' : '' ' ''' 0' ' wxw wyw z w − < < − < < < < There are a couple conclusions we can draw from these inequalities. They are listed in the following table. If − < w x ' ' then x ' is in the inside-halfspace of the left clipping plane If x ' ' < w then x ' is in the inside-halfspace of the right clipping plane If − < w y ' ' then y ' is in the inside-halfspace of the bottom clipping plane If y w ' ' < then y ' is in the inside-halfspace of the top clipping plane If 0 ' < z then z' is in the inside-halfspace of the near clipping plane If z w ' ' < then z' is in the inside-halfspace of the far clipping plane Now suppose that we wanted to test, if x ' is in the inside-halfspace of the left clipping plane. This is the case, if the following inequality is true: −w x ' ' < Using the information from the beginning of this chapter, we can rewrite this inequality as: 4 1 −( )( ) v v • <• col col This in turn is equal to: 4 1 0( )( ) < v v • +• col col And finally: 4 1 0( ) < v • + col col This already represents the plane equation for the left clipping plane of the 'untransformed' viewing frustum: 14 11 24 21 34 31 44 41 xm m ym m zm m wm m ( ) ( ) ( ) ( )0 ++ ++ ++ + = Since w =1, we can simplify this as follows: 14 11 24 21 34 31 44 41 xm m ym m zm m m m ( ) ( ) ( )( )0 ++ ++ ++ + = This gives a standard plane equation: ax by cz d + + += 0 , where 14 11 24 21 34 31 44 41 am m bm m cm m dm m =+ =+ =+ =+ , ,, We have now shown that the left clipping plane of the viewing frustum can be directly extracted from the projection matrix. It is important to note, that the resulting plane equation is not normalized (i.e., the plane’s normal vector is not a unit vector), and that the normal vector is pointing to the inside halfspace. This means, that 0 < +++ ax by cz d , if the vertex v is in the inside halfspace of the left clipping plane. We can now repeat the above steps for the other inequalities to determine the plane equations for the remaining clipping planes. Clipping plane Inequality Coefficients for the plane equation left 4 1 1 1 4 4 ' ' ( ) 0 ( )( ) 0 ( ) w x col col col col col col − < −• < • < • +• < • + v v v v v 14 11 24 21 34 31 44 41 am m bm m cm m dm m = + = + = + = + right 1 4 4 4 1 1 ' ' 0 ( )( ( ) ) 0 x w col co col col col co l l < •< • < • −• < • − v v v v v 14 11 24 21 34 31 44 41 am m bm m cm m dm m = − = − = − = − bottom 4 2 2 2 4 4 ' ' ( ) 0 ( )( ) 0 ( ) w y col col col col col col − < −• < • < • +• < • + v v v v v 14 12 24 22 34 32 44 42 am m bm m cm m dm m = + = + = + = + top 2 4 4 4 2 2 ' ' 0 ( )( ( ) ) 0 y w col co col col col co l l < •< • < • −• < • − v v v v v 14 12 24 22 34 32 44 42 am m bm m cm m dm m = − = − = − = − near 3 0 ' 0 o z c l < < • v 13 23 33 43 a m b m c m d m = = = = far 3 4 4 4 3 3 ' ' 0 ( )( ( ) ) 0 z w col co col col col co l l < •< • < • −• < • − v v v v v 14 13 24 23 34 33 44 43 am m bm m cm m dm m = − = − = − = − 2.1 Supporting Non-Identity World and View Matrices Until now we have assumed that both, the world and the view matrix are identity matrices. However, the goal obviously is to make the algorithm work for an arbitrary view. This is, in fact, so easy that it’s almost unbelievable. If you think about it for a moment then you’ll immediately understand it, and that’s why we are not going to explain this in any great detail. All we are giving to you is this: 1. If the matrixM is equal to the projection matrix P (i.e., M P = ), then the algorithm gives the clipping planes in view space (i.e., camera space). 2. If the matrixM is equal to the combined view and projection matrices, then the algorithm gives the clipping planes in world space (i.e., M VP = ⋅ , where V is the view matrix, and P is the projection matrix). 3. If the matrix M is equal to the combined world, view, and projection matrices, then the algorithm gives the clipping planes in object space (i.e., M WVP = ⋅ ⋅ , where W is the world matrix, V is the view matrix, and P is the projection matrix). 4. and so on... 3. Plane Extraction in OpenGL We start off by extracting the frustum planes from the projection matrix only. That is, the modelview matrix is an identity matrix. This means that the camera is located at the origin of the world coordinate system, and we are looking along the positive z-axis. Let ( ) 1 T v = = xyzw be a vertex and let ( ) M = mij be a 4 4 × projection matrix. Transforming the vertex v with the matrix M results in the transformed vertex ' ''' ' ( ) T v = x yzw . This can be written as: 11 12 13 14 1 21 22 23 24 2 31 32 33 34 3 41 42 43 44 4 ' ' ' ' x x m y m z m w m row y x m y m z m w m row z x m y m z m w m row w xm ym zm wm row     ⋅ +⋅ +⋅ + ⋅ •         ⋅ +⋅ +⋅ + ⋅ •     =⇒ = =   ⋅ +⋅ +⋅ + ⋅ •         ⋅ +⋅ +⋅ + ⋅ •       v v Mv v' v v where • denotes the dot product, and row m m m m i ii ii = ( 1234 ) denotes the vector represented by the th i row of matrix M . For example, 2 v • row denotes the dot product of v with the vector represented by the 2nd row of matrix M . After this transformation, the vertex v ' is in homogeneous clipping space. In this space the viewing frustum actually is an axis-aligned box whose size is API-specific. If vertex v ' is inside this box, then the untransformed vertex v is inside the 'untransformed' viewing frustum. Using OpenGL, the vertex v ' is inside the box, if the following inequalities are all true for the components of v ' : '' ' ''' '' ' wxw wyw wzw − < < − < < − < < There are a couple conclusions we can draw from these inequalities. They are listed in the following table. If − < w x ' ' then x ' is in the inside-halfspace of the left clipping plane If x ' ' < w then x ' is in the inside-halfspace of the right clipping plane If − < w y ' ' then y ' is in the inside-halfspace of the bottom clipping plane If y w ' ' < then y ' is in the inside-halfspace of the top clipping plane If − < w z ' ' then z' is in the inside-halfspace of the near clipping plane If z w ' ' < then z' is in the inside-halfspace of the far clipping plane Now suppose that we wanted to test, if x ' is in the inside-halfspace of the left clipping plane. This is the case, if the following inequality is true: −w x ' ' < Using the information from the beginning of this chapter, we can rewrite this inequality as: 4 1 −( )( ) v v • <• row row This in turn is equal to: 4 1 0( )( ) < v v • +• row row And finally: 4 1 0( ) < v • + row row This already represents the plane equation for the left clipping plane of the 'untransformed' viewing frustum: 41 11 42 12 43 13 44 14 xm m ym m zm m wm m ( ) ( ) ( ) ( )0 ++ ++ ++ + = Since w =1, we can simplify this as follows: 41 11 42 12 43 13 44 14 xm m ym m zm m m m ( ) ( ) ( )( )0 ++ ++ ++ + = This gives a standard plane equation: ax by cz d + + += 0 , where 41 11 42 12 43 13 44 14 am m bm m cm m dm m =+ =+ =+ =+ , ,, We have now shown that the left clipping plane of the viewing frustum can be directly extracted from the projection matrix. It is important to note, that the resulting plane equation is not normalized (i.e., the plane’s normal vector is not a unit vector), and that the normal vector is pointing to the inside halfspace. This means, that 0 < +++ ax by cz d , if the vertex v is in the inside halfspace of the left clipping plane. We can now repeat the above steps for the other inequalities to determine the plane equations for the remaining clipping planes. Clipping plane Inequality Coefficients for the plane equation left 4 1 1 1 4 4 ' ' ( ) 0 ( )( ) 0 ( ) w x row row row row row row − < −• < • < • +• < • + v v v v v 41 11 42 12 43 13 44 14 am m bm m cm m dm m = + = + = + = + right 1 4 4 4 1 1 ' ' 0 ( )( ( ) ) 0 x w row ro row row row ro w w < •< • < • −• < • − v v v v v 41 11 42 12 43 13 44 14 am m bm m cm m dm m = − = − = − = − bottom 4 2 2 2 4 4 ' ' ( ) 0 ( )( ) 0 ( ) w y row row row row row row − < −• < • < • +• < • + v v v v v 41 21 42 22 43 23 44 24 am m bm m cm m dm m = + = + = + = + top 2 4 4 4 2 2 ' ' 0 ( )( ( ) ) 0 y w row ro row row row ro w w < •< • < • −• < • − v v v v v 41 21 42 22 43 23 44 24 am m bm m cm m dm m = − = − = − = − near 4 3 3 3 4 4 ' ' ( ) 0 ( )( ) 0 ( ) w z row row row row row row − < −• < • < • +• < • + v v v v v 41 31 42 32 43 33 44 34 am m bm m cm m dm m = + = + = + = + far 3 4 4 4 3 3 ' ' 0 ( )( ( ) ) 0 z w row ro row row row ro w w < •< • < • −• < • − v v v v v 41 31 42 32 43 33 44 34 am m bm m cm m dm m = − = − = − = − 3.1 Supporting a Non-Identity Modelview Matrix Until now we have assumed that the modelview matrix is an identity matrices. However, the goal obviously is to make the algorithm work for an arbitrary view. This is, in fact, so easy that it’s almost unbelievable. If you think about it for a moment then you’ll immediately understand it, and that’s why we are not going to explain this in any great detail. All we are giving to you is this: 1. If the matrixM is equal to the projection matrix P (i.e., M P = ), then the algorithm gives the clipping planes in view space (i.e., camera space). 2. If the matrixM is equal to the combined projection and modelview matrices, then the algorithm gives the clipping planes in model space (i.e., VP M ⋅ = , where V is the modelview matrix, and P is the projection matrix). Appendix A – Plane Equations Given a fixed pointp and a non-zero normal vector n , a plane is defined as the set of all points x = (, ,) x y z , such that n xp •− = ( )0. In English that means, that the vector from p to any point x is orthogonal (perpendicular) to the plane’s normal vector n , if the point x is on the plane. Remember that the dot product of two vectors is zero, if the angle between those two vectors is 90 degrees. We are now going to bring the plane equation n xp • ( )0 − = into another form: n xp • ( )0 − = Is equal to: ( )( )0 nx np • −•= If n = ( ) abc , then we can rewrite the above as: ax by cz + +−•= ( )0 n p The term − • ( ) n p is constant and we now define that d :( ) = − • n p . Thus the plane equation becomes: ax by cz d + + += 0 It should be noted, that many people prefer to use the form ax by cz d + + = . In this case, the constant d is defined as d :( ) = • n p . For use in a computer language, we can store such a plane equation in a simple structure. struct Plane { float a, b, c, d; }; A.1 – Halfspaces A plane cuts three-dimensional space into two separate parts. These parts are called halfspaces. The halfspace the plane’s normals vector points into is called the positive halfspace, and the other halfspace is called the negative halfspace. A.2 – Normalizing the Plane Equation Normalizing a plane equation ax by cz d + + += 0 means to change the plane equation, such that the normal vector n = ( ) abc becomes a unit vector (i.e., n =1). All we need to do to normalize a vector is to divide each of its components by the magnitude of the vector. However, in order to normalize a plane equation it is not enough to simply divide the coefficients abc , , by the magnitude of the normal vector n . We also need to divide the constant d by the magnitude of the normal vector n : n xp • ( )0 − = First we normalize n by dividing its components by the magnitude of n : • ( )0 − = n x p n This is equal to: 0       • − •=    n n x p n n If n = ( ) abc , then we can rewrite the above as: 0 abc xyz   + + − •=     n p nnn n Now define d :( ) =− • n p : 0 ab cd xyz + + += nnnn Which is equal to: 222 222 222 222 0 ab cd xyz abc abc abc abc + + += ++ ++ ++ ++ Now simplify: 222 0 ax by cz d abc + + + = + + So in order to normalize a plane equation, we need to divide the coefficients abc , , and the constant d by the magnitude of the normal vector n = ( ) abc . It is useful to have a function that does this for you: void NormalizePlane(Plane & plane) { float mag; mag = sqrt(plane.a * plane.a + plane.b * plane.b + plane.c * plane.c); plane.a = plane.a / mag; plane.b = plane.b / mag; plane.c = plane.c / mag; plane.d = plane.d / mag; } A.3 – The Signed Distance From a Plane to a Point It is more than often necessary to compute the shortest signed distance from a given plane ax by cz d + + += 0 to a point p = (, ,) x y z . As it turns out, we can simply plug the coordinates of the point p into the plane equation, and the result is the signed distance from the plane top . This distance, however, is not necessarily a 'true' distance. Instead it is the signed distance in units of the magnitude of the plane’s normal vector n = ( ) abc . So in order to obtain a 'true' distance, you need to normalize the plane equation. float DistanceToPoint(const Plane & plane, const Point & pt) { return plane.a*pt.x + plane.b*pt.y + plane.c*pt.z + plane.d; } If the plane equation is not normalized, then we can still get some valuable information from the 'non-true' distance dist : 1. If dist < 0 , then the point p lies in the negative halfspace. 2. If dist = 0 , then the point p lies in the plane. 3. If dist > 0 , then the point p lies in the positive halfspace. This gives us another useful function that also works for non-normalized plane equations: enum Halfspace { NEGATIVE = -1 ON_PLANE = 0, POSITIVE = 1, }; Halfspace ClassifyPoint(const Plane & plane, const Point & pt) { float d; d = plane.a*pt.x + plane.b*pt.y + plane.c*pt.z + plane.d; if (d < 0) return NEGATIVE; if (d > 0) return POSITIVE; return ON_PLANE; } Appendix B – Implementation In the following we are going to give sample implementations of the plane extraction method for both, OpenGL and Direct3D. Note, that the normal vectors of the resulting planes are pointing to the inside of the viewing frustum. If you prefer planes where the normal vectors are pointing to the outside, then all you need to do is to negate the coefficients abc , , and the constant d of each plane equation. B.1 Plane Extraction for OpenGL struct Matrix4x4 { // The elements of the 4x4 matrix are stored in // column-major order (see "OpenGL Programming Guide", // 3rd edition, pp 106, glLoadMatrix). float _11, _21, _31, _41; float _12, _22, _32, _42; float _13, _23, _33, _43; float _14, _24, _34, _44; }; void ExtractPlanesGL( Plane * p_planes, const Matrix4x4 & comboMatrix, bool normalize) { // Left clipping plane p_planes[0].a = comboMatrix._41 + comboMatrix._11; p_planes[0].b = comboMatrix._42 + comboMatrix._12; p_planes[0].c = comboMatrix._43 + comboMatrix._13; p_planes[0].d = comboMatrix._44 + comboMatrix._14; // Right clipping plane p_planes[1].a = comboMatrix._41 - comboMatrix._11; p_planes[1].b = comboMatrix._42 - comboMatrix._12; p_planes[1].c = comboMatrix._43 - comboMatrix._13; p_planes[1].d = comboMatrix._44 - comboMatrix._14; // Top clipping plane p_planes[2].a = comboMatrix._41 - comboMatrix._21; p_planes[2].b = comboMatrix._42 - comboMatrix._22; p_planes[2].c = comboMatrix._43 - comboMatrix._23; p_planes[2].d = comboMatrix._44 - comboMatrix._24; // Bottom clipping plane p_planes[3].a = comboMatrix._41 + comboMatrix._21; p_planes[3].b = comboMatrix._42 + comboMatrix._22; p_planes[3].c = comboMatrix._43 + comboMatrix._23; p_planes[3].d = comboMatrix._44 + comboMatrix._24; // Near clipping plane p_planes[4].a = comboMatrix._41 + comboMatrix._31; p_planes[4].b = comboMatrix._42 + comboMatrix._32; p_planes[4].c = comboMatrix._43 + comboMatrix._33; p_planes[4].d = comboMatrix._44 + comboMatrix._34; // Far clipping plane p_planes[5].a = comboMatrix._41 - comboMatrix._31; p_planes[5].b = comboMatrix._42 - comboMatrix._32; p_planes[5].c = comboMatrix._43 - comboMatrix._33; p_planes[5].d = comboMatrix._44 - comboMatrix._34; // Normalize the plane equations, if requested if (normalize == true) { NormalizePlane(p_planes[0]); NormalizePlane(p_planes[1]); NormalizePlane(p_planes[2]); NormalizePlane(p_planes[3]); NormalizePlane(p_planes[4]); NormalizePlane(p_planes[5]); } } B.2 Plane Extraction for Direct3D struct Matrix4x4 { // The elements of the 4x4 matrix are stored in // row-major order. float _11, _12, _13, _14; float _21, _22, _23, _24; float _31, _32, _33, _34; float _41, _42, _43, _44; }; void ExtractPlanesD3D( Plane * p_planes, const Matrix4x4 & comboMatrix, bool normalize) { // Left clipping plane p_planes[0].a = comboMatrix._14 + comboMatrix._11; p_planes[0].b = comboMatrix._24 + comboMatrix._21; p_planes[0].c = comboMatrix._34 + comboMatrix._31; p_planes[0].d = comboMatrix._44 + comboMatrix._41; // Right clipping plane p_planes[1].a = comboMatrix._14 - comboMatrix._11; p_planes[1].b = comboMatrix._24 - comboMatrix._21; p_planes[1].c = comboMatrix._34 - comboMatrix._31; p_planes[1].d = comboMatrix._44 - comboMatrix._41; // Top clipping plane p_planes[2].a = comboMatrix._14 - comboMatrix._12; p_planes[2].b = comboMatrix._24 - comboMatrix._22; p_planes[2].c = comboMatrix._34 - comboMatrix._32; p_planes[2].d = comboMatrix._44 - comboMatrix._42; // Bottom clipping plane p_planes[3].a = comboMatrix._14 + comboMatrix._12; p_planes[3].b = comboMatrix._24 + comboMatrix._22; p_planes[3].c = comboMatrix._34 + comboMatrix._32; p_planes[3].d = comboMatrix._44 + comboMatrix._42; // Near clipping plane p_planes[4].a = comboMatrix._13; p_planes[4].b = comboMatrix._23; p_planes[4].c = comboMatrix._33; p_planes[4].d = comboMatrix._43; // Far clipping plane p_planes[5].a = comboMatrix._14 - comboMatrix._13; p_planes[5].b = comboMatrix._24 - comboMatrix._23; p_planes[5].c = comboMatrix._34 - comboMatrix._33; p_planes[5].d = comboMatrix._44 - comboMatrix._43; // Normalize the plane equations, if requested if (normalize == true) { NormalizePlane(p_planes[0]); NormalizePlane(p_planes[1]); NormalizePlane(p_planes[2]); NormalizePlane(p_planes[3]); NormalizePlane(p_planes[4]); NormalizePlane(p_planes[5]); } }