空间滤波本质是对像素点的邻域处理。
选取中心点(x,y)，仅对预先定义的关于点(x,y)的邻域内的像素执行操作，令运算结果为该点处的响应，对图像中每一点重复该处理。中心点移动会产生新的邻域，每个邻域对应输入图像的一个像素。若对邻域中像素执行的计算是线性的，则该操作称为线性空间滤波，否则为非线性空间滤波。
空间线性滤波可以视作原始图像与滤波器之间进行卷积运算，表示为：
g
(
x
,
y
)
=
∑
m
∑
n
f
(
m
,
n
)
h
(
x
+
m
,
y
+
n
)
g(x,y)=sum_msum_nf(m,n)h(x+m,y+n)
g(x,y)=m∑n∑f(m,n)h(x+m,y+n)
其中，h(x,y)代表滤波器，f(x,y)为原始图像，g(x,y)为滤波后的图像。numpy.convolve()函数可以实现卷积计算
补充：
大小为m×n的滤波器，m和n一般为奇数，m=2a+1, n=2b+1，常见的m=n
对图像f(x,y)中的某一点(x,y)的滤波结果g(x,y)可表示为：
g
(
x
,
y
)
=
∑
s
=
−
a
a
∑
t
=
−
b
b
h
(
s
,
t
)
f
(
x
+
s
,
y
+
t
)
g(x,y)=sum_{s=-a}^asum_{t=-b}^bh(s,t)f(x+s,y+t)
g(x,y)=s=−a∑at=−b∑bh(s,t)f(x+s,y+t)
在图像边界上，需要把图像分别拓展为a和b个像素。
目的：使图像模糊，减少图像中的噪声，可用于预处理和噪声消除
-
均值滤波
输出像素值是对应核窗口内像素的均值
如：3*3的均值滤波器掩模为： 1 9 [ 1 1 1 1 1 1 1 1 1 ] frac{1}{9}begin{bmatrix}1 & 1 & 1 \ 1 & 1 & 1 \ 1 & 1 & 1\ end{bmatrix} 91⎣⎡111111111⎦⎤
-
中值滤波
将图像中的每个像素用领域像素中的中值来代替，对椒盐噪声有很好的抑制作用
-
高斯滤波
对应窗口为归一化后的二维高斯分布矩阵，高斯滤波器对高斯分布的噪声效果比较好。二维高斯函数如公式（9）：
g ( x , y ) = 1 2 π σ 2 e − x 2 + y 2 2 σ 2 g(x,y)=frac{1}{2pisigma^2}e^{-frac{x^2+y^2}{2sigma^2}} g(x,y)=2πσ21e−2σ2x2+y2
中值滤波python代码实现如下：
# python image median filtering
import numpy as np
import PIL.Image as Image
import sys
import matplotlib.pyplot as plt
def MedianFilter( f, KernelSize ):
m = KernelSize.shape[0]
n = KernelSize.shape[1]
a = int(m/2)
b = int(n/2)
M = f.shape[0]
N = f.shape[1]
# 补边
padded = np.zeros( (M + 2*a, N + 2*b) )
padded[ a:-a, b:-b ] = f
indices = np.zeros(b, dtype=np.int32)
padded[a:-a, :b] = f[:,indices]
indices += N-1
padded[a:-a, -b:] = f[:,indices]
padded[0:a, :]=padded[a,:]
padded[-a:,:]=padded[-a-1,:]
g = np.zeros_like( f )
# 取中值
for r in range( M ):
for c in range( N ):
g[r,c] = np.median(padded[r:r+m,c:c+n])
return g
def main( fn1, fn2):
f = Image.open(fn1)
if f.mode == 'RGB' or f.mode=='RGBA':
f = f.convert('L')
f = np.array(f)
g = MedianFilter( f, 3 )
g = Image.fromarray( g )
g.save( fn2 )
if __name__=='__main__':
if len( sys.argv ) == 3:
main(sys.argv[1], sys.argv[2])
else:
print('Usage: ', sys.argv[0], "[image file 1] [image file 2]" )
高斯滤波核构造代码如下：
def Guassian1D(x,sigma):
g = 1./(2.*np.pi*np.square(sigma))*np.exp(-np.square(x)/(2.*np.square(sigma)))
return g
def Guassian2D(x,y,sigma):
g = 1./(2.*np.pi*np.square(sigma))*np.exp(-(np.square(x)+np.square(y))/(2.*np.square(sigma)))
return g
def CalcGaussianFilterKernelSize(sigma):
x = np.linspace(-25,25,51)
y = Guassian1D(x,sigma)
threshold = Guassian1D(3.*sigma,sigma)
mask = y > threshold
y = y[mask]
KernelSize = len(y)
if not KernelSize % 2:
KernelSize += 1
return KernelSize
def CalcGaussianFilterKernel(sigma):
KernelSize = CalcGaussianFilterKernelSize(sigma)
HalfKS = KernelSize/2
g = np.linspace(-HalfKS,HalfKS,KernelSize)
x,y = np.meshgrid(g,g)
kernel = Guassian2D(x,y,sigma)
kernel = kernel/np.sum(kernel)
return kernel
B. 锐化滤波
目的：突出图像中的微小细节、增强图像中被模糊的细节
可以通过空间域的差值/差分实现图像的锐化处理。
本质上讲, 在一点进行的微分运算的响应强度正比于该点处图像灰度值的不连续的程度, 因此, 图像微分将增强图像中的边缘及其它灰度不连续(如噪声), 而削弱灰度缓慢变化的区域.
一阶采用微分算子，二阶采用拉普拉斯算子
以拉普拉斯算子为例：可以得到两种滤波器掩模： [ 0 1 0 1 − 4 1 0 1 0 ] begin{bmatrix}0 & 1 & 0 \ 1 & -4 & 1 \ 0 & 1 & 0\ end{bmatrix} ⎣⎡0101−41010⎦⎤ [ 1 1 1 1 − 8 1 1 1 1 ] begin{bmatrix}1 & 1 & 1 \ 1 & -8 & 1 \ 1 & 1 & 1\ end{bmatrix} ⎣⎡1111−81111⎦⎤
在处理后的图像中只包含了边缘和其它灰度突变区域以及较暗的背景.
要使背景保持又突出边缘，可以将原图像与经过拉普拉斯算子处理后的图像叠加
g
(
x
,
y
)
=
{
f
(
x
,
y
)
−
∇
2
f
(
x
,
y
)
掩
模
中
心
为
负
f
(
x
,
y
)
+
∇
2
f
(
x
,
y
)
掩
模
中
心
为
正
g(x,y)=begin{cases} f(x,y)-nabla^2f(x,y) & 掩模中心为负 \ f(x,y)+nabla^2f(x,y) & 掩模中心为正 end{cases}
g(x,y)={f(x,y)−∇2f(x,y)f(x,y)+∇2f(x,y)掩模中心为负掩模中心为正
此时滤波器掩模为：
[
0
1
0
1
5
1
0
1
0
]
begin{bmatrix}0 & 1 & 0 \ 1 & 5 & 1 \ 0 & 1 & 0\ end{bmatrix}
⎣⎡010151010⎦⎤
[
1
1
1
1
9
1
1
1
1
]
begin{bmatrix}1 & 1 & 1 \ 1 & 9 & 1 \ 1 & 1 & 1\ end{bmatrix}
⎣⎡111191111⎦⎤
Python中的numpy提供了fft变换，包括：一维傅里叶变换fft.fft，一维傅里叶逆变换fft.ifft，二维傅里叶变换fft.fft2，二维傅里叶逆变换fft.ifft2；以及用于DC移位的fft.fftshift以及fft.ifftshift
补充：直流分量DC
当u=0，v=0时， F ( u , v ) = 1 M N ∑ x = 0 M − 1 ∑ y = 0 N − 1 f ( x , y ) F(u,v)=frac{1}{MN}sum_{x=0}^{M-1}sum_{y=0}^{N-1}f(x,y) F(u,v)=MN1∑x=0M−1∑y=0N−1f(x,y)是函数均值，也是DC分量；
频率滤波又可分为锐化滤波和平滑滤波。通过频率域的滤波器来抑制低频分量可以实现高通滤波（锐化），抑制高频可以实现低频滤波（平滑）。
滤波步骤：
-
变换前u给图像乘上 ( − 1 ) x + y (-1)^{x+y} (−1)x+y以使变换后的DC移至频率矩形的中心位置(M/2, N/2)
-
对图像计算离散傅里叶变换F(u,v)
-
将F(u,v)与一个滤波器函数H(u,v)相乘， G(u,v)=H(u,v)F(u,v)-逐点进行
-
求上述乘积的逆傅立叶变换,只保留其中的实部
-
给上述结果再乘以 ( − 1 ) x + y (-1)^{x+y} (−1)x+y就得到滤波后的图像
说明：滤波器H一般为实函数，F为负数，滤波时H同时与F的实部和虚部相乘，这一过程不改变F的相位，又称零相移滤波器
Python实现：
def FreqFilter(f, H):
f = np.fft.fft2(f) #先进行fft变换
fshift = np.fft.fftshift(f) #频偏，将DC分量搬移到频谱中心
g = np.multiply(fshift, H) #将搬移后的频谱与滤波器H相乘，等价于空域卷积
g = np.fft.ifftshift(g) #搬移频谱
g = np.fft.ifft2(g) #逆fft变换
g = np.real(g) #保留实部
g = np.uint8(g)
return g
A. 低通滤波器
理想低通滤波器：
H
(
u
,
v
)
=
{
1
D
(
u
,
v
)
≤
D
0
0
D
(
u
,
v
)
>
D
0
H(u,v)=begin{cases} 1 & D(u,v)le D_0 \ 0 & D(u,v)gt D_0 end{cases}
H(u,v)={10D(u,v)≤D0D(u,v)>D0
其中，
D
(
u
,
v
)
=
(
x
−
M
/
2
)
2
+
(
y
−
N
/
2
)
2
D(u,v)=sqrt{(x-M/2)^2+(y-N/2)^2}
D(u,v)=(x−M/2)2+(y−N/2)2
def IdealFilter(m, n, cutoff):
x = np.arange(n)
y = np.arange(m)
cx = int(n/2)
cy = int(m/2)
x2 = x[np.newaxis,:]
y2 = y[:,np.newaxis]
x2 -= cx
y2 -= cy
sxy = x2**2 + y2**2
mask = sxy < cutoff**2 #cutoff为截止频率半径
H = np.zeros((m,n))
H[mask] = 1
return H
B. 高斯滤波器
高斯平滑滤波器：
H
(
u
,
v
)
=
e
−
D
2
(
u
,
v
)
/
2
σ
2
H(u,v)=e^{-D^2(u,v)/2sigma^2}
H(u,v)=e−D2(u,v)/2σ2
其中，
D
(
u
,
v
)
=
(
x
−
M
/
2
)
2
+
(
y
−
N
/
2
)
2
D(u,v)=sqrt{(x-M/2)^2+(y-N/2)^2}
D(u,v)=(x−M/2)2+(y−N/2)2
参数
σ
sigma
σ决定了滤波器的形状，
σ
sigma
σ越大，滤波器的平滑效果越好。
def GuassFilter(m, n, sigma):
x = np.arange(n)
y = np.arange(m)
cx = int(n/2)
cy = int(m/2)
x2 = x[np.newaxis,:]
y2 = y[:,np.newaxis]
x2 -= cx
y2 -= cy
sxy = x2**2 + y2**2
H = np.exp(-sxy/(2*sigma**2))
return H
锐化滤波与平滑滤波满足：
H
h
p
(
u
,
v
)
=
1
−
H
l
p
(
u
,
v
)
H_{hp}(u,v)=1-H_{lp}(u,v)
Hhp(u,v)=1−Hlp(u,v)
参考
