图像复原程序已加G.pptVIP

数字图像处理 (Digital Image Processing） 计算机信息工程学院计算机信息工程学院 第五章 图像复原 5.1 退化及噪声 5.2 图像退化的数学模型 5.3 无约束恢复 5.4 有约束恢复 5.5 交互式恢复 5.1 退化及噪声 1. 图象复原的概念 1）图像复原的定义 图像复原也称图象恢复，是图象处理中的一大类技术。所谓图像复原， 是指去除或减轻在获取数字图像过程中发生的图像质量下降（退化）这些退 化包括由光学系统、运动等等造成图像的模糊，以及源自电路和光度学因素 的噪声。 图像复原的目标是对退化的图像进行处理，使它趋向于复原成没有退化 的理想图像。成像过程的每一个环节（透镜，感光片，数字化等等）都会引 起退化。 在进行图像复原时，既可以用连续数学，也可以用离散数学进行处理。 其次，处理既可在空间域，也可在频域进行。 2）图象恢复与图象增强的异同 相同点：改进输入图象的视觉质量 。 不同点：图象增强目的是取得较好的视觉结果(不考虑退化原因)； 图象 恢复根据相应的退化模型和知识重建或恢复原始的图象(考虑退化原因) 。 2. 图象退化的原因 图象退化指由场景得到的图象没能完全地反映场景的真实内容，产生了 失真等问题。其原因是多方面的。如： 透镜象差/色差 聚焦不准（失焦，限制了图象锐度） 模糊（限制频谱宽度） 噪声（是一个统计过程） 抖动（机械、电子） 3. 图象复原方法分类 按采用的技术可分为：无约束和有约束 按采用的策略可分为：自动和交互 按采用的处理所在域可分为：频域和空域 图像退化举例1： 图像退化举例2： 4.噪声及其特性 噪声是最常见的退化因素之一，对信号来说，噪声是一种外部干扰。 但噪声本身也是一种信号（携带了噪声源的信息）。 1）关于噪声的度量 人们常只关心噪声的强度 ，可用信噪比(signal-to-noise ratio,SNR)、 能量比(电压平方比) 等来描述。分别表示为： 2）常见噪声 热噪声：白噪声（频率覆盖整个频谱均匀） 高斯噪声（幅度符合高斯分布） 闪烁噪声：具有反比于频率（1/f）的频谱； 粉色噪声（在对数频率间隔内有相同的能量） 发射噪声：高斯分布（电子运动的随机性） 3）噪声的概率密度函数 噪声为随机变量，用概率密度来刻画。 (1) 高斯噪声 (2) 均匀噪声 其均值和方差为： (3) 脉冲噪声 噪声脉冲可以是正的或负的，一 般假设a和b。都是“饱和”值双极性脉 冲噪声也称椒盐噪声。 4）噪声的形成 高斯噪声：电子噪声、弱光照/温度条件下的传感器噪声 瑞利分布：深度成像、超声波图像 指数和Gamma分布：激光成像 椒盐噪声：快速瞬变、误切换 周期噪声：图像采集过程中的电子或电磁干扰 图像中噪声的概率密度函数举例1： 原图 图像和其直方图 原图 图像和其直方图 图像中噪声的概率密度函数举例2： 图像中噪声的概率密度函数举例3： 图像中的周期噪声 5）噪声参数估计 (1) 周期噪声的参数估计 一般可以通过图像的频谱进行估计；特殊情况下可以直接从图像中 噪声分量的周期性进行推断（简单情形）。 (2) 一般噪声参数的估计 可以根据所采用的传感器类型进行噪声分布的部分推断；通常通过 特定的成像安排进行估计当只有已采集到的图像时，一般通过图像中的 平滑区域进行PDF参数的估计。如下图： 5.2 图像退化的数学模型 1.退化模型示意图 其中H为退化过程，n(x, y)为加性噪声（统计特性已知）。 2. 系统H的基本定义 就一般而言，系统是某些元件或部件以某种方式构造而成的整体。系统 本身所具有的某些特性就构成了通过系统的输入信号与输出信号的某种联系 。系统的分类可有：线性系统和非线性系统，时变系统和非时变系统，集总 参数系统和分布参数系统，连续系统和离散系统。 1）线性系统：是具有均匀性和相加性的系统 2）时不变系统：满足各个参数不随时间变化。 3）空间不变系统：满足),(),( baba-=-yxgyxfH 实际上，大部分系统是非线性和空间变化的，但以这样的模型处理起 来困难很大，一般都简化为线性的非时变和非空间变化的近似模型进行处 理。这样近似的优点是使线性系统理论中的许多理论可以直接用来解决图 象复原问题。 3.连续函数的退化模型 设系统H对坐标为(,)处的冲激函数(x-,y-)的冲激响应为h(x,y,)， 则 此式说明，如果系统H对冲激函数的响应为已知，则对任意输入的响应可 用上式求得，即，线性系统H完全可以由冲激响应来表征。图像中冲激响应 也称为点扩散函数。 在有噪音的情况下： 4.离散函数的退化模型 对和进行均匀取样后，就可引伸出离散函数的退化模型。用一维的来 说明。如果f (x)和h(x)周期分别A和B的序列，为避免卷积周期重叠需要对它们进 行周期扩展为周期为M A + B 1。 那么它们的时域离散卷积可定义为下式： 显然，上式也是具有周期M的序列。 如果用矩阵来表示上述离散退化模型，可写成下式之形式: 用矩阵形式表示 根据周期性：he(x) = he(x+M) 由于的h(x)周期性，使得H成为一个循环矩阵。 推广到2-D 扩展： 则退化过程为： 由于的h(x，y)的周期性，使得H成为一个块循环矩阵。每一 块如下所示： (1) 长时间曝光下大气湍流造成的转移函数： 其中，C是与湍流性质有关的常数。 5.退化参数的估计 退化参数的估计包括噪声估计和点扩展函数的估计。下面主要介绍点扩 散函数的估计。 1）运用先验知识估计 大气湍流、光学系统散焦、照相机与景物相对运动等，根据导致模糊的 物理过程（先验知识）来确定h(x,y)或H(u,v)。 （2）光学散焦转移函数： 其中，d是散焦点扩展函数的直径, J1() 是第一类贝塞尔函数。 (3)照相机与景物相对运动 设T为快门时间，x0(t)，y0(t)是位移的x分量和y分量。 点扩散函数示图： (a) 光脉冲 (b) 退化的光脉冲 5.3无约束恢复 1.无约束和有约束复原 由退化模型得： 最小均方误差准则： 1) 无约束复原： 在最小二乘方意义上说，希望找到一个 使下式达到最小: 求 n 2最小等效于求 2 fHg- 最小，即求 2 ) (fHgfJ-= 的极小值问题。这里选择 f 除了要求 ) ( fJ 为最小外，不受任何 其它条件约束，因此称为非约束复原。求 ) ( fJ 的极小值方法就 是一般的求极值的方法。把 ) ( fJ 对 f 微分，并使结果为 0，即: 这种方法要求知道成象系统的表达式H。 2）有约束复原方法 在最小二乘方复原处理中，为了在数学上更容易处理，常常附加某 种约束条件。 因为 H 是一方阵，并且设H -1 存在，则可求得 :f 例如，可以令 Q 为 f 的线性算子，那么，最小二乘方复原问题可看 成是使形式为 2 f Q 的函数，服从约束条件 的最小化 2 2 fHgn-= 问题。而这种有附加条件的极值问题可用拉格朗日乘数法来处理。 其处理方法如下： 寻找一个f ，使下述准则函数为最小 ) ( )( 2 22 nfHgfQfJ-+= 式中 为一常数，是拉格朗日系数。加上约束条件后，就可以按一般 求极小值的方法进行求解。将上式对 微分，并使结果为零，则有: 2. 逆滤波 1）逆滤波原理 设M = N ，则： 退化函数H (u, v)与F (u, v)相乘为退化过程，用H (u, v)去除G (u, v) 是复 原过程，称其为逆滤波。可描述为： 记M (u, v)为复原转移函数，则其等于1 / H (u, v)。 2）分析/讨论 由逆滤波知： (1) H (u, v)在UV 平面上取零或很小，N (u, v) / H (u, v)就会使恢复结果与 预期的结果有很大差距; (2) 噪声带来更严重的问题（知道H也估计不准 f ） H (u, v)常随u，v与原点距离的增加而迅速减小，而噪声N (u, v)却一般 变化缓慢。在这种情况下，恢复只能在与原点较近（接近频域中心）的范围 内进行。 3）改进： 去除零点，改进为： 消除振铃现象，改进为： 4) 图象退化和恢复模型： 将点源图象看做单位脉冲函数（F (x, y) = 1）的近似则有: G(u, v) = H(u, v) F(u, v) H(u, v) 图象退化和恢复示例 退化图 滤波器 除去零点 减少振铃 3. 消除匀速直线运动模糊 匀速直线运动时： T: 采集时间长度。 x方向运动分量 y方向运动分量 对水平方向匀速直线运动 设x0(t) = ct / T ，y0(t) = 0 当n为整数时，H在u = n/c处为零。 当 f (x, y)在区间0 x L之外为零或已知时： 举例： 1.维纳（Wiener）滤波器 它一种最小均方误差滤波器。 设 Rf 是 f 的相关矩阵： Rf 的第 ij 个元素是Efi fj，代表 f 的第 i 和第 j 元素的相关。 设 Rn是n 的相关矩阵： 根据两个象素间的相关只是它们相互距离而不是位置的函数的假设，可 将 Rf 和 Rn 都用块循环矩阵表达，并借助矩阵W来对角化： 5.4有约束恢复 fe(x, y)的功率谱，记为Sf (u, v) ；ne(x, y)的功率谱，记为Sn(u, v) 。 D是1个对角矩阵，D(k, k) = (k) ，则有： 定义： 代入： 则有： 两边同乘以W 1，有： 原图 退化图像 全逆滤波 半径受限逆滤波 维纳滤波结果 复原举例1： 运动模糊 复原举例2： (a) 运动模糊和加性噪声图像 (b) (b) 逆滤波复原 (c) 维纳滤波复原 (d), (e), (f) 顺序同上，但其中的噪声幅值降低一个数量级 (g), (h), (i) 顺序同上，但其中的噪声幅值降低五个数量级 2.有约束最小平方恢复 只需有关噪声均值和方差的知识就可对每个给定图象得到最优结果 （仍需确定变换矩阵Q）： 建立基于平滑测度上的最优准则 ，设 f (x, y)在(x, y)处的二阶微分为： 其模板为： 扩展为： 最优准则: 循环矩阵为： 其中： 对之进行对角化为： 其中E是对角矩阵，它的元素为 : P(u, v)是pe(x, y)的2-D傅里叶变换；k / N代表不超过k/N的最大的整数； k mod N代表用N除k得到的余数 。 约束为： ，最优解为： 有约束最小平方恢复举例： 上图是最小平方恢复的结果 下图是维纳滤波的结果 %P118用wiener去除运动造成的模糊 I=imread(LENA256.bmp); subplot(2,2,1);imshow(I); title(LENA256原图); %设置运动位移为30个像素 LEN=30; %设置运动角度为75度 THETA=75; %建立二维仿真线性运动滤波器PSF PSF=fspecial(motion,LEN,THETA); %用PSF产生退化图像 MF=imfilter(I,PSF,circular,conv); subplot(2,2,2);imshow(MF); title(运动模糊的图像); %用Wiener滤波消除运动模糊的图像 wnr=deconvwnr(MF,PSF); subplot(2,2,4);imshow(wnr); title(消除模糊后图像); %P122最小二乘滤波复原程序 I=imread(LENA256.bmp); %读取原始图像 LEN=31; %图像的模糊化 THETA=11; PSF1=fspecial(motion,LEN,THETA); PSF2=fspecial(gaussian,10,5); Blurred1=imfilter(I,PSF1,circular,conv); Blurred2=imfilter(I,PSF2,conv); %模糊化图像加噪 V=.002; BlurredNoisy1=imnoise(Blurred1,gaussian,0,V); BlurredNoisy2=imnoise(Blurred2,gaussian,0,V); subplot(2,2,1);imshow(I);title(LENA256.bmp) %用真实PSF函数和噪声强度作为参数进行图像复原 NP=V*prod(size(I); reg1=deconvreg(BlurredNoisy1,PSF1,NP); reg2=deconvreg(BlurredNoisy2,PSF2,NP); subplot(2,2,2);imshow(reg1); title(Restored1 with NP) subplot(2,2,4) ;imshow(reg2); title(Restored2 with NP) %P123补Lucy_Richardson滤波复原 I=imread(board.tif);I=I(50+1:256,2+1:256,:); subplot(2,2,1);imshow(I);title(Original Image); PSF=fspecial(gaussian,5,5); Blurred=imfilter(I,PSF,symmetric,conv);V=.002; BlurredNoisy=imnoise(Blurred,gaussian,0,V); subplot(2,2,2);imshow(BlurredNoisy);title(Blurred luc1=deconvlucy(BlurredNoisy,PSF,5); subplot(2,2,3);imshow(luc1);title(Restored Image,NUMIT=5); luc2=deconvlucy(BlurredNoisy,PSF,10); subplot(2,2,4);imshow(luc2);title(Restored Image,NUMIT=10较好); figure; subplot(1,3,3) ;imshow(wnr2); title(Restored2,True PSF) BlurredNoisy1=imnoise(Blurred1,gaussian,0,V); BlurredNoisy2=imnoise(Blurred2,gaussian,0,V); wnr3=deconvwnr(Blurred1,PSF1); wnr3=deconvwnr(Blurred1,fspecial(motion,2*LEN,THETA); wnr4=deconvwnr(Blurred1,fspecial(motion,LEN,2*THETA); figure; subplot(1,2,1);imshow(wnr3); title(Restored1,True PSF); figure; subplot(1,2,2);imshow(wnr2); title(Restored2,True PSF) NP=(V*prod(size(I).2; NPOW=su