5
$\begingroup$

For a square pixel grid, the ideal 2-d low-pass filter with a horizontal and a vertical cut-off angular frequency $\omega_c$ in radians has an impulse response (kernel) $h_{\small\square}(x, y)$ that is the product of a horizontal and a vertical stretched and scaled sinc function, with $x$ and $y$ the integer horizontal and vertical pixel coordinates:

$$h_{\small\square}[x, y] = \frac{\omega_c\operatorname{sinc}\left(\frac{\omega_c x}{\pi}\right)}{\pi}\times\frac{\omega_c\operatorname{sinc}\left(\frac{\omega_c y}{\pi}\right)}{\pi} = \begin{cases}\frac{\omega_c^2}{\pi^2}&\text{if }x = y = 0,\\\frac{\sin(\omega_c x)\sin(\omega_c y)}{\pi^2 x y}&\text{otherwise.}\end{cases}\tag{1}$$

If $\omega_c = \pi$, the kernel of Eq. 1 is simply:

$$h_{\small\square} = [1]\quad \text{if}\quad\omega_c = \pi.\tag{2}$$

For real-valued $x$ and $y$, Eq. 1 would not describe a circularly symmetric kernel, and neither is the square-like frequency response $H_{\small\square}$ circularly symmetric. Sometimes an isotropic kernel may be desired, for example if we want to remove the effect of the choice of the direction of the image axes.

What is a kernel $h_\circ[x, y]$ such that the frequency response $H_{\circ}$ is $1$ inside a circle of radius $\omega_c$ (angular frequency in radians) and $0$ outside it?

$\endgroup$
3
$\begingroup$

Excerpted from Jae S.Lim 2D signal and image processing ch.1, as an example of $2$-D circularly symmetric lowpass filter with a cutoff frequency of $\omega_c$ radians per sample, whose impulse response is given by: $$h[n_1,n_2] = \frac{\omega_c}{2\pi \sqrt{n_1^2 + n_2^2} } J_1 \big( \omega_c \sqrt{n_1^2 + n_2^2} \big) $$

where $J_1$ is the Bessel function of the first kind and the first order...

Interested readers may consult the book for a derivation which is not untitive but nevertheless tractable; familiarity with Bessel functions is required but also provided as is.[Derivation is added.]

[Olli below]

At $n_1 = n_2 = 0$ the limiting value must be used:

$$h[0, 0] = \frac{\omega_c^2}{4\pi}$$

A slice through the middle of $h[n_1,n_2]$ with $\omega_c = \pi$:

enter image description here

Python source for 2-d filter kernel (you may want to apply a 2-d window function):

from scipy import special
import numpy as np

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

Example with $\omega_c = \pi$:

import matplotlib.pyplot as plt

kernelN = 11  # Horizontal size of the kernel, also its vertical size. Must be odd.
omega_c = np.pi  # Cutoff frequency in radians <= pi
kernel = circularLowpassKernel(omega_c, kernelN)
plt.imshow(kernel, vmin=-1, vmax=1, cmap='bwr')
plt.colorbar()
plt.show()

enter image description here

Example with $\omega_c = \pi/4$:

kernelN = 41  # Horizontal size of the kernel, also its vertical size. Must be odd.
omega_c = np.pi/4  # Cutoff frequency in radians <= pi
kernel = circularLowpassKernel(omega_c, kernelN)
plt.imshow(kernel, vmin=-np.max(kernel), vmax=np.max(kernel), cmap='bwr')
plt.colorbar()
plt.show()

enter image description here

[fat32 below]

Ok, for those who would benefit from a (not so intuitive) derivation, here I make a (n almost) verbatim copy of it from the same book.

First, lets simply write the 2D discrete-time inverse Fourier transform to define the impulse response as:

