基础
- 复原的目的:针对质量降低或失真的图像,试图恢复其原始的内容或质量
- 基本要点
- 图像复原总是试图寻找引起图像质量下降的客观原因,有针对性地进行“复原”处理
- 获得使图像质量下降的先验知识,建立退化模型是图像复原处理的前提与关键
- 图像恢复总是假定已知或可以通过估计得到引起图像降质的模型,而图像增强不需要知道图像降质模型
- 引起图像质量下降的客观因素
- 成象系统的象差、畸变、带宽有限等造成图像失真
- 几何失真: 由于成象器件拍摄姿态和扫描非线性引起
- 灰度失真: 光学系统或成象传感器本身特性不均匀,造成同样亮度景物成象灰度不同
- 运动模糊: 成象传感器与被拍摄景物之间的相对运动,引起所成图像的运动模糊
- 辐射失真: 由于场景能量传输通道中的介质特性如大气湍流效应、大气成分变化引起图像失真
- 噪声干扰: 图像在成象、数字化、采集和处理过程中引入的噪声等
- 针对上述客观因素,建立图像退化模型,是进行图像复原处理的第一步
- 图像退化与复原模型
- 降质图像可表示为:$g(x,y) = H[f(x,y)]+η(x, y)$
- $H[.]$ 是综合所有退化因素的函数
- $\eta(x,y)$ 是加性噪声
- 假定噪声干扰为零,进一步考察上述模型
- 若 $H[af_1(x,y) + bf_2 (x,y)] = aH[f_1(x,y)]+ bH[f_2(x,y)]$, 则该系统为线性系统
- 若 $H[f(x-\alpha,y-\beta)] = g(x-\alpha,y-\beta)$, 则该系统为位移不变系统
- 线性移不变退化模型
- 连续函数退化模型:$g(x,y)=f(x,y)*h(x,y)+\eta(x,y)$, 及频域表示:$G(u,v) = H(u,v)F(u,v)+N(u,v)$
- 降质图像可表示为:$g(x,y) = H[f(x,y)]+η(x, y)$
- 基本原理
- 根据退化原因,建立相应的数学模型,从被污染或畸变的图像信号中提取所需要的信息——沿着使图像降质的逆过程恢复图像本来面貌
- 实际的复原过程相当于设计一个滤波器,使其能从降质图像中计算得到真实图像的估值 ,使其根据预先规定的误差准则,最大程度地接近真实图像。
- 广义上讲,图像复原是一个求逆(反演)问题,逆问题经常存在非唯一解,甚至无解。为了得到逆问题的有用解,需要有先验知识以及对解的附加约束条件
- 引起图像退化的原因常常为非线性的,非线性以及与位置相关的退化处理技术,虽然可得到更加精确的结果,但在处理上将会非常困难并常常可能无解
- 对许多退化过程用线性系统近似不但可以方便求解,而且可得到基本满意的结果
- 从线性系统角度,图像的退化可看作为原始图像与退化函数的卷积,因此线性图像复原往往称之为“图像去卷积”,所采用的滤波器称之为“去卷积滤波器”
噪声单独干扰下的图像滤波复原
噪声模型
- 类型:
- 白噪声(White noise)
- 白噪声的傅里叶频谱为常数
- 假定白噪声与空间坐标系相互独立
- 假定白噪声与图像像素之间相互独立
- 周期性噪声
- 噪声分布与空间坐标系相关
- 大多数周期性噪声可通过频域滤波基本消除
- 脉冲噪声
- 在图像中引起黑、白点状的随机噪声
- 白噪声(White noise)
- 噪声模型
- 高斯噪声(Gaussian noise):$p(z)=\frac{1}{\sqrt {2\pi} \sigma} e^{-\frac{(z-\mu)^2}{2\sigma^2}}$
- 瑞利噪声(Rayleigh noise): $p(z)=\left\lbrace \begin{matrix} \frac ab(z-a)e^{-\frac{(z-a)^2}{b}} & z\ge a \\ 0 & z<a \end{matrix}\right.$, $\quad \mu=a+\sqrt {\frac{\pi b}{4}}, \quad \sigma^2=\frac{b(4-\pi)}{4}$
- 伽马噪声(Gamma noise): $p(z)=\left\lbrace \begin{matrix} \frac{a^bz^{b-1}}{(b-1)!}e^{-az} & z\ge 0 \\ 0 & z<0 \end{matrix}\right.$, $\quad \mu=\frac{b}{a}, \quad \sigma^2=\frac{b}{a^2}$
- 指数分布噪声(Exponential noise): $p(z)=\left\lbrace \begin{matrix} ae^{-az} & z\ge 0 \\ 0 & z<0 \end{matrix}\right.$, $\quad \mu=\frac{1}{a}, \quad \sigma^2=\frac{1}{a^2}$, 是伽马噪声的特殊形式(b=1)
- 均匀分布噪声(Uniform noise): $p(z)=\left\lbrace \begin{matrix} \frac{1}{b-a} & a\le z\le b \\ 0 & otherwise \end{matrix}\right.$, $\quad \mu=\frac{a+b}{2}, \quad \sigma^2=\frac{(b-a)^2}{12}$
- 脉冲噪声(椒盐噪声)(Impulse(Salt-and-Pepper)noise): $p(z)=\left\lbrace \begin{matrix} P_a & z=a \\ P_b & z=b \\ 0 & otherwise \end{matrix}\right.$
- 噪声参数的估计
- 对于Gaussian, uniform, Rayleigh, Erlang, exponential noises等噪声,关键是得到其均值与方差
- 通过传感器特性进行估计: 分析传感器成像器件特性, 使传感器针对均匀灰度图像成像
- 从图像本身进行估计: 在图像中截取具有恒定灰度值区域进行计算
图像滤波复原
- 原理:
- 引起图像降质的唯一原因是噪声干扰时,退化模型变为:$g(x,y)=f(x,y)+\eta(x,y)$
- 噪声函数$η(x,y)$一般很难得到
- 仅当噪声为周期性噪声时,可通过分析G(u,v)得到N(u,v),但难于建立分析规律
- 其他加性噪声存在时,通常选择空间滤波技术,与图像增强中滤波机理完全一样
- 结合噪声模型可建立一些特殊滤波器,其效果可优于在图像增强中常用方法
- 空间滤波器设计
- 算术均值滤波器
- $\hat f(x,y)=\frac{1}{mn} \sum_{(x,y)\in S_{xy}}g(s,t)$
- 几何均值滤波器
- $\hat f(x,y) = [\prod_{(x,y)\in S_{xy}}g(s,t)]^{\frac 1{mn}}$
- 平滑度与算术滤波器相当, 图像细节丢失更少
- 谐波均值滤波器
- $\hat f(x,y) = \frac{mn}{\sum_{(x,y)\in S_{xy}} \frac 1{g(s,t)}}$
- 对“盐”噪声较好
- 不适用于“胡椒”噪声
- 善于处理高斯噪声
- 顺序统计滤波器
- 中值滤波器, 最大值滤波器, 最小值滤波器, 中点滤波器, 修正阿尔法均值滤波器
- 算术均值滤波器
- 空间自适应滤波器
- 结合图像的局部特性自动修改滤波器参数或滤波策略
- 良好设计的自适应滤波器其性能要优于所有固定滤波器,而滤波器结构可能更为简单
- 局部噪声自适应滤波器
- 以局部区域随机变量的统计特性(均值、方差)为基础
- 设$g(x,y)=f(x,y)+\eta(x,y)$ 具有整体噪声方差:$\sigma^2_{\eta}$
- 在局部区域$S_{xy}$(为$g(x,y)$的子集),像素点具有局部均值$m_l$和方差$\sigma_l^2$
- 可有自适应滤波器 $\hat h(x,y)=g(x,y)-\frac{\sigma^2_\eta}{\sigma_l^2}[g(x,y)-m_l]$
- 该自适应滤波器基于如下的设想
- 如果噪声方差为零,则表明没有噪声; f(x,y)等于g(x,y)
- 如果局部方差高于噪声方差,表明该局部图像信息多于噪声,可能存在较多的边沿,应尽可能保留;则返回一个g(x,y)的近似值
- 如果局部方差与噪声方差相等,表明该区域噪声干扰严重;返回局部区域像素的算术平均值
- 自适应中值滤波器
- A 层:$A_1=z_{med}-z_{min}, \; A_2=z_{med}-z_{max}$,
- 如果 $A1>0$ 且 $A_2<0$ ,转到B层; 否则:增大窗口尺寸; 如果窗口尺寸 $\le S_{max}$ ,重复A层, 否则输出B层
- B层:$B_1=z_{xy}-z_{min}, \; B_2=z_{xy}-z_{max}$
- 如果 $B1>0$ 且 $B_2<0$, 输出$z_{xy}$; 否则:输出$z_{med}$
- 目标:滤出椒盐噪声(冲激噪声);平滑其他非冲激噪声;减少物体边界失真
- 设想:$z_{max}$、 $z_{min}$可粗略认为相当于冲激噪声的噪声成分
- 算法思路:
1)首先检测中值$z_{med}$是否为噪声脉冲
2)若不是,进一步判断中心像素$z_{xy}$是否为噪声脉冲,若不是,直接输出该点;若是,输出中值
3)若中值$z_{med}$为噪声脉冲,则增大窗口尺寸,直至在允许的范围内找到非噪声脉冲中值;否则输出中心像素$z_{xy}$值
- 周期性噪声频域滤波
- 带阻滤波器
- 理想带阻滤波器
- 巴特沃斯带阻滤波器
- 高斯带阻滤波器
- 陷波滤波器(Notch Filters)
- 在指定的中心频率点,阻止该中心频率点领域内的频率
- 由于傅里叶变换的对称性,陷波器必须以关于原点对称的形式出现,常常为一对或两对出现
- 理想陷波滤波器—中心在$(u_0,v_0)$并在$(-u_0,-v_0)$对称的理想陷波器
- 巴特沃斯陷波滤波器
- 高斯陷波滤波器
- 带阻滤波器
系统退化复原
- 假定退化图像中噪声干扰为零,退化模型变为:$g(x,y)=h(x,y)f(x,y)+\eta(x,y) \rightarrow g(x,y)=h(x,y)f(x,y)$
- 退化函数$H(u,v)$的估计,是进行系统退化复原的关键一步,对$H$的估计过程常称为“系统辨识过程”
- 常用方法为:观察法, 实验法, 数学建模法
线性系统退化函数的估计
- 图像观察估计法
- 通过图像自身结构信息进行估计
- 选择图像中具有强信号与强特征的局部区域图像$g_s(x,y)$,设法构建一个具有相同大小与特征、但没有退化的近似图像 $\hat f (x,y)$, 可有 $H_s(u,v)=\frac{G_s(u,v)}{\hat F_s(u,v)}$
- 利用从这一函数出发,进一步假设$H(u,v)$的数学表达形式,从而构建$H(u,v)$使其与$H_s(u,v)$具有基本相同的形状
- 试验估计法
- 使用或设计一个与图像退化过程相似的装置(过程),使其成像一个脉冲,可得到退化系统的冲激响应 $H(u,v)=\frac{G(u,v)}{A}$, $G(u,v)$为观察图像的傅里叶变换
- 模型估计法
- 从引起图像退化的基本原理进行推导,进而对原始图像进行模拟,在模拟过程中调整模型参数以获得尽可能精确的退化模型
- 运动模糊模型 — 由于物体向一个方向线性移动:
- $H(u,v)=\frac{T}{\pi (ua+ub)}\sin [\pi (ua+ub)]e^{-j(ua+ub)}$
- 运动模糊 — 照相机与景物相对运动
- 设$T$为快门时间,$x_0(t)$,$y_0(t)$是位移的$x$分量和$y$分量 $H(u,v)=\int^T_0 exp{-j2\pi(ux_0(t)+vy_0(t))}dt$
- 光学散焦模型: $H(u,v)=J_1(\pi d\rho)/\pi d \rho, \; \rho = (u^2+v^2)^{1/2}$
- $d$ 是散焦点扩展函数的直径, $J_1(•)$ 是第一类贝塞尔函数
- 点扩展函数法
- 相当一类退化函数可看作为对图像施加了某种类型的模糊操作,这种模糊操作中的退化函数常称之为点扩展函数(point spread function, PSF)
- 获得点扩展函数的方法
- 点光源法: 保持成像时的条件不变,将一个点光源输入成像系统,输出即为相应的PSF
- 频率响应法: 保持成像时的条件不变,将不同频率的光源输入成像系统,输出系列形成相应的PSF的傅里叶变换函数$H(u,v)$
- 图像像素测量法: 测量模糊图像中单个像素点发生变化的方向与强度,形成PSF模板,该模板作为该图像的退化函数。单个像素点发生变化的方向可能为:线性单一方向(水平、垂直、对角),所有方向
- 对于PSF模板,进一步的按照待处理图像的尺寸进行补零处理,使其具有相同的大小
- 对补零处理后的模板,进行循环移位,使得原始模板的中心位于图像的原点(0,0)
- 由此形成的退化函数,可采用傅里叶变换进行复原处理
逆滤波复原方法
- 若已知退化函数$H(u,v)$,则复原过程为 $F(u,v)=G(u,v)/H(u,v) \rightarrow f(x,y)$
- 当退化图像的噪声较小,即轻度降质时,采用逆滤波复原的方法可以获得较好的结果
- 逆滤波复原的病态问题
- 在$H(u,v)$等于零或非常小的数值点上,$F(u,v)$ 将变成无穷大或非常大的数
- 进一步考虑噪声影响 $F(u,v)=\frac{G(u,v)}{H(u,v)} - \frac{N(u,v)}{H(u,v)}$
- 由于噪声分布在很宽的频率空间,即使数值很小也会因为$H(u,v)$使得上式右侧第二项变得很大,噪声影响大大增强
- 改进方法: 伪逆滤波复原, 即对$H(u,v)$规定一个门限值: $H^{-1}(u,v)=\left\lbrace \begin{matrix} \frac1{H(u,v)} & |H(u,v)|>\delta \\ 0 & |H(u,v)|\le\delta \end{matrix} \right.$
- 伪逆滤波复原
- 施加圆形范围限制:$H^{-1}(u,v)=\left\lbrace \begin{matrix} \frac1{H(u,v)} & \sqrt{u^2+v^2}\le R \\ 0 & \sqrt{u^2+v^2} > R \end{matrix} \right.$
最小均方误差滤波(维纳滤波)
- 假定图像 $f(x,y)$ 和噪声 $η$ 均为随机信号,且相互之间互不相关
- 目标:寻求最佳复原图像,使得均方误差 $e^2=E{|f-\hat f|^2}$ 为最小
- 维纳滤波结果为: $\hat F(u,v)=\frac{|H(u,v)|^2}{|H(u,v)|^2+S_\eta (u,v)/S_f(u,v)}\cdot \frac{G(u,v)}{H(u,v)}$
- $H_{win}(u,v)=\frac{H^*(u,v)S_f(u,v)}{|H(u,v)|^2S_f(u,v)+S_\eta (u,v)}$
- $H(u,v)$ 为退化函数,$H^*(u,v)$ 为退化函数复共轭,$S_\eta(u,v)=E(|N(u,v)|^2)$ 为噪声功率谱,$S_f(u,v)=E(|F(u,v)|^2)$ 为未退化图像功率谱。
- 讨论:
- 维纳滤波为均方误差最小意义下的最佳滤波,可使具有噪声干扰图像的客观复原性能达到最佳
- 维纳滤波器假定图像与噪声均为平稳随机过程,若图像与噪声实际随机特性平稳性差距较大时,维纳滤波器难以得到最佳结果
- 维纳滤波器建立在最小化统计准则基础上,只是在平均意义上最优
- 维纳滤波需大量先验知识,实用中常常难以进行
- 尤其是噪声功率谱、未退化图像的功率谱很难获取
- 通过在滤波过程中调节K值以得到准最佳结果:$K=\frac{S_\eta(u,v)}{S_f(u,v)}= \frac{noise_power}{image_power}$
- 进一步假设噪声为零,K=0,维纳滤波器退化为逆滤波器
约束最小二乘方滤波器
- 问题的提出
- 已知系统退化函数H的条件下,系统退化复原的逆滤波函数理论上为1/H,但存在病态解和对噪声的高度敏感
- 维纳滤波实现了噪声干扰条件下均方意义上的最佳滤波,但存在需要大量先验知识尤其是噪声功率谱和未退化图像的功率谱的困难
- 进一步考察退化模型的矩阵形式 $g(x,y)=h(x,y)*f(x,y)+\eta (x,y) \Rightarrow g=Hf+\eta$, 应有 $g-Hf=\eta$
- 为了方便求解$f$,对上式采用泛函中范数表示 $g-HF=\eta \Rightarrow |g-Hf|^2=|\eta|^2$, 可有 $[g-Hf]^T[g-Hf]=e^2$
- 由于 $H$ 的问题,满足上式的解 $f$ 可能很多;对此可施一定的约束条件,以便从众多的 $f$ 中选择最佳的 $f$ 作为面临问题的最佳解——约束最小二乘方滤波器
- 正则化方法
- 考察无噪声条件下的复原处理模型,并表示为矩阵形式 $g(x,y)=h(x,y)*f(x,y) \Rightarrow g=Hf$
- 在没有退化系统 $H$ 先验知识的情况下,通过模型估计方法由 $g$ 求 $f$ ,原则上可归结为一种病态的反演,面临反演的不稳定性和多解性
- 对于这类不确定性问题的解决,正则化理论以某种途径提供了有效的方法,其原理是:定义一个稳定算子$q$,及$qf$的模
- 通过$q$ 的引入,$qf$ 描述了$f$ 的平滑程度;上式则描述了$f$ 在全部定义区域内的平滑程度,其值越小,$f$ 越平滑
- 在正则化方法下求解$f$ 一般有三种方法:
- 在约束$|qf|<C$的条件下求解$f$ ,使$|g-Hf|$达到最小
- 在约束$|g-Hf|<C$的条件下,求$f$使$|qf|$最小
- 求解$f$使下式达到最小: $E=|g-Hf|+\lambda |qf|$
- 根据退化模型可得 $|g-H\hat f|^2-|\eta|^2=0$
- 考虑对$\hat f$ 施加某种线性运算 $Q$,按照正则化方法,问题可归结为最小化如下目标函数 $W(\hat f)=|Q\hat f|^2+\lambda |g-H\hat f|^2-|\eta|^2$, $\lambda$ 为拉格朗日因子
- 令$W(\hat f)$ 对$\hat f$的导数为零, 可得 $\hat f=(H^TH+\gamma Q^TQ)^{-1}H^Tg$, $\gamma=1/\lambda$, 作为使等式 $|g-H\hat f|^2=|\eta|^2$ 成立的常数
- 上式称为约束最小二乘复原解的通用方程式
- 在约束最小二乘复原解通用方程式基础上,通过选择不同的$Q$,可得到不同类型的复原滤波器
- 伪逆滤波器:
- 令$Q=I$,$I$ 表示$MN×MN$单位矩阵,可有 $hat f=(H^TH+\gamma I)^{-1}H^Tg$
- 称之为伪逆滤波器
- 令$γ=0$,则上式成为标准的逆滤波器
- 参数化维纳滤波器
- 将$f、η$视为随机矢量,选择$Q$为 $\frac{\sqrt {R_\eta}}{\sqrt {R_f}}$
- 式中 $R_f=E(ff^T)$ 和 $R_\eta = E(\eta\eta^T)$ 分别为信号和噪声的协方差矩阵,可有 $\hat f=(H^TH+\gamma R_f^{-1}R_\eta)^{-1}H^Tg$
- 上式称之为参数化维纳滤波器;$γ$为可调节参数, 为 1 时上式成为经典维纳滤波器
- 平滑约束滤波器
- 对一幅带噪声模糊图像的逆滤波复原处理往往会放大噪声,产生预料之外的结果;
- 解决方法之一是选择适当的Q,对复原图像$\hat f$施加一定程度的光滑性约束;
- 由此可从约束最小二乘复原解通用方程式中导出一个光滑的、去模糊的、无噪声的复原滤波器
- 令Q对应于一个高通滤波运算(二阶拉普拉斯算子): $\nabla^2f(x,y)=[\frac{\partial^2}{\partial x^2} + \frac{\partial^2}{\partial y^2}]f(x,y)$
- 讨论:
- 与维纳滤波器相比,约束最小二乘方滤波器假定噪声为高斯分布,只需要噪声方差和均值的知识
- 噪声方差和均值常常可从退化图像中获得,因而该方法基本上可适用于大多数图像
- 使约束最小二乘方滤波器发挥最佳效果的重要因素,一是选择合适的$Q$;二是调整参量$γ$,使其满足约束条件
- 假定噪声与图像灰度值互不相关,可通过取样方式获得图像中的噪声均值与方差,由此即可进行最佳复原; 噪声参数获得的正确与否,是取得最佳效果的关键
- 虽然约束最小二乘方滤波器可导出经典维纳滤波器,但维纳滤波器在均方意义上的最佳滤波性能却不能在这里得到证明
- 参量$γ$的调整已研究了多种方法,如贪婪算法、模拟退火算法、共轭梯度法、牛顿-拉佛森算法等,但仍然是一个有待解决的问题
- 伪逆滤波器:
匀速运动模糊图像的复原
- 设图像$f(x,y)$有一个平面运动,令$x_0(t)$和$y_0(t)$分别为在$x$和$y$方向上运动的变化分量, $t$表示运动的时间,$T$表示照相机的快门时间
- 可有目标物与摄像机之间相对运动造成图像模糊的模型: $g(x,y)=\int^T_0 f[x-x_0(t),y-y_0(t)]dt$
- 对模型表达式两边进行傅立叶变换可得: $H(u,v)=\int^T_0 exp(-j2\pi[ux_0(t)+vy_0(t)])dt$
- 水平方向匀速直线运动时, 若原图像的宽度为$L$,总的运动位移量为$a$,总的运动时间为$T$,则运动的速率为 $x_0(t)=\frac{at}{T}$, 可有 $H(u,v)=\frac{T}{\pi ua} \sin (\pi ua)e^{-j\pi ua}$
- 对离散图像,可有如下的模型 $g(x,y)=\sum^{T-1}_{t=0} f[x-\frac{at}{T},y]\cdot B$, $B$表示单位时间内曝光度的影响因子,由于是匀速运动,$B=1/T$