I am trying to convolve an image with the code above using Convolution Theory and numpy's Fourier transform. However, my output seems to be slightly different than the result from scipy. I am not sure why my values are slightly off.

The inputs image is a numpy array, and the kernel input is a normalised 3x3 kernel.

Any help would be greatly appreciated.

def fft_convolution(image, kernel):
  image_padding = len(kernel)//2

  #padding image
  image_pad = np.pad(image_pad,(((image_padding+1)//2,image_padding//2),((image_padding+1)//2, image_padding//2)), mode='edge') 
  pad_x = image_pad.shape[0] - kernel.shape[0]
  pad_y = image_pad.shape[1] - kernel.shape[1]

  #pad kernel so it is same size as image
  kernel = np.pad(kernel, (((pad_x+1)//2,pad_x//2),((pad_y+1)//2,pad_y//2)), 'constant')

  #move the kernel so center is in the top left corner
  kernel = np.fft.ifftshift(kernel)

  #convert both to fft
  img_fft = np.fft.fft2(image_pad)
  kernel_fft = np.fft.fft2(kernel)

  #multiply 2 fourier matrices
  img = img_fft * kernel_fft

  #inverse fft
  img_inverse = np.fft.ifft2(img)

  #take the real numbers
  output = np.real(img_inverse)
  #slice array to get rid of padding and original size of image
  return output[image_padding:len(image_pad),image_padding:len(image_pad)]

I'm also finding increasing the kernel size (4x4,5x5,etc) results in a bigger change in my output compared to the image.

Edit1: So made some kind of a breakthrough. It was the slicing of my array that seemed vital to getting the right answer.

For 3x3 Kernel, output[1:len(image_pad)-1, 1:len(image_pad)-1] works perfect

4x4: output[2:len(image_pad)-1, 2:len(image_pad)-1]

However, I am unable to work out exactly what's going on and I can't find the right slicing for 5x5.

6x6: output[3:len(image_pad)-2, 3:len(image_pad)-2]

7x7 onwards: unsure


Got an ssd of 0 which seems to work for all square kernels

  • $\begingroup$ Also, I would like to ask whether rotating the kernel is necessary when performing fast convolution $\endgroup$
    – Toast
    Oct 26, 2022 at 22:11
  • $\begingroup$ So i did remove the rotation of the kernel and I did get a much closer result. The output comparison is the sum of squared differences (subtracting the scipy convolve output and output, raising it to the power of 2 and then adding). I currently have a value of 0.290 when I don't rotate the kernel. $\endgroup$
    – Toast
    Oct 26, 2022 at 23:24
  • $\begingroup$ When using a smaller testing np array of 10x10 and a test kernel of 3x3, looking at the last row, I can see a notable difference in the numbers (the other rows are out by less than 1.) I'm not sure why the bottom row is so different to my fft. I suspect maybe a padding issue. $\endgroup$
    – Toast
    Oct 27, 2022 at 0:12

1 Answer 1

  image_padding = len(kernel)//2
  image_pad = np.pad(image_pad,(((image_padding+1)//2,image_padding//2),((image_padding+1)//2, image_padding//2)), mode='edge') 

Here you are padding by half the amount your really want to pad. You should pad by image_padding on all sides. Furthermore, your input here is image_pad, which hasn't been defined yet. This line should be:

  image_pad = np.pad(image, image_padding, mode='edge')

Then, at the end, you crop those same amounts from all sides:

  return output[image_padding:-image_padding, image_padding:-image_padding]

(By the way, variable names image_pad and image_padding used together is quite confusing. I would name the padded image something like img_padded, padded, or simply image.)

The next issue is the padding of the kernel:

  kernel = np.pad(kernel, (((pad_x+1)//2,pad_x//2),((pad_y+1)//2,pad_y//2)), 'constant')

This is almost right. It works correctly if the input image is odd in size, and/or if the padding is even. But for an even-sized image with odd-sized padding (leading to an odd-sized image), the padding should be larger on the right, not on the left. Remember that the origin here always needs to be at shape//2, so simplest way to correctly compute the padding is as follows:

pad_0_low = image_pad.shape[0] // 2 - kernel.shape[0] // 2
pad_0_high = image_pad.shape[0] - kernel.shape[0] - pad_0_low
pad_1_low = image_pad.shape[1] // 2 - kernel.shape[1] // 2
pad_1_high = image_pad.shape[1] - kernel.shape[1] - pad_1_low
kernel = np.pad(kernel, ((pad_0_low, pad_0_high),(pad_1_low, pad_1_high)), 'constant')

Lastly, remember that the FFT introduces rounding errors, so don't expect an identical result. On an array that is 100x100 pixels, a sum of square differences of 0.523 is an average square difference of 5.23e-05 per pixel (assuming the errors are distributed across the whole image and not, say, concentrated at a row of pixels along the edge of the image). For pixels (I presume) in the range 0-255, this is a very small error.

  • $\begingroup$ Thank you, it works great now! One question I had was does uneven padding of the image matter at all? $\endgroup$
    – Toast
    Oct 27, 2022 at 13:32
  • $\begingroup$ @Bread it doesn’t. You can pad your image however you want to. It just needs to be at least the amount we pad here. The kernel padding has to correct though. That is the “sensitive” one. $\endgroup$ Oct 27, 2022 at 13:41
  • $\begingroup$ @Bread If the kernel is not padded correctly, the result will be shifted. $\endgroup$ Oct 27, 2022 at 13:45
  • $\begingroup$ I see thank you. So the centre of my kernel was not shape//2 but rather off slightly leading to a slightly shifted image. I also had a question about "the padding should be larger on the right, not on the left". Is this meant to be the other way round, where the left has more padding than the right? $\endgroup$
    – Toast
    Oct 27, 2022 at 13:50
  • $\begingroup$ @Bread that has to do with padding an even-sized kernel to an odd size, or an odd-sized kernel to an even size, and maintain the center at size//2. In the first case, you need to pad one more pixel on the right. In the second case you need to pads one more pixel on the left. $\endgroup$ Oct 27, 2022 at 14:01

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.