function PSF = psf_fit_fmincon(im, im_sharp, ksize, lambda)
	% psf_fit_fmincon(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:
	%   im		    sampled blurry image
	%	im_sharp    sampled sharp image (may be upsampled)
	%	ksize		kernel size
	%	lambda		regularisation strength parameter
	% 
	% Output:
	%	PSF		estimated PSF ksize x ksize

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

	dtype		= 'gpuArray';
	
	%% sanity check
	% size of particle images/sharp images
	sz_sharp	= size(im_sharp);
	sz_im		= size(im);
	n_img		= prod(sz_im(3:end));
	n_px		= prod(sz_im(1:2));
	n_px_sharp  = prod(sz_sharp(1:2));

	assert(prod(sz_sharp(3:end)) == n_img, 'number of particle images must match number of sharp particle images');
	% antialiasing factor
	aa_factor   = sz_sharp(1) / sz_im(1);
	assert(mod(aa_factor,1) == 0 && aa_factor >= 1, ...
		'size of sharp image must be an integer multiple of the size of the particle image');
	assert(sz_sharp(2) == aa_factor*sz_im(2), ...
		'antialiasing factor must be uniform in both dimensions of the sharp image');
	% size of the kernel
	assert(mod(ksize,2) == 1, 'kernel size must be odd');
	n_psf		= ksize^2;

	% get size of padded sharp image
	im_sharp    = double(im_sharp);
	im_sharp_pad= padarray(im_sharp(:,:,1), floor([ksize ksize]/2), 'replicate', 'both');
	sz_pad		= size(im_sharp_pad);

	%% build anti-aliasing matrix
	aa_mtx		= psf_aa_matrix(sz_im(1:2), aa_factor);

	%% build gradient penalty matrix
	Kx			= zeros(n_psf, n_psf);
	Ky			= zeros(n_psf, n_psf);
	C_idx		= zeros(n_px_sharp, n_psf);
	im_pad_idx  = reshape(1 : prod(sz_pad), sz_pad);
	
	for n = 1 : n_psf
		% represent convolution operation in matrix
		PSF		= zeros(ksize, ksize, dtype);
		PSF(n)  = 1;

		% create mapping for convolution operation
		im_conv_idx= conv2(im_pad_idx, PSF, 'valid');
		C_idx(:,n) = im_conv_idx(:);
		
		% represent gradient penalty
		[ky,kx] = gradient(PSF);
		Kx(:,n) = kx(:);
		Ky(:,n) = ky(:);
	end

	%% build convolution matrix
	% represent convolution operation as matrix multiplication
	% C * PSF(:) = im(:)
	CtC			= zeros(n_psf, n_psf, dtype);
	Cty			= zeros(n_psf, 1, dtype);
	if strcmp(dtype, 'gpuArray')
		 C_idx	= gpuArray(C_idx);
		aa_mtx  = gpuArray(aa_mtx);
	end
			
	for page = 1 : n_img
		%% pad sharp image
		im_sharp_pad= padarray(im_sharp(:,:,page), floor([ksize ksize]/2), 'both', 'replicate');
		if strcmp(dtype, 'gpuArray')
			im_sharp_pad= gpuArray(im_sharp_pad);
		end
		
		%% build convolution matrix
		conv_mtx	= im_sharp_pad(C_idx);
		% apply anti-aliasing operation
		C			= aa_mtx * conv_mtx;

		%% accumulate Jacobian
		CtC			= CtC + C.'*C;
		y			= im(:,:,page);
		Cty			= Cty + C.' * y(:);
	end
	CtC			= gather(CtC);
	Cty			= gather(Cty);
	
	%% estimate noise
	sigma		= estimate_noise(im);
	% mean over all images
	sigma		= sqrt(mean(sigma(:).^2));

	%% 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 * x - y|^2 = x' * C'C * x - 2*x'C'y + y'y
	% we evaluate the matrices C'C and C'y in the loop above
	% then add the regularisation factor below

	% regularisation factor
	gamma		= ksize^2 * n_px;	% normalisation for kernel size and number of observations
	Gamma		= sqrt(lambda*gamma*sigma^2);
	C_reg		= [Gamma*Kx; Gamma*Ky];
	CtC_reg		= C_reg.' * C_reg;
	% add 
	CtC			= CtC + CtC_reg;
	
	%% solve regularised LLS system
	% the linear least squares minimiser is
	% (C'C)\C'y
	PSF_lls		= CtC \ Cty;
	PSF_fg		= max(0, PSF_lls);
	
	%% solve constrained optimisation problem
	% we use fmincon to solve the optimisation where PSF>0
	% constraint is formulated as A*x <= b
	
	opts		= optimoptions('fmincon', ...
		'Display', 'final', ...
		'SpecifyObjectiveGradient', true, ...
		'HessianFcn', @(x,lambda) 2*CtC);
	PSF			= fmincon(@(x) quadratic_objective(x,CtC,Cty), ...
		PSF_fg, -eye(ksize^2), zeros(ksize^2,1), [], [], [], [], [], opts);
	[f,grad]=quadratic_objective(PSF_lls,CtC,Cty);

	%% reshape
	PSF			= reshape(PSF, ksize, ksize);
end

function [f, grad] = quadratic_objective(x, CtC, Cty)
	% we evaluate the objective function of the linear least squares
	% problem
	%	f = |C * x - y|^2 = x' * C'C * x - 2*x'C'y + y'y
	% we also evaluate the gradient df/dx and the hessian d2f/dx^2

	f	= x.' * CtC * x - 2*x.'*Cty;
	grad= 2*CtC * x - 2*Cty;
	%hess= 2*CtC;
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     = size(I,1);
	W     = size(I,2);
	% compute sum of absolute values of Laplacian
	M=[1 -2 1; -2 4 -2; 1 -2 1];
	sigma = sum(sum(abs(convn(I, M))));
	% scale sigma with proposed coefficients
	sigma = sigma*sqrt(0.5*pi)./(6*(W-2)*(H-2));
	sigma = squeeze(sigma);
end