# Efficient implementation of 2-d circularly symmetric low-pass filter

How can a two-dimensional ideal circularly symmetric low-pass filter or its approximation be efficiently implemented on data sampled on a square grid? I'm referring to an ideal filter with a spatial frequency response that equals $$1$$ inside radius $$\omega_c$$ (the cutoff frequency) and zero outside it, and with ideal impulse response:

$$h[x,y] = \frac{\omega_c}{2\pi \sqrt{x^2 + y^2} } J_1 \big( \omega_c \sqrt{x^2 + y^2} \big).$$

A "not efficient" approach is to convolve with a 2-d filter kernel that is the product of $$h[x, y]$$ and a circularly symmetric window function.

• I guess you also expect a computational efficiency due to the circular symmetry. A good deal of patience and circular math should be spent on it to see... May 21, 2019 at 18:42
• @Fat32 Maybe it can be approximated with the real part of a weighted sum of 2-d Gaussian filters with complex variance, each filter separable to a cascade of a vertical and a horizontal 1-d filter. May 21, 2019 at 18:50
• It seems you want to try ;-) So you would benefit from separability... May 21, 2019 at 18:53

# Approximation by the real part of a weighted sum of separable complex Gaussian component kernels

Figure 1. The proposed scheme illustrated as 1-d real convolutions ($$*$$) and additions ($$+$$), for cut-off frequency $$\omega_c = \pi/4$$ and kernel width $$N=41$$. Each the upper and the lower half of the diagram is equivalent to taking the real part of a 1-d horizontal and a 1-d vertical pass separable 2-d complex convolution. Except for truncation/cropping, the two final components summed and the resulting kernel are guaranteed circularly symmetrical due to use of complex-number-weighted complex-variance Gaussian 1-d kernels.

I have suggested this approach for this problem, but not tried it before. A discrete-space 2-d convolution kernel (impulse response) $$h[x, y]$$ that is circularly symmetric (see Eq. 4) about $$[0, 0]$$, can be approximated to arbitrary accuracy by the real part of a separable (see Eq. 3) complex number weighted sum of Gaussian functions of complex variance:

$$h[x, y] \approx \tilde h[x, y],\tag{1}$$

\begin{align}\tilde h[x, y] = &\sum_{m=0}^{M-1}\operatorname{Re}\left(a_m^2e^{b_m(x^2+y^2)}\right)\tag{2}\\ = &\sum_{m=0}^{M-1}\operatorname{Re}\left(a_m e^{b_m x^2}\times a_m e^{b_m y^2}\right)\tag{3}\\ = &\sum_{m=0}^{M-1}\operatorname{Re}\left(a_m^2 e^{b_m r^2}\right)\tag{4}\\ = &\sum_{m=0}^{M-1}e^{\operatorname{Re}(b_m)r^2}\Big({\big(\operatorname{Re}(a_m)^2 - \operatorname{Im}(a_m)^2\big)\cos\big(\operatorname{Im}(b_m)r^2\big)\\ - 2\operatorname{Re}(a_m)\operatorname{Im}(a_m)\sin\big(\operatorname{Im}(b_m)r^2\big)\Big),}\tag{5}\end{align}

where $$x$$ and $$y$$ are integer horizontal and vertical coordinates, $$r = \sqrt{x^2 + y^2}$$ is the radius or distance between $$[0, 0]$$ and $$[x, y]$$, $$M$$ is the approximation order, and $$a_m$$ and $$b_m$$ are complex amplitude and variance related coefficients.

## Optimized approximation, $$\omega_c = \pi$$

Eq. 5 with $$r^2 = x^2 + y^2$$ can be numerically optimized by global optimization to approximate a given $$h[x, y]$$ in spatial domain within the square $$-(N-1)/2 \le x \le (N-1)/2$$, $$-(N-1)/2 \le y \le (N-1)/2$$ with odd width $$N$$ of the square kernel. It suffices to do the optimization over a triangle $$0 \le y \le (N-1)/2$$, $$y \le x \le (N-1)/2$$. If the cost function is a sum of squares of errors or another sum of errors at each $$(x, y)$$, then cumulative $$2\times$$ error multipliers should be applied for each satisfied condition $$x > 0$$, $$y > 0$$, $$x > y$$.

