Linear convolution of two images in Matlab using fft2

Question

I would like to take two images and collapse them together in Matlab using 2D-FFT, without resorting to the conv2 function. However, I am not sure how matrices should be appropriately supplemented and prepared for convolution.

The mathematical operation is as follows:

A * B = C

In the above, * is the convolution operator ( Wikipedia link ).

The following Matlab program shows the difference between padding, not padding. I suspect that not filling the matrices leads to circular convolution, but I would like to perform linear convolution without smoothing.

If I overlap two matrices, then how do I crop the convolution output so that C is the same size as A and B

 A = rgb2gray(im2double(imread('1.png'))); % input A B = rgb2gray(im2double(imread('2.png'))); % kernel B figure; imagesc(A); colormap gray; title ('A') figure; imagesc(B); colormap gray; title ('B') [m,n] = size(A); mm = 2*m - 1; nn = 2*n - 1; C = (ifft2(fft2(A,mm,nn).* fft2(B,mm,nn))); figure; imagesc(C); colormap gray; title ('C with padding') C0 = (ifft2(fft2(A).* fft2(B))); figure; imagesc(C0); colormap gray; title ('C without padding') 

Here is the result of the program:

Answer 1

Without filling, the result will be equivalent to a circular convolution , as you specify. For linear convolution, when convolving 2 images (2D signals) A * B, the total output will have the size Ma+Mb-1 x Na+Nb-1 , where Ma x Na, Mb x Nb sizes of images A and B, respectively.

After adding to the expected size, multiplying and converting backwards, through ifft2 , you can save the central part of the resulting image (usually corresponding to the largest of A and B).

 A = double(imread('cameraman.tif'))./255; % image B = fspecial('gaussian', [15 15], 2); % some 2D filter function [m,n] = size(A); [mb,nb] = size(B); % output size mm = m + mb - 1; nn = n + nb - 1; % pad, multiply and transform back C = ifft2(fft2(A,mm,nn).* fft2(B,mm,nn)); % padding constants (for output of size == size(A)) padC_m = ceil((mb-1)./2); padC_n = ceil((nb-1)./2); % frequency-domain convolution result D = C(padC_m+1:m+padC_m, padC_n+1:n+padC_n); figure; imshow(D,[]); 

Now compare the above with performing spatial domain convolution using conv2D

  % space-domain convolution result F = conv2(A,B,'same'); figure; imshow(F,[]); 

The results are visually the same, and the total error between them (due to rounding) is of the order of e-10.