Frustum Clipping for PovRay

In normal rendering applications the potential exists for any piece of the geometry making up a scene to be visible or to have some effect on the visible part of the scene. For example, objects that can't be directly seen may be visible in reflective surfaces or they may cast shadows into the visible area.

It is however quite common in scientific visualisation projects that there is simply too much data for PovRay to render directly. One needs to get up to "tricks", for example, creating geometry at variable resolutions depending on the distance from the camera. Another trick that will discussed here is to prune away any geometry not directly visible. The particular rendering here involved a topology model of Mars containing over 130 million polygons. The polygonal approximation was adjusted depending on the distance of the camera to the surface to give about one polygon per pixel. In order to further reduce the polygons given to PovRay the polygons outside the view frustum were removed.

In order to compute which polygons are outside the frustum one needs to be able to define the 4 planes making up the view frustum against which each vertex of the polygons will be tested. If one defines the horizontal aperture as thetah then the vertical aperture thetav is given by the following.

The points p0, p1, p2, p3 making up the corners of frustum can be computed using the C source below.

   p0 = vp;
   p0.x += vd.x - right.x*tan(thetah/2) - vu.x*tan(thetav/2);
   p0.y += vd.y - right.y*tan(thetah/2) - vu.y*tan(thetav/2);
   p0.z += vd.z - right.z*tan(thetah/2) - vu.z*tan(thetav/2);
   p1 = vp;
   p1.x += vd.x + right.x*tan(thetah/2) - vu.x*tan(thetav/2);
   p1.y += vd.y + right.y*tan(thetah/2) - vu.y*tan(thetav/2);
   p1.z += vd.z + right.z*tan(thetah/2) - vu.z*tan(thetav/2);
   p2 = vp;
   p2.x += vd.x + right.x*tan(thetah/2) + vu.x*tan(thetav/2);
   p2.y += vd.y + right.y*tan(thetah/2) + vu.y*tan(thetav/2);
   p2.z += vd.z + right.z*tan(thetah/2) + vu.z*tan(thetav/2);
   p3 = vp;
   p3.x += vd.x - right.x*tan(thetah/2) + vu.x*tan(thetav/2);
   p3.y += vd.y - right.y*tan(thetah/2) + vu.y*tan(thetav/2);
   p3.z += vd.z - right.z*tan(thetah/2) + vu.z*tan(thetav/2);

Where vp is the view position vector, vd is the unit view direction vector, vu is the unit up vector, and right is the unit vector to the right (the cross product between vd and vu). The 4 frustum planes are (vp,p0,p1), (vp,p1,p2), (vp,p2,p3), (vp,p3,p0) from which the normal of each plane can be computed.

A simple C function that determines which side of a plane a vertex lies might be as follows.

/*
   Determine which side of a plane the point p lies.
   The plane is defined by a normal n and point on the plane vp.
   Return 1 if the point is on the same side as the normal, 
   otherwise -1.
*/
int WhichSide(XYZ p,XYZ n,XYZ vp)
{
   double D,s;

   D = -(n.x*vp.x + n.y*vp.y + n.z*vp.z);
   s = n.x*p.x + n.y*p.y + n.z*p.z + D;

   return(s<=0 ? -1 : 1);
}

The example below shows a portion of the landscape rendered on the left. Moving the camera back a bit shows that the geometry (of the whole planet) has been clipped to remove any polygons not within the view frustum. Note that as well as frustum clipping, back facing polygons have also been removed.