The following C++ program optimizes the approximate kernel $$\tilde h$$ in spatial domain least squares sense using Differential Evolution from an optimization library. The target kernel $$h$$ is windowed using a 2-d rotated cosine window. The source code includes pre-optimized approximations for cut-off frequency $$\omega_c = \pi$$, kernel size $$N=41$$, and approximation orders $$1 \le m \le 5$$.

// -*- compile-command: "g++ -Wno-unused-result -march=native -O3 -ffast-math optitest.cpp opti.cpp" -*-
// This file is optitest.cpp
#include <stdio.h>
#define _USE_MATH_DEFINES
#include <math.h>
#include "opti.hpp"
#include "keyboard.h"

double optimal0[4] = { // Dummy starting point for optimization
0,0,0,0
};

double optimal1[1*4] = {//bestcost=0.03593696209385285195 Least squares
0.90954095944946312,0.21237789645057972,-0.72663606661376057,0.76147736559107793
};

double optimal2[2*4] = {//bestcost=0.01072419992771576382 Least squares
-0.14935266179235707,0.28729650415029556,-0.08227954155725942,0.36475170816661134,1.03225137134955114,0.47158356759095016,-0.60039465413238513,0.47344404338750434
};

double optimal3[3*4] = {//bestcost=0.00322191268986909149 Least squares
-0.00646998371362690,-0.15983965263134517,-0.02321635125852370,0.20289085505437962,-0.30298121364781033,0.34395623806288950,-0.10323939264492392,0.32564190139614663,1.13148495464950427,0.66639187966234981,-0.54567361845465989,0.37519342198291905
};

double optimal4[4*4] = {//bestcost=0.00089888657916611925 Least squares
0.04481214216545543,0.08998210660870602,-0.01130470880820374,0.14003687599487286,0.02506113778123746,-0.24720956687380186,-0.03419309947921907,0.19643391040650712,1.27604032320787875,0.91159589755622883,-0.49585695598458995,0.29480666144138823,0.48259638667856242,-0.42280882153371496,-0.12730997070336811,0.29988730534029784
};

double optimal5[5*4] = {//bestcost=0.00017259634849055045 Least squares
-0.06413980110992069,-0.03272558679644168,-0.00736621171073370,0.10943881870260203,0.08213936888117918,0.16226536153011967,-0.01975126456281387,0.13958603391531316,-0.07270617897425770,0.37473199916953354,-0.04702984945995840,0.18925309227383197,1.43523098843984531,1.17056740170289952,-0.45755771004055446,0.24332707322808175,0.69700596730108921,-0.48854609666142051,-0.14144974773647198,0.27672159791886242
};

const int numOptimalKernels = 5;
double *optimalKernels[numOptimalKernels + 1] = {optimal0, optimal1, optimal2, optimal3, optimal4, optimal5};

