function [pimg_sharp, pimg_fit, BW] = pimg_fit_sharp(im_subpx, im_radius, pimg, PSF0, irel_aa, jrel_aa, aa_factor)
	% create an initial guess of a sharp particle image
	% based on a given sub-pixel position, radius and convolution kernel
	% 
	% Input:
	%	im_subpx	subpixel position 2 x n_pimg
	%	im_radius	particle radius   1 x n_pimg

	[a,n_pimg]  = size(im_subpx);
	assert(a==2, 'im_subpx must be 2 x  n_pimg');
	assert(isvector(im_radius) && numel(im_radius) == n_pimg, 'im_radius must be 1 x n_pimg');
	
	[n_support,a,b] = size(pimg);
	assert(a==n_support && b==n_pimg, 'pimg must be n_support x n_support x n_pimg');

	ksize		= size(PSF0);
	n_support	= size(pimg,1);

	%% create binarised particle image
	im_subpx	= reshape(im_subpx, [2 1 n_pimg]);
	im_radius	= reshape(im_radius, [1 1 n_pimg]);

	rmat		= sqrt((irel_aa(:) - im_subpx(1,1,:)).^2 + (jrel_aa(:).' - im_subpx(2,1,:)).^2);
	pimg_bw		= pimg_step(rmat, im_radius(1,1,:), 1);
	pimg_bw_pad = padarray(pimg_bw, floor(ksize/2), 'replicate', 'both');
	
	%% convolve with initial guess of PSF
	pimg_fg_white	= psf_downsample(convn( pimg_bw_pad, PSF0, 'valid'), aa_factor);
	pimg_fg_black	= 1 - pimg_fg_white;
	% fit a linear model 
	%	pimg = W * pimg_fg_white + B * pimg_fg_black
	% for unknown parameters B and W (white and black levels)
	% 
	C			= cat(2, reshape(pimg_fg_black, [n_support^2 1 n_pimg]), ...
				 	 	reshape(pimg_fg_white, [n_support^2 1 n_pimg]));
	y			= reshape(pimg, [n_support^2 1 n_pimg]);
	BW			= pagemldivide(C,y);
	% optimal fitted particle image
	pimg_fit	= pimg_fg_black.*BW(1,1,:) + pimg_fg_white.*BW(2,1,:);
	% optimal sharp image
	pimg_sharp  = pimg_bw.*BW(2,1,:) + (1-pimg_bw).*BW(1,1,:);
end