# Linear Filtering¶

Faisal Qureshi
Professor
Faculty of Science
Ontario Tech University
http://vclab.science.ontariotechu.ca

## Outline¶

• Linear Filtering in 1D
• Cross-correlation
• Convolution
• Gaussian filter
• Gaussian blurring
• Separability
• Relationship to Fourier transform
• Integral images
In [1]:
import numpy as np
import scipy as sp
from scipy import signal
import cv2 as cv
import matplotlib.pyplot as plt

In [2]:
# Load image of Van Gogh, convert from BGR to RGB image
I = cv.cvtColor(I, cv.COLOR_BGR2RGB)

In [3]:
plt.figure(figsize=(7,7))
plt.imshow(I)
plt.xticks([])
plt.yticks([]);


Let's plot a single row of the red channel of image $I$, say, row 200.

In [4]:
row = 200
rowI = I[row,:,0]
plt.figure()
plt.plot(rowI,'-')
plt.plot(rowI,'r.');


## Using kernel to compute mean¶

In [5]:
signal = rowI[0:5]
kernel = np.array([1,1,1,1,1]) / 5

print(rowI[0:5])
print('A = ', np.mean(rowI[0:5]))
print('B = ', np.dot(kernel, signal))

[133 127 120 119 125]
A =  124.8
B =  124.8


## Performing moving averages¶

Let's perform a moving average on this row. We assume window with half-width equal to $w$.

In [6]:
width = 5
signal = rowI
kernel = np.ones(2*width+1)/(2*width+1) # Array to represent the window (constant-valued)
result = np.ones(len(signal)-2*width)   # Array to store computed moving averages
for i in range(len(result)):
centre = i + width
result[i] = np.dot(kernel, signal[centre-width:centre+width+1])

In [7]:
# Generate a plot of the original 1D signal and its average
plt.figure()
plt.plot(signal,'r.', alpha=0.3)
plt.plot(signal,'-', alpha=0.3)
plt.plot(np.arange(width, len(signal)-width),result,'k-');


## Linear Filtering in 1D¶

• Signal: $\mathbf{f}$
• Kernel (also sometimes called a filter or a mask): $\mathbf{h}$
• Half-width of the kernel: $w$ (i.e., kernel width is $2 w + 1$)

### Cross-Correlation¶

$$CC(i) = \sum_{k \in [-w,w]} \mathbf{f}(i + k) \mathbf{h}(k)$$

### Convolution¶

$$(\mathbf{f} \ast \mathbf{k} )_i = \sum_{k \in [-w,w]} \mathbf{f}(i - k) \mathbf{h}(k)$$

### Flipping kernel¶

Convolution and cross-correlation give the same results for symmetric kernels

In [8]:
kernel = np.array([3,1,2,1,3])
flipped_kernel = np.flip(kernel)
print('Symmetric kernel')
print(f'kernel={kernel}, flipped={flipped_kernel}')

kernel = np.array([1,1,2,2,3])
flipped_kernel = np.flip(kernel)
print('Asymmetric kernel')
print(f'kernel={kernel}, flipped={flipped_kernel}')

Symmetric kernel
kernel=[3 1 2 1 3], flipped=[3 1 2 1 3]
Asymmetric kernel
kernel=[1 1 2 2 3], flipped=[3 2 2 1 1]


#### Exercise: computing a 1D convolution (from scratch)¶

Compute $\mathbf{f} \ast \mathbf{h}$ given

$$\mathbf{f} = \left[ \begin{array}{cccccccc} 1 & 3 & 4 & 1 & 10 & 3 & 0 & 1 \end{array} \right]$$

and

$$\mathbf{h} = \left[ \begin{array}{ccc} 1 & 0 & -1 \end{array}\right]$$
In [9]:
# Your solution here

In [10]:
# %load solutions/linear-filtering/solution_01.py
# Computing a 1D convolution from scratch
f = np.array([1,3,4,1,10,3,0,1])
h = np.array([1,0,-1])
width = 1
result = np.ones(len(f)-2*width)   # Array to store computed moving averages
for i in range(len(result)):
centre = i + width
result[i] = np.dot(h[::-1], f[centre-width:centre+width+1])
print(f'signal f:\t\t{f}')
print(f'kernel h:\t\t{h}')
print(f'convolution (f*h):\t{result}')

