function PSF = psf_estimate_nnls(x, y, im, im_sharp, ksize, lambda)
	% psf_estimate_local(x, y, im, im_sharp, ksize, lambda)
	%
	% estimate PSF based on the method outlined in 
	%	PSF Estimation using Sharp Edge Prediction
	% Neel Joshi† Richard Szeliski∗ David J. Kriegman
	%
	% Input:
	%	x			x coordinate of image
	%	y			y coordinate of image
	%   im		    sampled blurry image
	%	im_sharp    sampled sharp image
	%	ksize		kernel size
	%	lambda		regularisation strength parameter
	% 
	% Output:
	%	PSF		estimated PSF ksize x ksize

	if nargin < 6
		lambda		= 2;	% regularisation strength parameter
	end
	

	%% extract sharp reference image within valid region
	assert(mod(ksize,2) == 1, 'kernel size must be odd');
	i0			= ceil(ksize/2);
	PSF			= zeros(ksize, ksize);
	PSF(i0,i0)  = 1;
	im_valid    = convn(im, PSF, 'valid');
	n_px		= numel(im_valid);

	%% build convolution matrix
	% represent convolution operation as matrix multiplication
	% C * PSF(:) = im(:)
	n_psf		= ksize^2;
	C			= zeros(n_px, n_psf);
	Kx			= zeros(n_psf, n_psf);
	Ky			= zeros(n_psf, n_psf);
	im_sharp_pad= padarray(im_sharp, floor([ksize ksize]/2), 'replicate', 'both');

	for n = 1 : n_psf
		% represent convolution operation in matrix
		PSF		= zeros(ksize, ksize);
		PSF(n)  = 1;
		im_conv = convn(im_sharp, PSF, 'valid');
		C(:,n)  = im_conv(:);
		% represent gradient penalty
		[ky,kx] = gradient(PSF);
		Kx(:,n) = kx(:);
		Ky(:,n) = ky(:);
	end

	%% check
	% check convolution matrix constructed correctly
% 	PSF  = randn(ksize);
% 	im_conv_C   = C * PSF(:);
% 	im_conv_ref = convn(im_sharp, PSF, 'same');
% 	disp(norm(im_conv_ref(:)-im_conv_C(:)))

	%% estimate noise
	sigma		= estimate_noise(im);

	%% build NNLS system
	% express least squares minimisation of
	%	|im - convn(im_sharp, PSF)|^2 + sigma^2*lambda*gamma*|grad PSF|^2
	% i.e. equation 5 of [1]
	% in the form
	%	|C * k - d|^2

	% regularisation factor
	gamma		= ksize^2 * n_px;	% normalisation for kernel size and number of observations
	Gamma		= sqrt(lambda*gamma*sigma^2);
	C_full		= [C; Gamma*Kx; Gamma*Ky];
	d_full		= [im_valid(:); zeros(n_psf*2,1)];

	%% solve regularised NNLS system
	opts		= optimset('Display', 'final');
	PSF			= lsqnonneg(C_full, d_full, opts);
	%PSF			= lsqlin(C_full, d_full);
	PSF			= reshape(PSF, ksize, ksize);
end

function Sigma=estimate_noise(I)
	% by Tolga Birdal
	% This is an extremely simple m-file which implements the method described
	% in : J. Immerk�r, �Fast Noise Variance Estimation�, Computer Vision and
	% Image Understanding, Vol. 64, No. 2, pp. 300-302, Sep. 1996
	% 
	% The function inputs a grayscale image I and returns Sigma, the noise
	% estimate. Here is a sample use: 
	% 
	% I = rgb2gray(imread('sample.jpg'));
	% Sigma=estimate_noise(I);
	%         
	% The advantage of this method is that it includes a Laplacian operation 
	% which is almost insensitive to image structure but only depends on the 
	% noise in the image. 


	[H, W]=size(I);
	I=double(I);
	% compute sum of absolute values of Laplacian
	M=[1 -2 1; -2 4 -2; 1 -2 1];
	Sigma=sum(sum(abs(conv2(I, M))));
	% scale sigma with proposed coefficients
	Sigma=Sigma*sqrt(0.5*pi)./(6*(W-2)*(H-2));
end