function [M,U,xc] = frustum(cmodel)
	% [M,U,xc] = frustum(cmodel)
	% 
	% Calculate the viewing frustum for the pinhole camera model specified
	% by cmodel
	%
	% M is a 4x4 matrix whose rows provide the equation of plane for each
	% plane of the viewing frustum. That is to say, for each plane i,
	% points satisfying:
	%     M(i,1)*x + M(i,2)*y + M(i,3)*z + M(i,4) = 0
	% lie on that plane
	%
	% The plane normals are oriented so that (0,0,0) always corresponds to a
	% positive distance
	%
	% U is a 3x4 matrix whose columns describe the unit vectors of lines
	% describing the viewing frustum. Each edge i of the frustum is given by 
	% the equation of a line:
	%
	%     x(s) = xc + u(:,i)*s
	
	%% extract pinhole model parameters
	i_axis			= cmodel.i_axis;
	j_axis			= cmodel.j_axis;
	ic					= cmodel.calib.ic;
	kr					= cmodel.calib.kr;
	G					= cmodel.calib.G;
	R					= cmodel.calib.R;
	GR					= G*R;
	xc					= cmodel.calib.xc;
	
	%% calculate corners
	ij_corners		= [i_axis([1 1 end end]);
							j_axis([1 end end 1])];
	% account for distortion
	ij_corners		= inverse_radial_distortion(ij_corners, ic, kr);
	
	%% find unit vectors
	ij_tilde			= [ij_corners; ones(1,4)];
	U					= GR \ ij_tilde;
	U_norm			= sqrt(sum(U.^2,1));
	U					= U ./ ([1;1;1]*U_norm);
	% orientation is such that u . xc < 0
	flip				= -sign(U.' * xc)';
    flip(flip==0)       = 1; % no change if exactly zero
	U					= U .* ([1;1;1]*flip);
	
	%% find equations of plane
	M					= zeros(4,4);
	for c = 1 : 4
		m				= cross(U(:,c), U(:,mod(c,4)+1));
		m				= m / norm(m);
		d				= -dot(m, xc);
		% plane normals are oriented to ensure (0,0,0) is always on +ve side
		if d<0
			d			= -d;
			m			= -m;
		end
		M(c,:)		= [m', d];
	end
end