signal f:		[ 1  3  4  1 10  3  0  1]
kernel h:		[ 1  0 -1]
convolution (f*h):	[  3.  -2.   6.   2. -10.  -2.]


### Built-in Functions for 1D Convolution¶

NumPy, of course, has a builtin function numpy.convolve for convolution that wraps the slow loop in compiled code. Notice that, with the definition above, if the original array has length $N$ and the convolution kernel has width $2w+1$, the convolved array has length $N-2w$. In numpy.convolve, this behaviour is enforced using the keyword argument mode='valid'.

#### Exercise: computing a 1D convolution (with convolve)¶

Re-compute $\mathbf{f} \ast \mathbf{h}$ using numpy.convolve given

$$\mathbf{f} = \left[ \begin{array}{cccccccc} 1 & 3 & 4 & 1 & 10 & 3 & 0 & 1 \end{array} \right]$$

and

$$\mathbf{h} = \left[ \begin{array}{ccc} 1 & 0 & -1 \end{array}\right].$$
In [11]:
# Your solution here

In [12]:
# %load solutions/linear-filtering/solution_02.py
# Computing a 1D convolution using np.convolve
f = np.array([1, 3, 4, 1, 10, 3, 0, 1])
h = np.array([1, 0, -1])
r = np.convolve(f, h, mode='valid')
print(f'signal f:\t\t{f}')
print(f'kernel h:\t\t{h}')
print(f'convolution (f*h):\t{r}')

signal f:		[ 1  3  4  1 10  3  0  1]
kernel h:		[ 1  0 -1]
convolution (f*h):	[  3  -2   6   2 -10  -2]


### Example: Gaussian smoothing for noise removal¶

• Construct a 1D (noisy) signal
• Filter this signal with a 1D Gaussian kernel to suppress the noise and improve the signal-to-noise ratio
In [13]:
np.random.seed(0)
n = 20
x = np.linspace(0, np.pi, n)
y = np.sin(x)
mu = 0.0
sigma = 0.2
y_noisy = y + np.random.normal(mu, sigma, n)

plt.plot(x,y,'.-', label='data', alpha=0.5)
plt.plot(x,y_noisy,'.-r', label='noisy data')
plt.xlim(-.4,3.5)
plt.ylim(-1,2)
plt.xticks([])
plt.yticks([])
plt.legend(['data', 'noisy data']);


#### Exercise: lets try to recover original signal (data) from corrupted signal (noisy data) using smoothing¶

In [14]:
# %load solutions/linear-filtering/solution_03.py
f = y_noisy
width = 1
h = np.ones(2*width+1)/(2*width+1)
result = np.ones(len(f)-2*width)   # Array to store computed moving averages
for i in range(len(result)):
centre = i + width
result[i] = np.dot(h[::-1], f[centre-width:centre+width+1])

plt.plot(x,y,'.-', label='data', alpha=0.3)
plt.plot(x,y_noisy,'.-r', label='noisy data', alpha=0.5)
plt.plot(x[width:-width],result,'.-b', label='noisy data')
plt.xlim(-.4,3.5)
plt.ylim(-1,2)
plt.xticks([])
plt.yticks([])
plt.legend(['data', 'noisy data']);


Note that lengths of x and result are not the same?

## Gaussian in 1D¶

$$G(\mu,\sigma,x) = \frac{1}{\sigma \sqrt{2 \pi}} \exp ^ {- \frac{(x-\mu)^2}{2 \sigma^2}}$$

Here $\mu$ and $\sigma$ refer to the mean and standard deviation of this Gaussian.

## Aside: Mean and Standard Deviation¶

Given $n$ data points $\{x_1, x_2, \cdots, x_n \}$

$$\mu = \mathrm{E}[x] = \frac{1}{n} \sum_{i=1}^n x_i$$$$\sigma^2 = \mathrm{E}[(x-\mu)^2] = \frac{1}{n} \sum_{i=1}^n (x_i - \mu)^2$$