$$ h[n_1, n_2] = \frac{1}{(2\pi)^2} {\int \int}_{\omega_1^2+\omega_2^2< w_c^2} 1 \cdot e^{j(\omega_1 n_1 + \omega_2 n_2)} d\omega_1 d\omega_2 \tag{1} $$

Let's make a change of variables $\omega_1 = r \cos(\theta)$ and $\omega_2 = r \sin(\theta)$ (effectively write the integral in (1) in the cyclindrical (or polar) coordinates) :

$$ h[n_1, n_2] = \frac{1}{(2\pi)^2} \int_{r=0}^{\omega_c} \int_{\theta=a}^{a+2\pi} e^{j r (\cos(\theta) n_1 + \sin(\theta) n_2)} r ~ dr d\theta \tag{2} $$

Now, make a further change of variables $n_1 = n \cos(\phi)$ and $n_2 = n \sin(\phi)$, with $n = \sqrt{ n_1^2 + n_2^2 }$ and obtain

$$ h[n_1, n_2] = \frac{1}{(2\pi)^2} \int_{r=0}^{\omega_c} r ~dr \int_{\theta=a}^{a+2\pi} e^{j r n \big( \cos(\theta) \cos(\phi) + \sin(\theta) \sin(\phi) \big) } d\theta \tag{3} \\$$

which is now: $$ h[n_1, n_2] = \frac{1}{(2\pi)^2} \int_{r=0}^{\omega_c} r ~dr \int_{\theta=a}^{a+2\pi} e^{j r n \cos(\theta -\phi) } d\theta \tag{4} \\$$

defining $f(r) = \int_{\theta=a}^{a+2\pi} e^{ j r n \cos(\theta-\phi) } d\theta $ , then we get:

$$ h[n_1, n_2] = \frac{1}{(2\pi)^2} \int_{r=0}^{\omega_c} r f(r) dr \tag{5} $$

Now, exploring $f(r)$ we can apply the Euler's identity :

$$ f(r) = \int_{\theta=a}^{a+2\pi} \cos(r.n.\cos(\theta-\phi)) d\theta + j \int_{\theta=a}^{a+2\pi} \sin(r.n.\cos(\theta-\phi) ) d\theta \tag{6} $$

And we can notice that the imaginary part is zero (can check an integral table for that) and with $a = \phi$ then $f(r)$ becomes

$$f(r) = \int_{\theta=0}^{2\pi} \cos(r.n.\cos(\theta) ) d\theta \tag{7}$$

Now the integral in (7) is recognized as the the Bessel function of order zero, kind one $J_0(x)$ which is given as:

$$ J_0(x) = \frac{1}{2\pi} \int_{\theta =0}^{2\pi} \cos( x \cos( \theta) ) d\theta \tag{8} $$

from (7) and (8) we see that $f(r) = 2\pi J_0(r n) $...

And the last identity is given as is: $$ x J_1(x) |_a^b = \int_a^b x J_0(x) dx \tag{9}$$

Now putting the impulse response into the form

$$h[n_1,n_2] = \frac{1}{(2\pi)^2} \int_{r=0}^{\omega_c} r 2\pi J_0( r n) dr \tag{10}$$

applying (9) onto (10) with the recognition that $x = r n$ and $dr = dx/n$ and $n= \sqrt{n_1^2+n_2^2}$ yields the result:

$$\boxed{ h[n_1,n_2] = \frac{\omega_c}{2\pi \sqrt{n_1^2+n_2^2}} J_1( \omega_c \sqrt{n_1^2+n_2^2}) } \tag{11}$$

$\endgroup$
  • $\begingroup$ @OlliNiemitalo I see that you've added an efficient implentation for the (approximately) circular 2D filter. I will try to have a detailed look at it, when I can find some sufficient time... But it already seems great ! ;-) $\endgroup$ – Fat32 May 31 at 20:00
  • $\begingroup$ yes, here, and it it seems to work very well for lower cut-off frequencies. :-) $\endgroup$ – Olli Niemitalo Jun 6 at 16:51

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.