function [ R ] = deconv_sondhi( B, PSF )
% Reference: http://ieeexplore.ieee.org/xpls/abs_all.jsp?arnumber=1450713

% Hack to incooperate PSFs. Only linear uniform PSFs are acceptable and
% therefore produce valid results

a = length(PSF)-1;

% Copy blurred image to fill target with some values
R = B;

% Define commonly used vars
L = size(B,2);           % Length
N = fix(L/a);            % Def. of N

% Loop over all color channels
for c = 1:size(B,3)
    for y = 1:size(B,1)
        % Compute line average
        g = B(y,:,c);         % One image line, g
        avg = mean(g);      % Average of image line, <g>
        dg = gradient(g);   % derivate of first kind of g, g'

        for x = 1:size(B,2)-1

            % Compute first sum
            s1 = 0;
            for j = 0:(N-1)
                for k = 0:j
                    s1 = s1 + dg(x - fix(x/a)*a + j*a - k*a + 1);
                end
            end

            % Compute second sum
            s2 = 0;
            for k = 0:fix(x/a)
                s2 = s2 + dg(x - k*a + 1);
            end

            % Put all parts together
            px = avg - a/N*s1 + a*s2;
            R(y,x,c) = px;
        end
    end
end

% Perform shift of image
for i = 1:size(B,3)
    R(:,:,i) = (circshift(R(:,:,i)',fix(length(PSF)/2)))';
end
