Why quadrant shift using numpy ifft2 on log-gabor filter

Question

I build a log-gabor filter in the frequency domain following Peter Kovesi's tutorial. My results match the plots in the tutorial as can be seen here:

Results

By my understanding, in order to translate the filter to spatial domain, I have to call ifftshift (to match ifft2 input format) and a consecutive ifft2 (frequency -> spatial domain). However, this leads to a shifted image in spatial domain where quadrants are swapped. Calling an additional fftshift/ifftshift yields the expected image.

log_gabor_filter in spatial domain with filter located at the corners (swapped quadrants)

log_gabor_filter after calling an additional ifftshift in spatial domain

Why does ifft2 not yield the expected output? Why do I have to call an additional fftshift? You find my code below.

import numpy as np
from numpy.fft import ifft2, ifftshift, fftshift
import matplotlib.pyplot as plt

def plot_filter(filter, title):
    plt.clf()
    plt.title(title)
    plt.pcolormesh(filter, cmap='binary_r')
    plt.axis('off')
    plt.show()

def log_gabor_kernel(N, freq0, theta0, log_sigma_freq, sigma_theta):

    # Create Frequency Grid
    center = N//2
    domain = np.arange(-center+1, center+1)/N if N%2 else np.arange(-center, center)/N
    xs, ys = np.meshgrid(domain, domain, indexing='ij')
    freq = np.sqrt(xs**2+ys**2)
    freq[center, center] = 1

    # Calculation of radial factor
    radial = np.exp(-np.log(freq/freq0)**2/(2*log_sigma_freq**2))
    radial[center, center] = 0
    plot_filter(radial, '1) radial component in frequency domain')

    # Calculation of angular factor
    theta = np.arctan2(ys, xs)
    sin_theta = np.sin(theta)
    cos_theta = np.cos(theta)
    # (Using trigonometric addition formulas for angular wrap around problem)
    ds = sin_theta*np.cos(theta0)-cos_theta*np.sin(theta0)
    dc = cos_theta*np.cos(theta0)+sin_theta*np.sin(theta0)
    dtheta = np.abs(np.arctan2(ds, dc))

    angular = np.exp(-dtheta**2/(2*sigma_theta**2))
    plot_filter(angular, '2) angular component in frequency domain')

    # Caclulate log gabor kernel from components
    kernel = radial*angular
    return kernel

if(__name__ == "__main__"):
    kernel_frequency_domain = log_gabor_kernel(N=128, freq0=0.115, theta0=np.pi/4, log_sigma_freq=0.45, sigma_theta=0.65)
    plot_filter(kernel_frequency_domain, '3) log-gabor filter in frequency domain')

    kernel_spatial_domain = ifft2(ifftshift(kernel_frequency_domain))
    plot_filter(kernel_spatial_domain.real, '4i) real part log-gabor filter in spatial domain')
    plot_filter(kernel_spatial_domain.imag, '4ii) imaginary part log-gabor filter in spatial domain')

    kernel_spatial_domain = fftshift(kernel_spatial_domain)
    plot_filter(kernel_spatial_domain.real, '5i) real part log-gabor filter shifted in spatial domain')
    plot_filter(kernel_spatial_domain.imag, '5ii) imaginary part log-gabor filter shifted in spatial domain')

Thanks for reading, any comment is very much appreciated :)