0

I am trying to implement myself a 2D gaussian filter (to debug some C++ code, outside the scope of this question). I focus on the first step: multiplication of the image and the Gaussian kernel in the Fourier space.

The input data and its TF:

import numpy as np
import scipy
sigma = 2
size = 8
np.random.seed(42)
image = np.random.rand(size, size)
imageTf = np.fft.fft2(image)

First version of the multiplication:

imageBluredTfScipy = scipy.ndimage.fourier_gaussian(imageTf, sigma)

Second version with own implementation:

def makeGaussian(size, sigma):
    kernel = np.zeros((size,size))
    sigma2 = sigma**2
    tot = 0.
    for j in range(size):
        y = j - (size - 1.)/2.
        for i in range(size):
            x = i - (size - 1.)/2. 
            kernel[i,j] = np.exp(-(x**2/(2*sigma2) + y**2/(2*sigma2)))
            tot += kernel[i,j]
    kernel /= tot
return kernel

kernel = makeGaussian(size, sigma)
kernelTf = np.fft.fft2(kernel)
imageBluredTf = kernelTf*imageTf

Now, if I compare imageBluredTf and imageBluredTfScipy, it's not the same for both real and imaginary part (expect at [0,0] where both methods give the same value). For this example, the difference of the real part within [:2,:2] is:

array([[-7.10542736e-15,  2.01834288e-01],
   [ 7.12344007e-01,  7.59622995e-02]])

I have searched for the definition of the kernel used by scipy.ndimage.fourier_gaussian without success.

dagnic
  • 137
  • 10
  • Two things are wrong here: `j - (size - 1.)/2.` this needs to be an integer. And the origin of the FFT is on the top-left pixel, see `ifftshift` for how to get your kernel to match. I would suggest putting a `np.ones()` array through `scipy.ndimage.fourier_gaussian` and examining the result. – Cris Luengo Feb 24 '23 at 14:08
  • Another wrong thing: don’t normalize the kernel! The value at the origin must be 1. – Cris Luengo Feb 24 '23 at 14:09
  • Oh, you are FFT-ing the kernel, the above assumed you create the Gaussian directly in frequency domain. Similar considerations apply though, the origin has to be at the top-left pixel. – Cris Luengo Feb 24 '23 at 14:11
  • Hi. I tried `scipy.ndimage.fourier_gaussian(np.ones)`, this is the Gaussian in the frequency space. Indeed, the origin is 1. I reorganized the four quadrants with `ifftshift`, and clearly, the Gaussian is not centred, and 1 is found at [int(size/2), int(size/2)]. Thanks. – dagnic Feb 24 '23 at 17:26
  • Now, I would like to find an expression for the sigma in the frequency space, to avoid the TF of the kernel. I found different expressions from various sources, but no one corresponds to what I can see with `scipy.ndimage.fourier_gaussian(np.ones)`. – dagnic Feb 24 '23 at 18:05
  • 1
    See here: https://stackoverflow.com/a/74539736/7328782 – Cris Luengo Feb 24 '23 at 18:20

0 Answers0