Using @Shai's relevant suggestion, I can now give the detailed answer.
Options 2 and 3 that I suggested are in fact the same and are apparently the right way to go. It is also what was done in the E-step of the suggested paper linked by @Shai. Posing the over constrained problem is in fact quite simple.
To correctly pose these equations we use the fact that the dot product of every block the size of the kernel centered around some pixel in G
with the 180 degrees rotated version of h
should equal the corresponding pixel in B
. This is directly derived from the fact that B
and G
are related by convolution and so blocks in G
relate to pixels in B
by cross-correlation (and hence the 180-degree rotation).
The MATLAB code now becomes:
%inputs: B,G - gray level blurred and sharp images respectively (double)
% szKer - 2 element vector specifying the size of the required kernel
%outputs: mKer - the recovered kernel,
% imBsynth - the sharp image convolved with the recovered kernel
%
%example usage: mKer = calcKer(B, G, [11 11]);
function [mKer, imBsynth] = calcKer(B, G, szKer)
%get the "valid" pixels from B (i.e. those that do not depend
%on zero-padding or a circular assumption
imBvalid = B(ceil(szKer(1)/2):end-floor(szKer(1)/2), ...
ceil(szKer(2)/2):end-floor(szKer(2)/2));
%get a matrix where each row corresponds to a block from G
%the size of the kernel
mGconv = im2col(G, szKer, 'sliding')';
%solve the over-constrained system using MATLAB's backslash
%to get a vector version of the cross-correlation kernel
vXcorrKer = mGconv \ imBvalid(:);
%reshape and rotate 180 degrees to get the convolution kernel
mKer = rot90(reshape(vXcorrKer, szKer), 2);
if (nargout > 1)
%if there is indeed a convolution relationship between B and G
%the following will result in an image similar to B
imBsynth = conv2(G, mKer, 'valid');
end
end
I also found that some constraints over the solution may be necessary for practical scenarios. Examples are enforcing the kernel to be positive, smooth, or symmetric. The exact method of incorporating these is out of the scope of this question but will usually come in the shape of a linear constraint or a regularization element when solving for vXcorrKer
.