## Evaluating Gaussian at a particular value.¶

In [15]:
mu = 2
sigma = 3

In [16]:
def g(mu, s, x):
exponent = -((x - mu)**2)/(2*s*s)
constant = 1/(s * np.sqrt(2 * np.pi))
return constant * np.exp( exponent )

In [17]:
g_ = np.empty(20)
for i in range(-10,10):
g_[i-10] = g(mu, sigma, i)
print(f'Gaussian with mu={mu} and sigma={sigma} evaluated at {i} is {g_[i-10]}')

Gaussian with mu=2 and sigma=3 evaluated at -10 is 4.461007525496179e-05
Gaussian with mu=2 and sigma=3 evaluated at -9 is 0.0001600902172069401
Gaussian with mu=2 and sigma=3 evaluated at -8 is 0.0005140929987637022
Gaussian with mu=2 and sigma=3 evaluated at -7 is 0.001477282803979336
Gaussian with mu=2 and sigma=3 evaluated at -6 is 0.003798662007932481
Gaussian with mu=2 and sigma=3 evaluated at -5 is 0.008740629697903166
Gaussian with mu=2 and sigma=3 evaluated at -4 is 0.017996988837729353
Gaussian with mu=2 and sigma=3 evaluated at -3 is 0.03315904626424957
Gaussian with mu=2 and sigma=3 evaluated at -2 is 0.05467002489199788
Gaussian with mu=2 and sigma=3 evaluated at -1 is 0.0806569081730478
Gaussian with mu=2 and sigma=3 evaluated at 0 is 0.10648266850745075
Gaussian with mu=2 and sigma=3 evaluated at 1 is 0.12579440923099774
Gaussian with mu=2 and sigma=3 evaluated at 2 is 0.1329807601338109
Gaussian with mu=2 and sigma=3 evaluated at 3 is 0.12579440923099774
Gaussian with mu=2 and sigma=3 evaluated at 4 is 0.10648266850745075
Gaussian with mu=2 and sigma=3 evaluated at 5 is 0.0806569081730478
Gaussian with mu=2 and sigma=3 evaluated at 6 is 0.05467002489199788
Gaussian with mu=2 and sigma=3 evaluated at 7 is 0.03315904626424957
Gaussian with mu=2 and sigma=3 evaluated at 8 is 0.017996988837729353
Gaussian with mu=2 and sigma=3 evaluated at 9 is 0.008740629697903166

In [18]:
plt.plot(g_)

Out[18]:
[<matplotlib.lines.Line2D at 0x139ee7e50>]

Vectorized code for evaluating Gaussian over a range of values

In [19]:
def gaussian1d(mu, sig, n, normalized=True):
s = n//2
x = np.linspace(-s,s,n)
exponent_factor = np.exp(-np.power(x - mu, 2.) / (2 * np.power(sig, 2.)))
normalizing_factor = 1. if not normalized else 1./(sig * np.sqrt(2*np.pi))
return normalizing_factor * exponent_factor

In [20]:
n = 9
width = n // 2
sigma = 1.4

g1 = gaussian1d(0.0, sigma, n, normalized=True)
g2 = gaussian1d(0.0, sigma, n, normalized=False)

plt.plot(np.linspace(-width,width,n), g1, 'r-.', label='Normalized')
plt.plot(np.linspace(-width,width,n), g2, 'b-o', label='Not normalized')
plt.title('Gaussian kernels')
plt.legend(loc='center left');

In [21]:
# %load solutions/linear-filtering/solution_05.py
y_smoothed = np.convolve(y_noisy, g1, mode='valid')

plt.plot(x,y,'b', label='data', alpha=0.3)
plt.plot(x,y_noisy,'r', label='noisy data', alpha=0.5)
plt.plot(x[width:-width],y_smoothed,'k', label='smoothed data')
plt.xlim(-.4,3.5)
plt.ylim(-1,2)
plt.xticks([])
plt.yticks([])
plt.legend(loc='upper left')

Out[21]:
<matplotlib.legend.Legend at 0x139f8c390>

## Linear Filter in 2D¶

