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.
