第五章图像复原与重建
1.前言
前章的图像增强主要是主观,而本章的图像复原主要是客观的,通过图像的退化现象先验知识复原已退化的图像,简单说就是用模型进行逆过程恢复图像
2.图像退化/复原处理的一个模型
本章令图像退化建模为算子 H \mathcal{H} H,这个算子与加性噪声项与原图像 f ( x , y ) f(x,y) f(x,y)运算,得到退化图像 g ( x , y ) g(x,y) g(x,y),图像复原目的就是通过上面三个得到原图像估计 f ‾ ( x , y ) \overline f(x,y) f(x,y)
空间域中退化图像如下
g
(
x
,
y
)
=
(
h
★
f
)
(
x
,
y
)
+
η
(
x
,
y
)
g(x,y)=(h\bigstar f)(x,y)+\eta(x,y)
g(x,y)=(h★f)(x,y)+η(x,y)在频率域中如下
G
(
u
,
v
)
=
H
(
u
,
v
)
F
(
u
,
v
)
+
N
(
u
,
v
)
G(u,v)=H(u,v)F(u,v)+N(u,v)
G(u,v)=H(u,v)F(u,v)+N(u,v)这两个式子是本章的重要基础
3.噪声模型
数字图像的噪声来源于图像获取、传输过程中,传感器温度、光照水平等都有影响,而频率性质指傅立叶域中噪声的频率含量,当噪声的傅立叶谱为常量时,称为白噪声
一些重要的噪声概率密度函数
下面是常见的几种噪声DFT:
- 高斯噪声 p ( z ) = 1 2 π σ e − ( z − z ‾ ) 2 / 2 σ 2 , − ∞ < z < ∞ p(z)=\frac{1}{\sqrt{2\pi}\sigma}e^{-(z-\overline z)^2/2\sigma^2},-\infty<z<\infty p(z)=2πσ1e−(z−z)2/2σ2,−∞<z<∞ z z z表灰度 z ‾ \overline z z是 z z z的平均值, σ \sigma σ是 z z z的标准差
- 瑞利噪声 p ( z ) = { 2 b ( z − a ) e − ( z − a ) 2 / b , z ≥ a 0 , z < a p(z)=\begin{cases}\frac{2}{b}(z-a)e^{-(z-a)^2/b},&z\geq a\\0,&z<a\end{cases} p(z)={b2(z−a)e−(z−a)2/b,0,z≥az<a其中 z ‾ = a + π b / 4 \overline z=a+\sqrt{\pi b/4} z=a+πb/4, σ 2 = b ( 4 − π ) 4 \sigma^2=\frac{b(4-\pi)}{4} σ2=4b(4−π)
- 爱尔兰(伽马)噪声 p ( z ) = { a b z b − 1 ( b − 1 ) ! e − a z , z ≥ 0 0 , z < 0 p(z)=\begin{cases}\frac{a^bz^{b-1}}{(b-1)!}e^{-az},&z\geq 0\\0,&z<0\end{cases} p(z)={(b−1)!abzb−1e−az,0,z≥0z<0其中 z ‾ = b a \overline z=\frac{b}{a} z=ab, σ 2 = b a 2 \sigma^2=\frac{b}{a^2} σ2=a2b
- 指数噪声 p ( z ) = { a e − a z , z ≥ 0 0 , z < 0 p(z)=\begin{cases}ae^{-az},&z\geq 0\\0,&z<0\end{cases} p(z)={ae−az,0,z≥0z<0其中 z ‾ = 1 a \overline z=\frac{1}{a} z=a1, σ 2 = 1 a 2 \sigma^2=\frac{1}{a^2} σ2=a21
- 均匀噪声 p ( z ) = { 1 ( b − a ) , a ≤ z ≤ b 0 , 其他 p(z)=\begin{cases}\frac{1}{(b-a)},&a\leq z\leq b\\0,&其他\end{cases} p(z)={(b−a)1,0,a≤z≤b其他其中 z ‾ = a + b 2 \overline z=\frac{a+b}{2} z=2a+b, σ 2 = ( b − a ) 2 12 \sigma^2=\frac{(b-a)^2}{12} σ2=12(b−a)2
- 椒盐噪声
p
(
z
)
=
{
P
s
,
z
=
2
k
−
1
P
p
,
z
=
0
1
−
(
P
s
+
P
p
)
,
z
=
V
p(z)=\begin{cases}P_s,&z=2^k-1\\P_p,&z=0\\1-(P_s+P_p),&z=V\end{cases}
p(z)=⎩
⎨
⎧Ps,Pp,1−(Ps+Pp),z=2k−1z=0z=V其中
V
V
V是
0
~
2
k
−
1
0~2^k-1
0~2k−1间的任意整数,
z
‾
=
(
0
)
P
p
+
K
(
1
−
P
S
−
P
p
)
+
(
2
k
−
1
)
P
s
\overline z=(0)P_p+K(1-P_S-P_p)+(2^k-1)P_s
z=(0)Pp+K(1−PS−Pp)+(2k−1)Ps,
σ
2
=
(
0
−
z
‾
)
2
P
p
+
(
K
−
z
‾
)
2
(
1
−
P
S
−
P
p
)
+
(
2
k
−
1
)
2
P
s
\sigma^2=(0-\overline z)^2P_p+(K-\overline z)^2(1-P_S-P_p)+(2^k-1)^2P_s
σ2=(0−z)2Pp+(K−z)2(1−PS−Pp)+(2k−1)2Ps
下面是一些噪声概率密度函数的图像
以下是一副灰度层次图进行不同噪声添加后的效果
可以看到高斯噪声使原图增加了很多的黑色小点点,亮度和对比度没有什么变化;瑞利噪声使图像深色部分的灰度级降低了,在白色区域的噪声影响比较小;爱尔兰噪声将整个图像都用白色噪声染白了;指数噪声增加了许多灰色点点,颜色比瑞利噪声要深;均匀噪声比高斯噪声产生的噪声少一些;椒盐噪声在深色产生白色点点,白色产生黑色点点,即在不同灰度产生相反色度的噪声点
周期噪声及估计噪声参数
- 周期噪声是在获取图像时由电气或机电干扰产生,空间域正弦波幅度很强导致的,可以通过频率域滤波降低,如下图所示,周期噪声对图像会产生呈现有条理周期性的噪声影响
- 噪声参数可以通过检测图像的傅立叶谱来估计,或者直接通过检测频率尖峰;由图像的一小片背景灰度可以估计DFT参数,即通过图像直方图能够计算灰度级的均值和方差,或者直方图形状分析需要得到的其他参数
4.只存在噪声的复原——空间滤波
首先介绍一下加性噪声,是指在原始图像的像素值上添加的随机干扰。加性噪声以加法的方式作用于图像,将噪声值与原始像素值相加。这些噪声会对图像质量和可视化效果产生不利影响,比如提高图像亮度、将低对比度、颜色偏差、破坏细节边缘等。
均值滤波器
只存在加性随机噪声时,就可以用空间滤波的方法对有噪声的图像 g ( x , y ) g(x,y) g(x,y)估计原来的图像 f ( x , y ) f(x,y) f(x,y),下面是一些均值滤波器的相关介绍,其中 S x y S_{xy} Sxy表示中心在 ( x , y ) (x,y) (x,y),大小为 m ∗ n m*n m∗n的领域一组坐标, g ( r , c ) g(r,c) g(r,c)表示有噪声的图像, f ^ ( x , y ) \widehat f(x,y) f (x,y)表示去噪后的图像
滤波器名称 | 复原图像 f ^ ( x , y ) \widehat f(x,y) f (x,y) | 作用效果 | 备注 |
---|---|---|---|
算术平均滤波器 | f ^ ( x , y ) = 1 m n ∑ ( r , c ) ∈ S x y g ( r , c ) \widehat f(x,y)=\frac{1}{mn}\sum_{(r,c)\in S_{xy}}g(r,c) f (x,y)=mn1∑(r,c)∈Sxyg(r,c) | 平滑局部变化,降低噪声,但会模糊图像 | |
几何均值滤波器 | f ^ ( x , y ) = [ ∏ ( r , c ) ∈ S x y g ( r , c ) ] 1 m n \widehat f(x,y)=\left[\prod_{(r,c)\in S_{xy}} g(r,c) \right]^{\frac{1}{mn}} f (x,y)=[∏(r,c)∈Sxyg(r,c)]mn1 | 实现的平滑效果可与均值滤波相比,损失图像细节更少 | |
谐波滤波器 | f ^ ( x , y ) = m n ∑ ( r , c ) ∈ S x y 1 g ( r , c ) \widehat f(x,y)=\frac{mn}{\sum_{(r,c)\in S_{xy}}\frac{1}{g(r,c)}} f (x,y)=∑(r,c)∈Sxyg(r,c)1mn | 能处理盐粒噪声和高斯类噪声,但不能处理胡椒噪声 | |
反谐波滤波器 | f ^ ( x , y ) = ∑ ( r , c ) ∈ S x y g ( x , y ) Q + 1 ∑ ( r , c ) ∈ S x y g ( x , y ) Q \widehat f(x,y)=\frac{\sum_{(r,c)\in S_{xy}}g(x,y)^{Q+1}}{\sum_{(r,c)\in S_{xy}}g(x,y)^{Q}} f (x,y)=∑(r,c)∈Sxyg(x,y)Q∑(r,c)∈Sxyg(x,y)Q+1 | 适用于降低消除椒盐噪声,Q为正时消除胡椒噪声,为负时消除盐粒噪声,Q=0时为均值滤波,Q=-1时为谐波滤波器 | Q为滤波器的阶数 |
不同的滤波器适合于不同的噪声,以算数平均滤波和几何均值滤波为例子,有如下代码操作
import cv2
import numpy as np
import matplotlib.pyplot as plt
def add_gaussian_noise(image, mean, stddev):
noise = np.random.normal(mean, stddev, image.shape)
noisy_image = image + noise
noisy_image = np.clip(noisy_image, 0, 255).astype(np.uint8)
return noisy_image
def geometric_mean_filter(image, kernel_size):
kernel = np.ones((kernel_size, kernel_size), dtype=np.float32) / (kernel_size**2)
filtered_image = cv2.filter2D(image, -1, kernel)
return filtered_image
def arithmetic_mean_filter(image, kernel_size):
kernel = np.ones((kernel_size, kernel_size), dtype=np.float32) / (kernel_size**2)
filtered_image = cv2.filter2D(image, -1, kernel)
return filtered_image
# 读取输入图像
image = cv2.imread('./img/cat1.jpeg', cv2.IMREAD_GRAYSCALE)
# 添加高斯噪声
mean = 0
stddev = 100
noisy_image = add_gaussian_noise(image, mean, stddev)
# 应用几何均值滤波器
geometric_kernel_size = 3
geometric_filtered_image = geometric_mean_filter(noisy_image, geometric_kernel_size)
# 应用算术平均滤波器
arithmetic_kernel_size = 3
arithmetic_filtered_image = arithmetic_mean_filter(noisy_image, arithmetic_kernel_size)
# 创建一个2x2的图像网格,用于显示原始图像、噪声图像、几何均值滤波器结果和算术平均滤波器结果
fig, axs = plt.subplots(2, 2)
# 显示原始图像
axs[0, 0].imshow(image, cmap='gray')
axs[0, 0].set_title('Original Image')
# 显示添加高斯噪声后的图像
axs[0, 1].imshow(noisy_image, cmap='gray')
axs[0, 1].set_title('Noisy Image')
# 显示几何均值滤波器结果
axs[1, 0].imshow(geometric_filtered_image, cmap='gray')
axs[1, 0].set_title('Geometric Filtered Image')
# 显示算术平均滤波器结果
axs[1, 1].imshow(arithmetic_filtered_image, cmap='gray')
axs[1, 1].set_title('Arithmetic Filtered Image')
# 调整子图之间的间距
plt.tight_layout()
# 显示图像
plt.show()
在代码中先是调用添加高斯噪声的函数对原图像添噪,再通过算数平均滤波和几何均值滤波进行去噪,两者都有一定的去噪效果,差别不是很大,但是几何均值比算数平均的模糊效果稍微好了一些
下面是谐波和反谐波滤波的相关实验过程,代码仅展示重要部分
def harmonic_mean_filter(image, kernel_size):
inverted_image = 1.0 / (image + 1e-6)
kernel = np.ones((kernel_size, kernel_size), dtype=np.float32) / (kernel_size**2)
filtered_image = cv2.filter2D(inverted_image, -1, kernel)
filtered_image = 1.0 / (filtered_image + 1e-6)
return filtered_image
def contra_harmonic_mean_filter(image, kernel_size, q):
numerator = np.power(image, q+1)
denominator = np.power(image, q)
kernel = np.ones((kernel_size, kernel_size), dtype=np.float32) / (kernel_size**2)
numerator_filtered = cv2.filter2D(numerator, -1, kernel)
denominator_filtered = cv2.filter2D(denominator, -1, kernel)
filtered_image = numerator_filtered / (denominator_filtered + 1e-6)
return filtered_image
# 读取输入图像
image = cv2.imread('input_image.jpg', cv2.IMREAD_GRAYSCALE)
# 添加冲激噪声
salt_prob = 0.02 # 盐噪声的概率
pepper_prob = 0.02 # 椒噪声的概率
noisy_image = add_impulse_noise(image, salt_prob, pepper_prob)
# 应用谐波滤波器
harmonic_kernel_size = 3
harmonic_filtered_image = harmonic_mean_filter(noisy_image, harmonic_kernel_size)
# 应用反谐波滤波器
contra_harmonic_kernel_size = 3
q_value = 1.5 # 这里可以调整q的值
contra_harmonic_filtered_image = contra_harmonic_mean_filter(noisy_image, contra_harmonic_kernel_size, q_value)
# 显示图像
plt.show()
同之前的实验步骤一样都是先添加噪声再进行滤波处理,在上图中可以看到,谐波和反谐波滤波对亮暗噪声的不同效果,因为图像中暗噪声更多,反谐波的去噪效果会更好一些
统计排序滤波器
滤波器名称 | 复原图像 f ^ ( x , y ) \widehat f(x,y) f (x,y) | 作用效果 | 备注 |
---|---|---|---|
中值滤波器 | f ^ ( x , y ) = m e d i a n ( r , c ) ∈ S x y { g ( r , c ) } \widehat f(x,y)=median_{(r,c)\in S_{xy}}\{g(r,c)\} f (x,y)=median(r,c)∈Sxy{g(r,c)} | 与线性平滑滤波器比,能降低某些随机噪声,且模糊度小得多 | |
最大值滤波器 | f ^ ( x , y ) = m a x ( r , c ) ∈ S x y { g ( r , c ) } \widehat f(x,y)=max_{(r,c)\in S_{xy}}\{g(r,c)\} f (x,y)=max(r,c)∈Sxy{g(r,c)} | 用于找图像最亮点,或削弱亮色相邻的暗色区域,降低胡椒噪声 | |
最小值滤波器 | f ^ ( x , y ) = m i n ( r , c ) ∈ S x y { g ( r , c ) } \widehat f(x,y)=min_{(r,c)\in S_{xy}}\{g(r,c)\} f (x,y)=min(r,c)∈Sxy{g(r,c)} | 与最大值滤波相反,降低盐粒噪声 | |
中点滤波器 | f ^ ( x , y ) = 1 2 [ m a x ( r , c ) ∈ S x y { g ( r , c ) } + m i n ( r , c ) ∈ S x y { g ( r , c ) } ] \widehat f(x,y)=\frac{1}{2}[max_{(r,c)\in S_{xy}}\{g(r,c)\}+min_{(r,c)\in S_{xy}}\{g(r,c)\}] f (x,y)=21[max(r,c)∈Sxy{g(r,c)}+min(r,c)∈Sxy{g(r,c)}] | 适合于随机分布的噪声,如高斯噪声货均匀噪声 | |
修正 α \alpha α均值滤波器 | f ^ ( x , y ) = 1 m n − d ∑ ( r , c ) ∈ S x y { g R ( r , c ) } \widehat f(x,y)=\frac{1}{mn-d}\sum_{(r,c)\in S_{xy}}\{g_R(r,c)\} f (x,y)=mn−d1∑(r,c)∈Sxy{gR(r,c)} | d=0时为均值滤波器,d=mn-1时为中值滤波器,其他时适合处理混合噪声,如高斯、椒盐噪声 | 设在 S x y S_{xy} Sxy删除 g ( r , c ) g(r,c) g(r,c)的 d / 2 d/2 d/2个最低和最高灰度值 |
以最大和最小值滤波为例,进行实验,仅展示相关代码
#定义两个滤波函数
def max_filter(image, kernel_size):
kernel = np.ones((kernel_size, kernel_size), dtype=np.uint8)
filtered_image = cv2.dilate(image, kernel)
return filtered_image
def min_filter(image, kernel_size):
kernel = np.ones((kernel_size, kernel_size), dtype=np.uint8)
filtered_image = cv2.erode(image, kernel)
return filtered_image
# 读取输入图像
image = cv2.imread('./img/cat1.jpeg', cv2.IMREAD_GRAYSCALE)
# 应用最大值滤波器
max_kernel_size = 3
max_filtered_image = max_filter(noisy_image, max_kernel_size)
# 应用最小值滤波器
min_kernel_size = 3
min_filtered_image = min_filter(noisy_image, min_kernel_size)
# 显示图像
plt.show()
实验步骤与之前的类似,主要分析滤波的效果。可以看到最大值和最小值滤波对亮暗噪声的响应效果,都是把噪声点弄成了最大或者最小的灰度级,导致产生了黑白点,在最小值滤波中,对于猫深色区域的效果会好,而最大值滤波则有相反的效果
自适应滤波器
自适应滤波器顾名思义,就是适应图像中不同点的特征变化的滤波器,能够根据输入信号的统计特性自动调整滤波器参数,因此滤波性能很好,但相应的滤波器复杂度加大,有下面两种自适应滤波器
- 自适应局部降噪滤波器
以中心为 ( x , y ) (x,y) (x,y)的领域 S x y S_{xy} Sxy进行操作,令噪声的方差为 σ η 2 \sigma_{\eta }^2 ση2,局部平均灰度为 z ‾ S x y \overline z_{S_{xy}} zSxy,局部方差为 σ S x y 2 \sigma_{S_{xy}}^2 σSxy2,则表达式为 f ^ ( x , y ) = g ( x , y ) − σ η 2 σ S x y 2 [ g ( x , y ) − z ‾ S x y ] \widehat f(x,y)=g(x,y)-\frac{\sigma_{\eta}^2}{\sigma_{S_{xy}}^2} [g(x,y)-\overline z_{S_{xy}}] f (x,y)=g(x,y)−σSxy2ση2[g(x,y)−zSxy] - 自适应中值滤波器
自适应滤波器优点是能处理更大概率的噪声,且保留图像细节同时平滑非冲激噪声,其基本原理如下:- 定义一个滑动窗口,包含待滤波的像素及其邻域像素。
- 将窗口内的像素按照灰度值进行排序,找到中间值作为中值。
- 计算中值与当前像素之间的差值(称为差值度量)。
- 如果差值度量小于预设的阈值,认为当前像素不是噪声,将其保留。
- 如果差值度量大于等于阈值,认为当前像素受到噪声影响,用中值替代当前像素。
- 滑动窗口移动到下一个像素位置,重复步骤2~5。
- 对整个图像重复以上步骤,直到所有像素都被处理。
可以看到,与前面几个滤波相比,经过自适应中值滤波器处理的图像去噪效果更加明显,
5.使用频率域滤波降低周期噪声
陷波滤波深入介绍
陷波滤波器的传递函数一般形式为
H
N
R
(
u
,
v
)
=
∏
k
=
1
Q
H
k
(
u
,
v
)
H
−
k
(
u
,
v
)
H_{NR}(u,v)=\prod_{k=1}^{Q}H_k(u,v)H_{-k}(u,v)
HNR(u,v)=k=1∏QHk(u,v)H−k(u,v)其中三种常见的滤波器之前已经有了解,下面是理想、巴特沃斯和高斯陷波带阻滤波器传递函数的三维图(透视图)
陷波滤波器能降低图像的噪声消除干扰
最优陷波滤波
最优滤波器比一般的陷波滤波更能适应多个干扰分量的存在产生的影响,并能检测有干扰信息的宽裙摆,主要过程如下:
- 确定需要抑制的干扰频率或频率范围。
- 计算干扰频率在频率域上的位置。
- 根据干扰频率位置设计一个陷波滤波器,使其在干扰频率处有较高的衰减,并在其他频率处保持较低的衰减。
- 将设计好的陷波滤波器应用于信号,以抑制干扰频率的干扰成分。
- 得到经过最优陷波滤波后的信号,其中干扰频率的成分被显著减弱或消除。
6.估计退化函数
估计图像复原的退化函数主要有三种:(1)观察法(2)试验法(3)数学建模法,这三种过程称为盲去卷积,强调真正的退化函数很少完全已知,以下的公式都基于基于基本公式 G ( u , v ) = H ( u , v ) F ( u , v ) + N ( u , v ) G(u,v)=H(u,v)F(u,v)+N(u,v) G(u,v)=H(u,v)F(u,v)+N(u,v)
- 采用观察法估计退化函数
一般是寻找图像中信号很强的区域(比如高对比度区域)再处理图像,令 g s ( x , y ) g_s(x,y) gs(x,y)表示观察子图像 f ^ s ( x , y ) \widehat f_s(x,y) f s(x,y)表示处理后的字图像,则有公式 H s ( u , v ) = G s ( u , v ) F ^ s ( u , v ) H_s(u,v)=\frac{G_s(u,v)}{\widehat F_s(u,v)} Hs(u,v)=F s(u,v)Gs(u,v)根据图像位置不变的假设来推断完整的退化函数,在过程中可以通过锐化降噪等操作辅助处理 - 采用试验法估计退化函数
即用相同技术的已知原图来推测退化函数,例如冲激图像中通过一步步贴近要复原的图像位置,再把过程中使用的系统设置技术施加在一个冲激成像上,得到退化的冲激响应,公式如下 H ( u , v ) = G ( u , v ) A H(u,v)=\frac{G(u,v)}{A} H(u,v)=AG(u,v) - 采用建模法估计退化函数
一种是直接建立退化模型,如湍流建模;另一种是根据基本原理推导一个数学模型,例如运动模糊图像的建模
7.图像复原的其他方法:
逆滤波
逆滤波是最简单的复原方法,即用退化函数的傅立叶变换除以退化函数,来估计原图像,但往往知道退化函数也不能准确复原图像,且性能较差
最小均方误差(维纳)滤波
通常在退化函数是零或者非常小的值时直接使用逆滤波效果并不好,而维纳滤波器以使未污染图像和估计图像的均方误差最小为基础,复原效果较好,在平均意义上是最优的,但要求未退化图像和噪声的功率谱已知
以下是对一张图片进行模糊处理后再用维纳滤波复原的代码实例
import cv2
import numpy as np
import matplotlib.pyplot as plt
# 读取图像
im = cv2.imread("./img/cat1.jpeg", cv2.IMREAD_GRAYSCALE)
# 调整图像尺寸
M = max(im.shape)
im_resized = cv2.resize(im, (M, M))
# 显示原图
plt.subplot(231)
plt.imshow(im, cmap='gray')
plt.title('Original image')
# 傅里叶变换
F = np.fft.fftshift(np.fft.fft2(im_resized))
# 计算退化模型
[M, _] = im_resized.shape
T = 1
a = 0.1
b = 0.1
v = np.arange(-M / 2, M / 2)
u = v[:, np.newaxis]
A = np.tile(a * u, (1, M)) + np.tile(b * v, (M, 1))
H = T / np.pi / A * np.sin(np.pi * A) * np.exp(-1j * np.pi * A)
H[A == 0] = T
# 计算退化后的图像
SDF = F * H
im_smb = np.real(np.fft.ifft2(np.fft.ifftshift(SDF)))
plt.subplot(232)
plt.imshow(im_smb, cmap='gray')
plt.title('motion blur.')
# 加噪
noise = np.random.normal(loc=0, scale=0.00001, size=(M, M))
FNoise = np.fft.fftshift(np.fft.fft2(noise))
SDFN = F * H + FNoise
im_smbn = np.real(np.fft.ifft2(np.fft.ifftshift(SDFN)))
plt.subplot(233)
plt.imshow(im_smbn, cmap='gray')
plt.title('motion blur + Gaussian noise.')
# 无噪直接逆滤波
F_N = SDF / H
im_N = np.real(np.fft.ifft2(np.fft.ifftshift(F_N)))
plt.subplot(234)
plt.imshow(im_N, cmap='gray')
plt.title('Direct inverse filtering with no noise')
# 有噪直接逆滤波
F_NN = SDFN / H
im_NN = np.real(np.fft.ifft2(np.fft.ifftshift(F_NN)))
plt.subplot(235)
plt.imshow(im_NN, cmap='gray')
plt.title('Direct inverse filtering with noise')
# 截止半径R内逆滤波
R = 5
FD = np.zeros_like(F)
for i in range(M):
for j in range(M):
if np.sqrt((i - M/2) ** 2 + (j - M/2) ** 2) < R:
FD[i, j] = SDFN[i, j] / H[i, j]
im_NNR = np.real(np.fft.ifft2(np.fft.ifftshift(FD)))
plt.subplot(236)
plt.imshow(im_NNR, cmap='gray')
plt.title('Inverse filtering with noise and cutoff radius of 5')
# 维纳滤波
plt.figure()
plt.subplot(221)
plt.imshow(im, cmap='gray')
plt.title('Original image')
plt.subplot(222)
plt.imshow(im_smbn, cmap='gray')
plt.title('motion blur with Gaussian noise.')
# 维纳滤波公式
buf = np.abs(H) ** 2
SNR = (np.abs(FNoise) ** 2) / (np.abs(F) ** 2)
F_WN = SDFN / H * buf / (buf + SNR)
im_wf = np.real(np.fft.ifft2(np.fft.ifftshift(F_WN)))
plt.subplot(223)
plt.imshow(im_wf, cmap='gray')
plt.title('Theoretical Wiener filtering')
# k = 0.02的维纳滤波
bestK = 0.02
F_WN = SDFN / H * buf / (buf + bestK)
im_wf = np.real(np.fft.ifft2(np.fft.ifftshift(F_WN)))
plt.subplot(224)
plt.imshow(im_wf, cmap='gray')
plt.title('k=0.02 Wiener filtering.')
plt.show()
从图上可以看出维纳滤波对于模糊图像的复原效果还是比较有作用的,但是与理想化的维纳滤波相比还是有一定的差距,可以通过k值来调整,但是经过其他参数的测试,上图中的k值还是比较好的效果,代码中并没有显示
约束最小二乘方滤波
约束最小二乘方的显著特点是仅要求噪声的方差和均值,这些参数可通过给定的退化图像计算出来,并且它对每幅图像都能产生理论上的最优结果。约束最小二乘方的核心是退化函数对噪声敏感度的问题,因此需要使用问题的参数约束复原
8.由投影重建图像
投影重建图像是一种通过从不同角度或方向上获取投影数据,并利用这些投影数据重建原始图像的技术。这种技术常用于医学影像学领域,如计算机断层扫描(CT扫描)和正电子发射断层扫描(PET扫描)。
理解投影重建图像的概念可以通过以下步骤来解释:
- 获取投影数据:在投影重建过程中,通过向物体或人体传递一束射线(例如X射线或正电子)并在另一侧接收射线的强度信息,从而获取到一系列投影数据。这些投影数据是通过在不同位置或角度上对射线进行测量得到的。
- 重建过程:在投影重建图像中,通过将投影数据反投影回原始图像空间来重建图像。简单来说,就是将投影数据的信息根据它们的位置和强度值重新映射回图像空间。
- 反投影和插值:反投影过程涉及将投影数据的强度值根据其在图像中的位置进行插值,并将其放置在相应的位置。这样,通过将所有投影数据反投影回图像空间,并将它们相加,就可以逐渐重建出原始图像。
- 过滤和重建:在反投影过程中,为了提高图像质量和减少伪影,常常使用合适的滤波器对反投影数据进行处理。经过滤波后,将得到的数据进行累加,得到最终的重建图像。
傅立叶切片定理
傅立叶切片定理描述了二维图像和其在频域中的傅立叶变换之间的关系,其表达式如下 G ( ω , θ ) = [ F ( u , v ) ] u = ω c o s θ ; v = ω s i n θ = F ( ω c o s θ , ω s i n θ ) G(\omega,\theta)=[F(u,v)]_{u=\omega \ cos\theta;v=\omega \ sin\theta}=F(\omega \ cos\theta,\omega \ sin\theta) G(ω,θ)=[F(u,v)]u=ω cosθ;v=ω sinθ=F(ω cosθ,ω sinθ)
如果一个二维函数在空域中进行了平移,那么它的傅立叶变换的切片(沿着一条平行于平移方向的直线截取傅立叶变换的值)等于平移后的函数的傅立叶变换。简单来说,傅立叶切片定理说明了在频域中,二维函数的平移操作等价于其傅立叶变换的切片。
根据傅立叶切片定理,我们可以采用以下步骤进行图像重建:
- 获取投影数据:通过从不同角度或方向上获取物体的投影数据,例如通过X射线扫描或正电子发射断层扫描(PET扫描)。
- 执行傅立叶变换:对每个投影数据进行傅立叶变换,将其转换到频域。
- 平移操作:根据扫描过程中的角度和位置信息,将频域数据按照平移操作进行调整。
- 沿平行线提取切片:沿着一条平行于平移方向的直线,截取平移后的傅立叶变换数据的值,形成一个一维频谱切片。
- 反傅立叶变换:对频谱切片进行反傅立叶变换,将其转换回空域。
- 重建图像:重建出的图像即为原始物体的二维表示。
傅立叶切片定理为图像重建提供了一种有效的方法,它允许我们通过在频率域上的切片来获取信号在时域上的部分信息,并使得可以仅通过一维的频谱切片数据进行图像的重建,从而减少了计算量和存储需求。
滤波反投影重建图像
直接使用反投影会导致图像的模糊,在投影前进行一些滤波操作能降低模糊效果,比如平行射线束滤波和扇形射线束滤波,两者的基本步骤都很类似,在上述投影重建图像的过程中已经显示了,在区别上,扇形射线束滤波更接近于实际的CT扫描设备,具有更广的角度覆盖范围,一般使用的范围广;平行射线束滤波则是简化的版本,适用于一些特定场景
9.章总结
本章中对于一些基本噪声都有了了解,对于它们的特点以及合适与不同噪声的线性、非线性、自适应空间滤波都有了对应的关系,例如几何均值滤波适合去除高斯噪声,自适应中值滤波适合于椒盐噪声,陷波滤波适合于周期性滤波等等,正逆滤波是一个相反的过程,也是根据条件的合适度使用,但逆滤波恢复方法对噪声极为敏感,通常难以直接使用得到较好的结果。维纳滤波可在有噪声条件下,从退化图像复原出去噪图像的估计值,使得该估计值符合一定的准则。反投影也是图像复原的一个技术,通过平行射线束滤波和扇形射线束滤波都可以得到估计的原图像,扇形射线束滤波一般更好