• Signal: $\mathbf{f}$
• Kernel (also sometimes called a filter or a mask): $\mathbf{h} \in \mathbb{R}^{(2w+1)\times(2h+1)}$

### Cross-Correlation¶

$$CC(i,j) = \sum_{\substack{k \in [-w,w] \\ l \in [-h,h]}} \mathbf{f}(i + k, j+l) \mathbf{h}(k,l)$$

### Convolution¶

$$(\mathbf{f} \ast \mathbf{k} )_{i,j} = = \sum_{\substack{k \in [-w,w] \\ l \in [-h,h]}} \mathbf{f}(i - k, j - l) \mathbf{h}(k,l)$$

Convolution is equivalent to flipping the filter (kernel) in both directions. Convolution and cross-correlation are the same for symmetric filters.

In [22]:
f = np.linspace(0,15,16).reshape((4,4))
h = np.ones((3,3))

print(f'f = \n{f}')
print(f'h = \n{h}')

f =
[[ 0.  1.  2.  3.]
[ 4.  5.  6.  7.]
[ 8.  9. 10. 11.]
[12. 13. 14. 15.]]
h =
[[1. 1. 1.]
[1. 1. 1.]
[1. 1. 1.]]


### Exercise¶

Compute $(\mathbf{f} \ast \mathbf{k} )_{i,j}$ where $(i,j) = (1,2)$.

In [23]:
# Your answer goes here

foo = f[:3,1:]
print(foo)
print(np.reshape(foo, (9)))
print(np.reshape(h, (9)))

print(np.dot(np.reshape(foo, (9)), np.reshape(h, (9))))

[[ 1.  2.  3.]
[ 5.  6.  7.]
[ 9. 10. 11.]]
[ 1.  2.  3.  5.  6.  7.  9. 10. 11.]
[1. 1. 1. 1. 1. 1. 1. 1. 1.]
54.0


### 2D Averaging Kernels¶

A 3-by-3 summing kernel

In [24]:
average_3_by_3 = np.array([[1,1,1],
[1,1,1],
[1,1,1]], dtype='float32')


### Exercise¶

Construct a 3-by-3 averaging kernel.

In [25]:
# Your answer goes here

average_3_by_3 = average_3_by_3/np.sum(average_3_by_3)
print(f'average_3_by_3 = \n{average_3_by_3}')

average_3_by_3 =
[[0.11111111 0.11111111 0.11111111]
[0.11111111 0.11111111 0.11111111]
[0.11111111 0.11111111 0.11111111]]


### Convolution with averaging kernels¶

In [26]:
average_3_by_3 = np.ones(9).reshape(3,3) / 9
average_5_by_5 = np.ones(25).reshape(5,5) / 25
average_7_by_7 = np.ones(49).reshape(7,7) / 49

In [27]:
print(f'average_5_by_5 = \n{average_5_by_5}')

average_5_by_5 =
[[0.04 0.04 0.04 0.04 0.04]
[0.04 0.04 0.04 0.04 0.04]
[0.04 0.04 0.04 0.04 0.04]
[0.04 0.04 0.04 0.04 0.04]
[0.04 0.04 0.04 0.04 0.04]]

In [28]:
I_gray = cv.cvtColor(I, cv.COLOR_BGR2GRAY)

I3 = sp.signal.convolve2d(I_gray, average_3_by_3, mode='same')
I5 = sp.signal.convolve2d(I_gray, average_5_by_5, mode='same')
I7 = sp.signal.convolve2d(I_gray, average_7_by_7, mode='same')

plt.figure(figsize=(20,10))
plt.subplot(241)
plt.imshow(I_gray, cmap='gray')
plt.subplot(242)
plt.imshow(I3, cmap='gray')
plt.subplot(243)
plt.imshow(I5, cmap='gray')
plt.subplot(244)
plt.imshow(I7, cmap='gray')
plt.subplot(245)
plt.imshow(I_gray[215:290,175:250], cmap='gray')
plt.subplot(246)
plt.imshow(I3[215:290,175:250], cmap='gray')
plt.subplot(247)
plt.imshow(I5[215:290,175:250], cmap='gray')
plt.subplot(248)
plt.imshow(I7[215:290,175:250], cmap='gray');