function y = inverse_radial_distortion(yprime, p, k)
	% y = inverse_radial_distortion(yprime, p, k)
	%
	% reverse radial distortion, to map from
	% (distorted)   real world image coordinates y'
	% (undistorted) pinhole image coordinates y
	% 
	% these are related according to
	% y' - p = (y-p)(1 + k(1)*|y-p|^2 + k(2)*|y-p|^4)
	%
	% yprime is 2 x N matrix
	% p is 1x2 vector, representing the principal point
	% k is 1x2 vector, representing distortion coefficients
	
	% let z' = y' - p
	%     z  = y - p
	%     r' = |z'|
	%     r  = |z|
	% then solve
	% r' = r*(1 + k(1)r^2 + k(2)*r^4)
	% for r, for the root nearest r'
	zprime	= bsxfun(@minus, yprime, p(:));
	rprime	= sqrt(sum(zprime.^2, 1));
	
	% use newton-raphson method, just a few iterations are sufficient
	r			= rprime;
	f			= @(r) k(2)*r.^5 + k(1)*r.^3 + r - rprime;
	fprime	= @(r) 5*k(2)*r.^4 + 3*k(1)*r.^2 + 1;
	for it = 1 : 3
		r		= r - f(r) ./ fprime(r);
	end
	
	% now calculate y
	z			= bsxfun(@rdivide, zprime, 1 + k(1)*r.^2 + k(2)*r.^4);
	y			= bsxfun(@plus, z, p(:));
end