Gaussian Image filtering using FFT

Question

For image segmentation I use Difference of Gaussian features using OpenCV's GaussianBlur (ranging from 0.8 to 8.43 with exponential step size of 1.4). My images are of size 4096 x 2160 so this takes quite some time (on one core 8 seconds, which is quite long when processing a video).

Can you give me any suggestions on how to speed things up? Currently I am trying to implement Gaussian filtering in FFT. I have the following code so far:

ftimage = np.fft.fft2(image)
ftimage = np.fft.fftshift(ftimage)
kernel = cv2.getGaussianKernel(11, 3)
kernel = kernel * kernel.T
ftkernel = np.fft.fftshift(np.fft.fft2(kernel, (h, w)))
ftimagep = ftimage * gmask
imagep = np.fft.ifft2(ftimagep)
imageq = cv2.GaussianBlur(image, (11,11), 3))

The problem here is that imagep and imageq are shifted versions of each other. Secondly, as the Fourier of a Gaussian is also a Gaussian, how can I compute ftkernel in a straight-forward manner?

---

To the answer below: I have implemented approximate filtering:

 def approx_g(image, sigma_g, n=5):
      w = np.sqrt(12*sigma_g**2/n + 1)
      wu = np.ceil(w) if np.ceil(w) % 2 == 1 else np.ceil(w)+1
      wl = np.floor(w) if np.floor(w) % 2 == 1 else np.floor(w)-1
      if w == w//1:
          wl -= 2
          wu += 2
      m = round((12*sigma_g**2 - n*wl**2 - 4*n*wl - 3*n) / (-4*wl - 4))
      wl = int(wl)
      wu = int(wu)
      for num in range(0,int(m)):
          image = cv2.blur(image, (wl, wl))
      for num in range(0,int(n-m)):
          image = cv2.blur(image, (wu, wu))
      return image

The L2 pixel difference looks quite good for n=4:

I have also done the speed comparison for different sigma's:

Answer 1

A Gaussian filter can be approximated by a cascade of box (averaging) filters, as described in section II of Fast Almost-Gaussian Filtering. This method requires using the Integral Image, and allows faster application of (near) Gaussian filtering, especially for high blur cases.

The code below demonstrates how one might do this using the steps from the paper linked above.

Let the filter radius be 8.43 as in the question.

sigma_g = 8.43

The number of successive applications of the box filter determines the level of approximation. For this example I set it to 5:

n = 5

First, find the ideal width of the box filter using equation 3:

w = np.sqrt(12*sigma_g**2/n + 1)

As discussed in the paper, using two different size box filters works better. The filters need to be of odd length for symmetry, with lengths that differ by two. The code below takes w and finds the nearest odd integers. (it could probably be written better):

wu = np.ceil(w) if np.ceil(w) % 2 == 1 else np.ceil(w)+1
wl = np.floor(w) if np.floor(w) % 2 == 1 else np.floor(w)-1
if w == w//1:
    wl -= 2
    wu += 2

If n successive applications are desired, then m are performed using the first filter with width wu and (n-m) are performed with the second, with width wl. Equation 5 shows how to calculate m:

m = round((12*sigma_g**2 - n*wl**2 - 4*n*wl - 3*n) / (-4*wl - 4))

Next, the functions to calculate the integral image for both horizontal and vertical:

def integral_image_1d_hor(image):
    ''' Calculated the 1d horizontal integral
    image of an image.'''
    n1, n2 = np.shape(image)
    int_im = np.zeros((n1, n2))
    for row in range(0,n1):
        int_im[row,0] = image[row,0]

    for row in range(0,n1):
        for col in range(1,n2):
            int_im[row,col] = image[row,col] + int_im[row,col-1]

    return int_im


def integral_image_1d_ver(image):
    ''' Calculated the 1d vertical integral
        image of an image.'''
    n1, n2 = np.shape(image)
    int_im = np.zeros((n1, n2))
    for col in range(0,n2):
        int_im[0,col] = image[0,col]

    for col in range(0,n2):
        for row in range(1,n1):
            int_im[row,col] = image[row,col] + int_im[row-1,col]

    return int_im

To filter using the integral images I have these functions:

def box_1d_filter_hor(int_im_1d, width):
    w = int((width-1)/2)
    fil_im = np.zeros(np.shape(int_im_1d))
    pad = w
    int_im_1d = np.pad(int_im_1d, pad, 'constant')
    n1 = np.shape(int_im_1d)[0]
    n2 = np.shape(int_im_1d)[1]
    for row in range(pad, n1-pad):
        for col in range(pad, n2-pad):
            fil_im[row-pad,col-pad] = (int_im_1d[row,col+w]
                                    - int_im_1d[row,col-w-1])/width
    return fil_im


def box_1d_filter_ver(int_im_1d, width):
    w = int((width-1)/2)
    fil_im = np.zeros(np.shape(int_im_1d))
    pad = w
    int_im_1d = np.pad(int_im_1d, pad, 'constant')
    n1 = np.shape(int_im_1d)[0]
    n2 = np.shape(int_im_1d)[1]
    for col in range(pad, n2-pad):
        for row in range(pad, n1-pad):
            fil_im[row-pad,col-pad] = (int_im_1d[row+w,col]
                                    - int_im_1d[row-w-1,col])/width
    return fil_im

Then I define two more functions, for processing the image in the horizontal and vertical directions:

def process_hor(image, w):
    int_im = integral_image_1d_hor(image)
    fil_im = box_1d_filter_hor(int_im, w)
    return fil_im

def process_ver(image, w):
    int_im = integral_image_1d_ver(image)
    fil_im2 = box_1d_filter_ver(int_im, w)
    return fil_im2

Finally, using all these previous functions, to approximate gaussian filtering use the following function:

def approximate_gaussian(image, wl, wu, m, n):
    for num in range(0,int(m)):
        image = process_hor(image, wl)
        image = process_ver(image, wl)
    for num in range(0,int(n-m)):
        image = process_hor(image, wu)
        image = process_ver(image, wu)
    return image

I didn't really handle the edges of the image but that can be adjusted by modifying the functions above. This should be faster, especially for cases where the Gaussian blurring radius is very high.