class CircularLPProblem : public Opti::Problem {
private:
double *minimum;
double *maximum;
int numComponents;
int numSamples;
int *r2s;
int *multipliers;
double *target;
public:

int getNumDimensions() {
return numComponents*4;
}

double *getMin() {
return minimum;
}

double *getMax() {
return maximum;
}

double costFunction(double *params, double compare) {
for (int m = 0; m < numComponents; m++) {
params[4*m + 2] = -fabs(params[4*m + 2]);
params[4*m + 3] = fmod(fabs(params[4*m + 3]), M_PI);
}
for (int m = 0; m < numComponents - 1; m++) {
if (params[4*m + 3] > params[4*(m + 1) + 3]) {
for (int k = 0; k < 4; k++) {
double temp = params[4*m + k];
params[4*m + k] = params[4*(m + 1) + k];
params[4*(m + 1) + k] = temp;
}
}
}
double cost = 0;
for (int k = 0; k < numSamples; k++) {
double r2 = r2s[k];
double trial = 0;
for (int m = 0; m < numComponents; m++) {
//  trial += exp(params[4*m + 2]*r2)*(params[4*m + 0]*cos(params[4*m + 3]*r2) + params[4*m + 1]*sin(params[4*m + 3]*r2)); // Max absolute error
trial += exp(params[4*m + 2]*r2)*((params[4*m + 0]*params[4*m + 0] - params[4*m + 1]*params[4*m + 1])*cos(params[4*m + 3]*r2) - 2*params[4*m + 0]*params[4*m + 1]*sin(params[4*m + 3]*r2)); // Least squares
}
/*      if (fabs(trial - target[k]) > cost) { // Max absolute error
cost = fabs(trial - target[k]);
}*/
cost += (trial - target[k])*(trial - target[k])*multipliers[k]; // Least squares
if (cost > compare) {
return cost;
}
}
return cost;
}

// numComponents = number of components
// N = kernel width
// omega_c = cutoff frequency (radians)
// Parameter vector: Re(a_0), Im(a_0), Re(b_0), Im(b_0), Re(a_1), ...
CircularLPProblem(int numComponents, int N, double omega_c = M_PI): numComponents(numComponents) {
numSamples = 0;
for (int y = 0; y < (N-1)/2 + 1; y++) {
numSamples += (N-1)/2 + 1 - y;
}
r2s = new int[numSamples];
multipliers = new int[numSamples];
target = new double[numSamples];
int k = 0;
for (int y = 0; y < (N-1)/2 + 1; y++) {
for (int x = y; x < (N-1)/2 + 1; x++) {
r2s[k] = x*x + y*y;
target[k] = omega_c*j1(omega_c*sqrt(x*x + y*y))/(2*M_PI*sqrt(x*x + y*y));
double window = cos(M_PI/2*sqrt(pow(x/((N - 1)/2 + 1.0), 2) + pow(y/((N - 1)/2 + 1.0), 2)));
if (window < 0) {
target[k] = 0;
} else {
target[k] *= window;
}
multipliers[k] = ((x > 0) ? 2 : 1) * ((y > 0) ? 2 : 1) * ((x > y) ? 2 : 1);
k++;
}
}
target[0] = omega_c*omega_c/(4*M_PI);
minimum = new double[4*numComponents];
maximum = new double[4*numComponents];
k = 0;
for (int i = 0; i < 4*numComponents; i++) {
minimum[i] = optimalKernels[numComponents - 1][k]-pow(0.1, numComponents - 1);
maximum[i] = optimalKernels[numComponents - 1][k]+pow(0.1, numComponents - 1);
k++;
if (k >= (numComponents-1)*4) {
k -= 4;
if (k < 0) {
k = 0;
}
}
}
}

~CircularLPProblem() {
delete[] minimum;
delete[] maximum;
delete[] r2s;
delete[] multipliers;
delete[] target;
}
};

int main()
{
INITKEYBOARD;

CircularLPProblem problem(1, 41, M_PI); // Parameterize this!

Opti::Strategy *optimizer;

//optimizer=new Opti::G3(&problem, 2000);
optimizer=new Opti::DE(&problem, 2000);

printf("\nOptimizing...\n\n");

for(int t = 0;; t++) {
double bestcost = optimizer->evolve();
if (!(t % 1000)) {
printf("gen=%d, bestcost=%.20f, average=%.20f\n", t, bestcost, optimizer->averageCost());
if (kbhit()) {
printf("Parameter vector printout:\n");
problem.print(optimizer->best());
if (getch() == 27) break;
getch();
}
}
}
delete optimizer;

DEINITKEYBOARD;
return 0;
}


Figure 2. Sum-of-squares cost of the best found approximate kernel $$\tilde h$$ as function of approximation order $$M$$, for $$\omega_c = \pi$$, $$\omega_c = \pi/2$$, and $$\omega_c = \pi/4$$, and $$N=41$$. A good indicator that these are the globally optimal parameterizations is the steady decrease of error as $$M$$ is incremented, except for $$\omega_c = \pi/4$$, $$M=3$$ for which we might (or not) have a sub-optimal solution.

