Image Convolution Using DFT (FFT)

Question

I just can't get my head around Fourier transform and convolution in 2D. I am trying to implement image convolution using fast Fourier transform (in julia). So the first thing I need to do is to pad the image with zeros. Let's assume that the image is $N\times P$ and the kernel $s \times q$, so I would pad the image so that it is of size $(N+s-1)\times (P+q-1)$.

(N, P) = size(image)
(s, q) = size(kernel)
image_padded = zeros(N+s-1, P+q-1)
image_padded[1:N, 1:P] = image

Then I need to make the kernel the same size as the image, so I pad it with zeros, too, but this time, I have to put the kernel in the center of the huge block of zeros.

kernel_padded = zeros(N+s-1, P+s-1)
Npul = Int((N+s-2)/2) #both image and kernel have odd sizes
Ppul = Int((P+s-2)/2)
spul = Int((s-1)/2)
qpul = Int((q-1)/2)
kernel_padded[Npul-spul+1:Npul+spul+1, Ppul-spul+1:Ppul+spul+1] = k

Then I take fft of both padded image and padded kernel, take their Hadamard product and use inverse fft to get the result back to image domain.

real.(ifft(fft(image_padded).*fft(kernel_padded)))

But the result does not look as I would expect, there is like a cross in the middle of my resulting image. Actually, it looks like there is one quater of the image in each corner but rotated so that the border effect creates the dark cross in the middle. The only way I can get rid of it is to pad the image and kernel to $(N+N+1)\times(P+P+1)$.

So do I understand the process wrong or is something wrong with my code? Thanks for any advice.

Answer 1

First, You use fft() which is 1D DFT.
For 2D convolution you need 2D DFT which is implemented in fft2().

A reference implementation is given in Replicate MATLAB's conv2() in Frequency Domain.

You seem to be after the full case which means you need to pad signals as you did.