2D Convolution in MATLAB Causes Artifacts (Boundary Issues)

Question

I`m trying to do a 2D fast convolution in Matlab of large matrices. If I use FFT version based on convolution theorem ( https://en.wikipedia.org/wiki/Convolution_theorem ), there are some artefacts in the image. Only imfilter produces correct results, but it is working ~100x slower, than the FFT version (conv2fft).

What can be done to avoid these artefacts? I do know that in 1D we need to zero-pad, but here both the image and the kernel have the same size.

% generate image
len = 2^10;
CICcut = zeros (len);
CICcut = imnoise (CICcut, 'salt & pepper', 0.0001);
CICcut = CICcut.*(rand(len)).^2;
gauss = fspecial('gaussian', round(sqrt(len)), sqrt(sqrt(len)));
CICcut = imfilter (CICcut, gauss, 'replicate', 'conv');

% generate kernel
g = zeros(len);
lenMone = len-1;
for i = 1:len
    for j = 1:len
        g(i, j) = ((i-1)/lenMone - 0.5)^2 + ((j-1)/lenMone - 0.5)^2;
    end
end
g = -log(sqrt(g));

% convolution
tic
filtered    = imfilter (g, CICcut, 'replicate', 'conv');
toc
tic
filteredFFT = conv2fft (g, CICcut, 'same');
toc
tic
filteredN   = convn (g, CICcut, 'same');
toc

% display
figure('units', 'normalized', 'outerposition', [0 0.25 1 0.5])
subplot 151, imshow (CICcut, []); title ('Mass density')
subplot 152, imshow (g, []); title ('Green`s function')
subplot 153, imshow (filtered, []); title ({'Gravitational potential' 'imfilter'})
subplot 154, imshow (filteredFFT, []); title ({'Gravitational potential' 'conv2fft'})
subplot 155, imshow (filteredN, []); title ({'Gravitational potential' 'convn'})

Best regards, Alex

Answer 1

If I understood right this is a duplicate question but I cannot find the link let me re-type the answer. Given that you want to implement a 2D discrete convolution of two images of sizes $N_1 \times N_2$ and $M_1 \times M_2$ by the method of DFT multiplication, then in order to avoid time (space) domain aliasing due to circular convolution implied by the inverse DFT stage, the sizes of the forward DFTs should be properly selected (in other words original time/space domain signals should be padded with zeros to increase thier sizes to overcome the aliasing).

Particularly the 2D forward and inverse DFT sizes should be selected as: $$ L_1 \geq N_1 + M_1 -1$$ and $$ L_2 \geq N_2 + M_2 -1$$ in order to avoid circular artifacts and get the exact convolution.

Select the central portion of the resulting convolution to get the final image the same size of the original image.