This Python script generates the "exact" (windowed) kernel for comparison and implements the approximations:

import matplotlib.pyplot as plt
from scipy import special
import numpy as np
import scipy.ndimage
import skimage
import time

def circularLowpassKernel(omega_c, N):  # omega = cutoff frequency in radians (pi is max), N = horizontal size of the kernel, also its vertical size, must be odd.
kernel = np.fromfunction(lambda x, y: omega_c*special.j1(omega_c*np.sqrt((x - (N - 1)/2)**2 + (y - (N - 1)/2)**2))/(2*np.pi*np.sqrt((x - (N - 1)/2)**2 + (y - (N - 1)/2)**2)), [N, N])
kernel[(N - 1)//2, (N - 1)//2] = omega_c**2/(4*np.pi)
return kernel

def rotatedCosineWindow(N):  # N = horizontal size of the targeted kernel, also its vertical size, must be odd.
return np.fromfunction(lambda y, x: np.maximum(np.cos(np.pi/2*np.sqrt(((x - (N - 1)/2)/((N - 1)/2 + 1))**2 + ((y - (N - 1)/2)/((N - 1)/2 + 1))**2)), 0), [N, N])

N = 41  # Kernel width
M = 1  # Number of components
omega_c = np.pi  # Cutoff frequency <= np.pi
kernel = circularLowpassKernel(omega_c, N)*rotatedCosineWindow(N)

def saveKernel(name, kernel):
plt.imsave(name+'.png', plt.cm.bwr(plt.Normalize(vmin=-kernel.max(), vmax=kernel.max())(skimage.transform.rescale(kernel, 4, 0))))
absF = np.abs(np.fft.fftshift(np.fft.fft2(kernel)))
plt.imsave(name+'_dft.png', plt.cm.Greys(plt.Normalize(vmin=0, vmax=absF.max())(skimage.transform.rescale(absF, 4, 0))))

saveKernel('exactpi', kernel)

plt.imsave('bwrkey.png', plt.cm.bwr(np.repeat([(np.arange(41*4)/(41*4-1))], 16, 0)))
plt.imsave('Greyskey.png', plt.cm.Greys(np.repeat([(np.arange(41*4)/(41*4-1))], 16, 0)))

def reComponentKernel(N, param):
x = np.arange(N)-(N - 1)//2
return np.exp(param[2]*x**2)*(param[0]*np.cos(param[3]*x**2) - param[1]*np.sin(param[3]*x**2))

def imComponentKernel(N, param):
x = np.arange(N)-(N - 1)//2
return np.exp(param[2]*x**2)*(param[1]*np.cos(param[3]*x**2) + param[0]*np.sin(param[3]*x**2))

optimal1pi = [0.90954095944946312,0.21237789645057972,-0.72663606661376057,0.76147736559107793]
optimal2pi = [-0.14935266179235707,0.28729650415029556,-0.08227954155725942,0.36475170816661134,1.03225137134955114,0.47158356759095016,-0.60039465413238513,0.47344404338750434]
optimal3pi = [-0.00646998371362690,-0.15983965263134517,-0.02321635125852370,0.20289085505437962,-0.30298121364781033,0.34395623806288950,-0.10323939264492392,0.32564190139614663,1.13148495464950427,0.66639187966234981,-0.54567361845465989,0.37519342198291905]
optimal4pi = [0.04481214216545543,0.08998210660870602,-0.01130470880820374,0.14003687599487286,0.02506113778123746,-0.24720956687380186,-0.03419309947921907,0.19643391040650712,1.27604032320787875,0.91159589755622883,-0.49585695598458995,0.29480666144138823,0.48259638667856242,-0.42280882153371496,-0.12730997070336811,0.29988730534029784]
optimal5pi = [-0.06413980110992069,-0.03272558679644168,-0.00736621171073370,0.10943881870260203,0.08213936888117918,0.16226536153011967,-0.01975126456281387,0.13958603391531316,-0.07270617897425770,0.37473199916953354,-0.04702984945995840,0.18925309227383197,1.43523098843984531,1.17056740170289952,-0.45755771004055446,0.24332707322808175,0.69700596730108921,-0.48854609666142051,-0.14144974773647198,0.27672159791886242]

class SeparableCircularLowpassFilter:
def __init__(self, N, coefs):
self.N = N
self.reKernels = []
self.imKernels = []
for i in range(len(coefs)//4):
self.reKernels.append(np.array([reComponentKernel(N, coefs[i*4:])]))
self.imKernels.append(np.array([imComponentKernel(N, coefs[i*4:])]))
def filter(self, x):
reZ = scipy.ndimage.convolve(scipy.ndimage.convolve(x, self.reKernels[0].transpose()), self.reKernels[0]) - scipy.ndimage.convolve(scipy.ndimage.convolve(x, self.imKernels[0].transpose()), self.imKernels[0])
for i in range(1, len(self.reKernels)):
reZ += scipy.ndimage.convolve(scipy.ndimage.convolve(x, self.reKernels[i].transpose()), self.reKernels[i]) - scipy.ndimage.convolve(scipy.ndimage.convolve(x, self.imKernels[i].transpose()), self.imKernels[i])
return reZ

filter1pi = SeparableCircularLowpassFilter(N, optimal1pi) # 1 component kernel
filter2pi = SeparableCircularLowpassFilter(N, optimal2pi) # 2 component kernel
filter3pi = SeparableCircularLowpassFilter(N, optimal3pi) # 3 component kernel
filter4pi = SeparableCircularLowpassFilter(N, optimal4pi) # 4 component kernel
filter5pi = SeparableCircularLowpassFilter(N, optimal5pi) # 5 component kernel

x = np.zeros([N, N])  # Input image
x[N//2, N//2] = 1     # (unit impulse)
# x = plt.imread('sample.tif').astype(float)[:,:,1]  # (green channel of some image from file)

t0 = time.time()
Z = scipy.ndimage.convolve(x, kernel)  # Exact
t1 = time.time()
print(t1-t0)

t0 = time.time()
reZ1pi = filter1pi.filter(x)  # 1 component kernel
t1 = time.time()
print(t1-t0)
saveKernel('reZ1pi', reZ1pi)

t0 = time.time()
reZ2pi = filter2pi.filter(x)  # 2 component kernel
t1 = time.time()
print(t1-t0)
saveKernel('reZ2pi', reZ2pi)

t0 = time.time()
reZ3pi = filter3pi.filter(x)  # 3 component kernel
t1 = time.time()
print(t1-t0)
saveKernel('reZ3pi', reZ3pi)

t0 = time.time()
reZ4pi = filter4pi.filter(x)  # 4 component kernel
t1 = time.time()
print(t1-t0)
saveKernel('reZ4pi', reZ4pi)

t0 = time.time()
reZ5pi = filter5pi.filter(x)  # 5 component kernel
t1 = time.time()
print(t1-t0)
saveKernel('reZ5pi', reZ5pi)


$$M=1$$
$$M=2$$
$$M=3$$
$$M=4$$
$$M=5$$

exact
Figure 3. Kernels and the absolute value of their discrete Fourier transform (DFT), enlarged by a factor of 4 to make individual pixels visible. Top to bottom: Optimal approximating kernels $$\tilde h$$ with $$M = 1$$, $$2$$, $$3$$, $$4$$ and $$5$$ complex separable components, ideal kernel $$h$$ for $$\omega_c = \pi$$ and $$N = 41$$. Color keys: kernel: blue negative, white zero, red positive (normalized); abs DFT: white zero, black maximum.

## Optimized approximation, $$\omega_c = \pi/2$$

Approximation of lower-cutoff kernels works better with this approach. I think this is because the number of ripples that fit in the windowed target kernel are reduced. For a high-cutoff filter then perhaps the kernel width $$N$$ could be reduced, as this would give the same number of ripples and frequency-domain approximation accuracy. But that would then favor a direct 2-d implementation (see Fig. 6). The C++ source code in the above was modified (not shown) for $$\omega_c = \pi/2$$ (and later for $$\omega_c = \pi/4$$) and the obtained parameters were used in a Python implementations of the approximate isotropic low-pass filters: (continued from the previous script)

omega_c = np.pi/2
kernelpi2 = circularLowpassKernel(omega_c, N)*rotatedCosineWindow(N)
saveKernel('exactpi2', kernelpi2)

optimal1pi2 = [0.44103810622146067,0.08998875769710178,-0.17010258583392401,0.19960767673288432]
optimal2pi2 = [-0.07233719880423649,0.14289630144713414,-0.02205699413927855,0.09165233018125875,0.51013245392078410,0.22632987351129516,-0.15215724343836151,0.12131467270512424]
optimal3pi2 = [-0.03867959516035375,-0.07754823299868645,-0.00905504286350691,0.05573851697983074,0.60460672854618647,0.36747347219951876,-0.13550005930232881,0.08061493799161984,-0.15528642640407436,0.23307152837452039,-0.03561076440525033,0.08259420496099962]
optimal4pi2 = [0.06323093028956613,0.02993129946018375,-0.00937470404526276,0.03826932634049100,0.78321346985185014,0.58823521191007977,-0.11611698422922974,0.05181454568143690,-0.07541120991444331,-0.16653642294319324,-0.01587217217852400,0.05400812761661748,-0.30454011473801174,0.38069299325271988,-0.04600795167685135,0.07149320869042795]

filter1pi2 = SeparableCircularLowpassFilter(N, optimal1pi2) # 1 component kernel
filter2pi2 = SeparableCircularLowpassFilter(N, optimal2pi2) # 2 component kernel
filter3pi2 = SeparableCircularLowpassFilter(N, optimal3pi2) # 3 component kernel
filter4pi2 = SeparableCircularLowpassFilter(N, optimal4pi2) # 4 component kernel

reZ1pi2 = filter1pi2.filter(x)
saveKernel('reZ1pi2', reZ1pi2)
reZ2pi2 = filter2pi2.filter(x)
saveKernel('reZ2pi2', reZ2pi2)
reZ3pi2 = filter3pi2.filter(x)
saveKernel('reZ3pi2', reZ3pi2)
reZ4pi2 = filter4pi2.filter(x)
saveKernel('reZ4pi2', reZ4pi2)


$$M=1$$
$$M=2$$
$$M=3$$
$$M=4$$

exact
Figure 4. Kernels and the absolute value of their discrete Fourier transform (DFT), enlarged by a factor of 4 to make individual pixels visible. Top to bottom: Optimal approximating kernels $$\tilde h$$ with $$M = 1$$, $$2$$, $$3$$, and $$4$$ complex separable components, ideal kernel $$h$$ for $$\omega_c = \pi/2$$ and $$N = 41$$. Color keys: kernel: blue negative, white zero, red positive (normalized); abs DFT: white zero, black maximum.

## Optimized approximation, $$\omega_c = \pi/4$$

In Python: (continued from the first Python script)

omega_c = np.pi/4
kernelpi4 = circularLowpassKernel(omega_c, N)*rotatedCosineWindow(N)
saveKernel('exactpi4', kernelpi4)

optimal1pi4 = [0.22210250507512549,0.04738076096442791,-0.04565916667208288,0.04924268271256500]
optimal2pi4 = [0.28445762468654334,0.14050386040493884,-0.03988634664177378,0.02123143861749057,-0.01717623293554685,0.11158938041615812,-0.01191733840600387,0.02614013272136841]
optimal3pi4 = [-0.00332428227606697,0.01722626235358831,-0.00516799783527921,0.00727825629264402,0.30296628929626013,0.16328004086794692,-0.03813290200362767,0.01805103148816284,-0.02296006483369628,0.12768470633404483,-0.01299438181561595,0.02558273151301529]

filter1pi4 = SeparableCircularLowpassFilter(N, optimal1pi4) # 1 component kernel
filter2pi4 = SeparableCircularLowpassFilter(N, optimal2pi4) # 2 component kernel
filter3pi4 = SeparableCircularLowpassFilter(N, optimal3pi4) # 3 component kernel

reZ1pi4 = filter1pi4.filter(x)
saveKernel('reZ1pi4', reZ1pi4)
reZ2pi4 = filter2pi4.filter(x)
saveKernel('reZ2pi4', reZ2pi4)
reZ3pi4 = filter3pi4.filter(x)
saveKernel('reZ3pi4', reZ3pi4)


$$M=1$$
$$M=2$$
$$M=3$$

exact
Figure 5. Kernels and the absolute value of their discrete Fourier transform (DFT), enlarged by a factor of 4 to make individual pixels visible. Top to bottom: Optimal approximating kernels $$\tilde h$$ with $$M = 1$$, $$2$$ and $$3$$ complex separable components, ideal kernel $$h$$ for $$\omega_c = \pi/4$$ and $$N = 41$$. Color keys: kernel: blue negative, white zero, red positive (normalized); abs DFT: white zero, black maximum.

## Computational complexity

Complexity as the number of 1-d dot products of length $$N$$ per pixel is $$4M$$ for the approximations and $$N$$ for the naive 2-d convolution implementation. The approximations suffer from an additional performance overhead due to the use of additional memory for storing intermediate results.

Figure 6. Complexity as the number of real multiplications per pixel vs. kernel width $$N$$. In reality the approximations are not quite as efficient, see Fig. 7. Possible advantages from symmetry of the kernels has not been taken into account.

Figure 7. Actual measured run time in seconds for 2-d convolution of a 5424×3636-pixel single-channel image by the ideal kernel $$h$$ and by separable approximating kernels $$\tilde h$$ with $$M = 1$$, $$2$$, $$3$$, $$4$$ and $$5$$ complex separable components, for $$N = 41$$, using Python's scipy.ndimage.convolve.

## Further ideas

In order to adjust $$\omega_c$$ without kernel parameter re-optimization, a nearly continuous approximating kernel could be optimized by using a small $$\omega_c$$ and a large $$N$$. The result could then be sampled using interpolation.

It would be interesting to try direct optimization of the complex kernels without a parameterization that enforces circular symmetry. Perhaps the approximation could even be optimized as a sum of $$M$$ real separable kernels, similar to this Mathematics Stack Exchange question. Also minimization of (least squares) error in the frequency domain is something to try.

• Ok. That's another great answer from Olli :-) ! Jun 18, 2019 at 20:11
• Nevertheless, I must confess the disappointment I had when I saw the numerical optimization for the approximation of the circular impulse response by a sum of complex weighted complex variance Gaussian kernels. I was expecting some analytical expressions for computing those complex weights and complex variances ;-) apart from kidding, the effort went into the numerical method is also breathtaking ! Jun 18, 2019 at 20:14
• The only objection I would have would be on the justification of using complex weighted complex variance Gaussian kernels: and then taking their real parts anyway? Why cannot you directly start with plain real weighted real Gaussians; may be with more real parameters though. Furthermore, in your justificaiton of the circularity of the approximation, i think it's only about the term x^2+y^2 and not about the complexness of the approximating kernels? Am I right ? Would the approximation still be circularly symmetric for a complex-kernel of $b \cdot e^{a(x+y)}$, for complex $a,b$ ? Jun 18, 2019 at 20:21
• @Fat32 According to my tests, you would not get convergence of the approximation with real-weighted real Gaussians. They do give circular symmetry though. $b \cdot e^{a(x+y)}$ does not give circular symmetry. Jun 19, 2019 at 7:19
• Hmm, the following type of kernels $a e^{b(x^2+y^2)} \cos(b(x^2+y^2))$, $a e^{b(x^2+y^2)} \sin(b(x^2+y^2))$ or their linear combination with real $a,b$ could possibly be converging... Jun 19, 2019 at 8:26