Image convolution and deconvolution can be performed quickly in the frequency domain. In the frequency domain, pointwise multiplication is convolution, while pointwise division is deconvolution. The pseudo-code of the procedure is:

1. Convert your input image into the frequency domain, via fast Fourier transform.
2. Pad your convolution kernel with zeros, until it reaches the same size as your input image.
3. Convert your convolution kernel into the frequency domain, via fast Fourier transform.
4. Pointwise multiply the frequency domain input image and the frequency domain convolution kernel together.
5. Convert the product from pointwise multiplication back into the spatial domain via inverse fast Fourier transform.

If you want to perform deconvolution, rather than perform pointwise multiplication in step 4, you do pointwise division. You have to be aware divide-by-zero problem, when you are performing deconvolution. You might want to replace all the zeros with eps().

Below are some Matlab examples. Note that Matlab takes care of the zero-padding inside the fft2() function call.

function [ out ] = get_laplacian( in )
%GET_LAPLACIAN Obtain the Laplacian of an image
% This is an implementation of the algorithm described in page 32-35 of 
% Lightness and Brightness Computation by Retinex-like Algorithms by Eran 
% Borenstein
 
if size(in, 3) == 3
    for i = 1:3
        tmp = in(:,:,i);
        out(:,:,i) = get_laplacian_real(tmp);
    end
else
    out = get_laplacian_real(in);
end
 
end
 
function [ out ] = get_laplacian_real(in)
%% Fourier transform input image
a = fft2(in);
b = fft2(fspecial('laplacian'), size(in, 1), size(in, 2));
c = a.*b;
out = real(ifft2(c));
end

function [ out ] = inv_laplacian(in)
%INV_LAPLACIAN Convert a Laplacian image back to its original
 
if size(in, 3) == 3
    for i = 1:3
        tmp = in(:,:,i);
        out(:,:,i) = inv_laplacian_real(tmp);
    end
else
    out = inv_laplacian_real(in);
end
 
end
 
function [out] = inv_laplacian_real(in)
%% Fourier transform input image
a = fft2(in);
b = fft2(fspecial('laplacian'), size(in, 1), size(in, 2));
b(b == 0) = eps;
c = a./b;
out = real(ifft2(c));
out = out + abs(min(out(:)));
end