function mask = calcmask(camset, c0)
	% mask = calcmask(camset, c)
	%
	% for each pixel in camera c of camset, determine the number
	% of other cameras its epipolar ray is visible in
	
	%% get things
	n_cameras= length(camset);
	cam0		= camset(c0);
	Ni			= length(cam0.i_axis);
	Nj			= length(cam0.j_axis);
	Npx		= Ni*Nj;
	
	%% dewarp
	[im,jm]	= ndgrid(cam0.i_axis, cam0.j_axis);
	ui_pos	= inverse_radial_distortion([im(:),jm(:)]', cam0.calib.ic, cam0.calib.kr);
	
	%% get equation of ray
	% we write the equation of ray as
	% x = x0 + s*e_s
	ytilde	= [ui_pos; ones(1,Npx)];
	GR			= cam0.calib.G * cam0.calib.R;
	e_ray		= GR \ ytilde;
	x0_ray	= cam0.calib.xc;
	
	%% determine mask
	mask				= zeros(Ni, Nj);
	for c1 = 1 : n_cameras
		%% skip if same camera
		if c0 == c1
			continue;
		end
		
		%% find linear inequality for bounds test
		% the test for whether a point y is within the bounds of image space
		% of this camera is described as an inequality Ay <= b. If satisfied,
		% the point is within bounds, if not, it is outside bounds. we
		% construct this inequality based on borders of image in undistorted
		% image space
		cam1			= camset(c1);
		i_axis		= cam1.i_axis(:);
		j_axis		= cam1.j_axis(:);
		Ni1				= length(i_axis);
		Nj1			= length(j_axis);
		ir				= [1 : floor(Ni1/4) : Ni1, Ni1];
		jr				= [1 : floor(Nj1/4) : Nj1, Nj1];
		i_corners	= [i_axis(ir([1 1 1 1 1]));
							i_axis(ir             );
							i_axis(ir([5 5 5 5 5]));
							i_axis(ir(end:-1:1   ))];
		j_corners	= [j_axis(jr             );
							j_axis(jr([5 5 5 5 5]));
							j_axis(jr(end:-1:1   ));
							j_axis(jr([1 1 1 1 1]))];
		di_corners	= [i_corners'; j_corners'];
		ui_corners	= inverse_radial_distortion(di_corners, cam1.calib.ic, cam1.calib.kr);
		[A,b]			= bounding_polygon(ui_corners);

		%% get equation of ray in image coordinates
		% we want an equation of the form
		% y = y0 + s * n
		% that describes the equation of each ray from camera 0 in the
		% (undistorted) image space of camera 1
		
		GR1					= cam1.calib.G * cam1.calib.R;
		z0						= GR1(1:2,:)*(x0_ray - cam1.calib.xc);
		zprime				= GR1(1:2,:)*e_ray;
		lambda0				= GR1(  3,:)*(x0_ray - cam1.calib.xc);
		lprime				= GR1(  3,:)*e_ray;
		
		y0						= bsxfun(@rdivide, zprime, lprime);
		e_n					= bsxfun(@minus, z0, lambda0*y0);
		e_n					= bsxfun(@rdivide, e_n, sqrt(sum(e_n.^2,1))); % normalise
		
		%% determine existence of solutions for inequality
		% subsituting in y = y0 + s*e_n
		% into our inequality Ay <= b
		% yields the inequality
		%		L*s <= R
		% where L = A*e_n and R = b - A*y0
		% if this inequality has solutions for any s, then there is at least
		% one point on the ray which is visible in camera 1
		
		L						= A*e_n;
		R						= bsxfun(@minus, b, A*y0);
		gamma_lim			= R ./ L;
		is_ub					= L >= 0;
		is_lb					= L <  0;
		
		% find smallest upper bound
		gamma_ub				= gamma_lim;
		gamma_ub(is_lb)	= inf;	% ignore all those lower bounds
		gamma_ub				= min(gamma_ub, [], 1);
		
		% find largest lower bound
		gamma_lb				= gamma_lim;
		gamma_lb(is_ub)	= -inf;	% ignore all the upper bounds
		gamma_lb				= max(gamma_lb, [], 1);
		
		b_visible			= gamma_lb <= gamma_ub;
		
		%% add to sum
		b_visible			= reshape(b_visible, [Ni Nj]);
		mask					= mask + b_visible;
	end
end
	