5图像复原程序已加G[1]解析课件.ppt

1、数字图像处理数字图像处理(Digital Image Processing）第五章第五章 图像复原图像复原 5.1 退化及噪声退化及噪声 5.2 图像退化的数学模型图像退化的数学模型 5.3 无约束恢复无约束恢复 5.4 有约束恢复有约束恢复 5.5 交互式恢复交互式恢复5.1 退化及噪声退化及噪声1.图象复原的概念图象复原的概念1）图像复原的定义）图像复原的定义 图像复原也称图象恢复，是图象处理中的一大类技术。所谓图像复原，是指去除或减轻在获取数字图像过程中发生的图像质量下降（退化）这些退化包括由光学系统、运动等等造成图像的模糊，以及源自电路和光度学因素的噪声。图像复原的目标是对退化的图像进

2、行处理，使它趋向于复原成没有退化的理想图像。成像过程的每一个环节（透镜，感光片，数字化等等）都会引起退化。在进行图像复原时，既可以用连续数学，也可以用离散数学进行处理。其次，处理既可在空间域，也可在频域进行。2）图象恢复与图象增强的异同）图象恢复与图象增强的异同 相同点：相同点：改进输入图象的视觉质量。不同点：不同点：图象增强目的是取得较好的视觉结果(不考虑退化原因)；图象恢复根据相应的退化模型和知识重建或恢复原始的图象(考虑退化原因)。2.图象退化的原因图象退化的原因 图象退化指由场景得到的图象没能完全地反映场景的真实内容，产生了失真等问题。其原因是多方面的。如：透镜象差/色差 聚焦不准（失

3、焦，限制了图象锐度）模糊（限制频谱宽度）噪声（是一个统计过程）抖动（机械、电子）3.图象复原方法分类图象复原方法分类 按采用的技术技术可分为：无约束和有约束 按采用的策略策略可分为：自动和交互 按采用的处理所在域处理所在域可分为：频域和空域 图像退化举例图像退化举例1：图像退化举例图像退化举例2：4.噪声及其特性噪声及其特性 噪声是最常见的退化因素之一，对信号来说，噪声是一种外部干扰。但噪声本身也是一种信号（携带了噪声源的信息）。1 1）关于噪声的度量）关于噪声的度量 人们常只关心噪声的强度，可用信噪比(signal-to-noise ratio,SNR)、能量比(电压平方比)等来描述。分别表

4、示为：2）常见噪声）常见噪声 热噪声：白噪声（频率覆盖整个频谱均匀）高斯噪声（幅度符合高斯分布）闪烁噪声：具有反比于频率（1/f）的频谱；粉色噪声（在对数频率间隔内有相同的能量）发射噪声：高斯分布（电子运动的随机性）2210log10nsVVSNR22ob噪声均方差灰度对比度CSNR3）噪声的概率密度函数）噪声的概率密度函数 噪声为随机变量，用概率密度来刻画。(1)高斯噪声(2)均匀噪声其均值和方差为：222)(exp21)(zzp其他0)/(1)(bzaabzp如果2/)(ba 12/)(22ab(3)脉冲噪声 噪声脉冲可以是正的或负的，一般假设a和b。都是“饱和”值双极性脉冲噪声也称椒盐噪

5、声。4）噪声的形成）噪声的形成 高斯噪声：电子噪声、弱光照/温度条件下的传感器噪声 瑞利分布：深度成像、超声波图像 指数和Gamma分布：激光成像 椒盐噪声：快速瞬变、误切换 周期噪声：图像采集过程中的电子或电磁干扰其他0)(bzPazPzpba如果如果图像中噪声的概率密度函数举例图像中噪声的概率密度函数举例1：原图原图图像和其直方图图像和其直方图原图原图图像和其直方图图像和其直方图图像中噪声的概率密度函数举例图像中噪声的概率密度函数举例2：图像中噪声的概率密度函数举例图像中噪声的概率密度函数举例3：图像中的周期噪声图像中的周期噪声5）噪声参数估计）噪声参数估计(1)周期噪声的参数估计 一般可

6、以通过图像的频谱进行估计；特殊情况下可以直接从图像中噪声分量的周期性进行推断（简单情形）。(2)一般噪声参数的估计 可以根据所采用的传感器类型进行噪声分布的部分推断；通常通过特定的成像安排进行估计当只有已采集到的图像时，一般通过图像中的平滑区域进行PDF参数的估计。如下图：5.2 图像退化的数学模型图像退化的数学模型1.1.退化模型示意图退化模型示意图其中H为退化过程，n(x,y)为加性噪声（统计特性已知）。2.系统系统H的基本定义的基本定义 就一般而言，系统是某些元件或部件以某种方式构造而成的整体。系统本身所具有的某些特性就构成了通过系统的输入信号与输出信号的某种联系。系统的分类可有：线性系

7、统和非线性系统，时变系统和非时变系统，集总参数系统和分布参数系统，连续系统和离散系统。1）线性系统：是具有均匀性和相加性的系统 2）时不变系统：满足各个参数不随时间变化。3）空间不变系统：满足x,yH()x,yf()n()x,yg),(),(babayxgyxfH 实际上，大部分系统是非线性和空间变化的，但以这样的模型处理起来困难很大，一般都简化为线性的非时变和非空间变化的近似模型进行处理。这样近似的优点是使线性系统理论中的许多理论可以直接用来解决图象复原问题。3.连续函数的退化模型连续函数的退化模型 设系统H对坐标为(a,b)处的冲激函数(x-a,y-b)的冲激响应为h(x,a,y,b)，则

8、此式说明，如果系统H对冲激函数的响应为已知，则对任意输入的响应可用上式求得，即，线性系统H完全可以由冲激响应来表征。图像中冲激响应也称为点扩散函数。在有噪音的情况下：bababaddyxhfyxg),(),(),(),(),(),(),(yxnddyxhfyxg bababa4.4.离散函数的退化模型离散函数的退化模型 对和进行均匀取样后，就可引伸出离散函数的退化模型。用一维的来说明。如果f(x)和h(x)周期分别A和B的序列，为避免卷积周期重叠需要对它们进行周期扩展为周期为M A+B 1。那么它们的时域离散卷积可定义为下式：显然，上式也是具有周期M的序列。如果用矩阵来表示上述离散退化模型，可

9、写成下式之形式:1010)()(eMxAAxxfxffHg 1010)()(eMxBBxxhxh10eee1,1 ,0)()()(MmMxmxhmfxg用矩阵形式表示用矩阵形式表示根据周期性：he(x)=he(x+M)1()1()0()0()2()1()2()0()1()1()1()0()1()1()0(eeeeeeeeeeeeeeeMfffhMhMhMhhhMhhhMgggHfg)0()2()1()2()0()1()1()1()0(eeeeeeeeehMhMhhhhhMhhH由于的h(x)周期性，使得H成为一个循环矩阵。推广到2-D 扩展：则退化过程为：1101010),(),(eNyBMx

10、AByAxyxfyxf或和1101010),(),(eNyDMxCDyCxyxhyxh或和1,1 ,01,1 ,0),(),(),(1010eee NyMxnymxhnmfyxgMmNn)1()1()0()1()1()0(eeeeee021201110MNnnnMNfffMMMHHHHHHHHHnHfg)1()1()0()1()1()0(eeeeee021201110MNnnnMNfffMMMHHHHHHHHHnHfg)0,()2,()1,()2,()0,()1,()1,()1,()0,(eeeeeeeeeihNihNihihihihihNihihiH 由于的h(x，y)的周期性，使得H成为一

11、个块循环矩阵。每一块如下所示：(1)长时间曝光下大气湍流造成的转移函数：exp),(6/522vucvuH其中，C是与湍流性质有关的常数。5.5.退化参数的估计退化参数的估计 退化参数的估计包括噪声估计和点扩展函数的估计。下面主要介绍点扩散函数的估计。1）运用先验知识估计 大气湍流、光学系统散焦、照相机与景物相对运动等，根据导致模糊的物理过程（先验知识）来确定h(x,y)或H(u,v)。（2）光学散焦转移函数：2/1221)()(),(vuddJvuH其中，d是散焦点扩展函数的直径,J1()是第一类贝塞尔函数。(3)照相机与景物相对运动 设T为快门时间，x0(t)，y0(t)是位移的x分量和y

12、分量。点扩散函数示图：dttvytuxjvuHT)()(2exp),(000(a)光脉冲(b)退化的光脉冲5.3无约束恢复无约束恢复 Hfgn)()(T2T2fHgfHgfHgnnn1.无约束和有约束复原 由退化模型得：最小均方误差准则：1)无约束复原：无约束复原：在最小二乘方意义上说，希望找到一个 使下式达到最小:f22fHgn求n2最小等效于求2fHg最小，即求2)(fHgfJ的极小值问题。这里选择f除了要求)(fJ为最小外，不受任何其它条件约束，因此称为非约束复原。求)(fJ的极小值方法就是一般的求极值的方法。把)(fJ对 f微分，并使结果为 0，即:0)(2)(fHgHffJgHfHH

13、gHHHf1)(gHgHHHf111)(这种方法要求知道成象系统的表达式H。2）有约束复原方法）有约束复原方法 在最小二乘方复原处理中，为了在数学上更容易处理，常常附加某种约束条件。因为 H 是一方阵，并且设H-1存在，则可求得 :f 例如，可以令 Q 为 f 的线性算子，那么，最小二乘方复原问题可看成是使形式为2fQ的函数，服从约束条件 的最小化22fHgn问题。而这种有附加条件的极值问题可用拉格朗日乘数法来处理。其处理方法如下：0)(22)(fHgHfQQffJgHfHHfQQ1gHQQHHf1)1(寻找一个f，使下述准则函数为最小)()(222nfHgfQfJ式中 为一常数，是拉格朗日系

14、数。加上约束条件后，就可以按一般求极小值的方法进行求解。将上式对 微分，并使结果为零，则有:f2.逆滤波逆滤波 1）逆滤波原理设M=N，则：退化函数H(u,v)与F(u,v)相乘为退化过程，用H(u,v)去除G(u,v)是复原过程，称其为逆滤波。可描述为：gWWDgWDWgHf11111)(gWDfW1111,1 ,0,),(),(),(),(11MyxvuHvuGvuFyxfFF1,1 ,0,),(),(),(MvuvuHvuGvuF记M(u,v)为复原转移函数，则其等于1/H(u,v)。2 2）分析）分析/讨论讨论由逆滤波知：(1)H(u,v)在UV 平面上取零或很小，N(u,v)/H(u

15、,v)就会使恢复结果与预期的结果有很大差距;(2)噪声带来更严重的问题（知道H也估计不准 f）H(u,v)常随u，v与原点距离的增加而迅速减小，而噪声N(u,v)却一般变化缓慢。在这种情况下，恢复只能在与原点较近（接近频域中心）的范围内进行。3）改进：）改进：去除零点，改进为：消除振铃现象，改进为：1,1 ,0,),(),(),(),(MvuvuHvuNvuFvuF20222022 1 ),(1),(wvuwvuvuHvuM如如其它如),(1),(),(vuHdvuHkvuM4)图象退化和恢复模型：图象退化和恢复模型：将点源图象看做单位脉冲函数（F (x,y)=1）的近似则有:G(u,v)=H

16、(u,v)F(u,v)H(u,v)图象退化和恢复示例图象退化和恢复示例 退化图 滤波器 除去零点 减少振铃fM()(x,yn()x,yH()u,vf()x,yg()x,yu,v3.消除匀速直线运动模糊消除匀速直线运动模糊 匀速直线运动时：),(),()d()(j2exp),(d dd )(j2)exp(),(-d)d(j2)exp,(),(000000vuFvuHttvytuxvuFtyxvyuxtyytxxfyxvyuxyxgvuGTT Tttyytxxfyxg000d )(),(),(T:采集时间长度。x方向运动分量 y方向运动分量 对水平方向匀速直线运动对水平方向匀速直线运动 设x0(t

17、)=ct/T，y0(t)=0 当n为整数时，H在u=n/c处为零。当 f(x,y)在区间0 x L之外为零或已知时：举例：ucucucTtTctuvuHTjexp)sin(d j2exp)(0,LxftTctxfxgxcxT0)d(d )(01.维纳（维纳（Wiener）滤波器）滤波器 它一种最小均方误差滤波器。设 Rf 是 f 的相关矩阵：Rf 的第 ij 个元素是Efi fj，代表 f 的第 i 和第 j 元素的相关。设 Rn是n 的相关矩阵：根据两个象素间的相关只是它们相互距离而不是位置的函数的假设，可将 Rf 和 Rn 都用块循环矩阵表达，并借助矩阵W来对角化：gHRRHHgHQQHH

18、fT1 1TT1 TT nfssTffRfETnnRnE5.4有约束恢复有约束恢复1WAWRf1WBWRnfe(x,y)的功率谱，记为Sf(u,v)；ne(x,y)的功率谱，记为Sn(u,v)。D是1个对角矩阵，D(k,k)=(k)，则有：定义：代入：则有：两边同乘以W 1，有：1WDWHgHQQHHfT1 TT snfRRQQ1T1WDWH1TW*WDHgHRRHHfnfT11T)(sgWWDBWWADWWDf11111*)*(sgWDBADDfW1111*)*(s1WAWRf1WBWRn221|(,)|(,)(,)(,)|(,)|(,)/(,)fH u vF u vG u vH u vH

19、u vS u vSu v原图 退化图像 全逆滤波 半径受限逆滤波 维纳滤波结果复原举例复原举例1：运动模糊复原举例复原举例2：(a)运动模糊和加性噪声图像 (b)(b)逆滤波复原 (c)维纳滤波复原(d),(e),(f)顺序同上，但其中的噪声幅值降低一个数量级(g),(h),(i)顺序同上，但其中的噪声幅值降低五个数量级2.有约束最小平方恢复有约束最小平方恢复 只需有关噪声均值和方差的知识就可对每个给定图象得到最优结果（仍需确定变换矩阵Q）：建立基于平滑测度上的最优准则，设 f(x,y)在(x,y)处的二阶微分为：其模板为：扩展为：2222(1,)(1,)(,1)(,1)4(,)ffxyf x

20、yf xyf x yf x yf x ygHQQHHfT1 TT s010(,)141010p x y131302020),(),(eNyMxyxyxpyxp或和最优准则:循环矩阵为：其中：22222),(),(minyyxfxyxf021201110CCCCCCCCCCMMM)0,()2,()1,()2,()0,()1,()1,()1,()0,(eeeeeeeeejpNjpNjpjpjpjpjpNjpjpjC对之进行对角化为：对之进行对角化为：其中E是对角矩阵，它的元素为:P(u,v)是pe(x,y)的2-D傅里叶变换；k/N代表不超过k/N的最大的整数；k mod N代表用N除k得到的余数

21、。约束为：，最优解为：CWWE1kikiNkNkPikE如如0mod ,),(22 nfHg1,1 ,0,),(),(),(),(),(22*MvuvuGvuPsvuHvuHvuFgWDEEDDfW1*1*1)(sgWWDEWWEDWWDgHCCHHf1*11*1*T1TT)()(ss有约束最小平方恢复举例：有约束最小平方恢复举例：上图是最小平方恢复的结果下图是维纳滤波的结果%P118用用wiener去除运动造成的模糊去除运动造成的模糊I=imread(LENA256.bmp);subplot(2,2,1);imshow(I);title(LENA256原图原图);%设置运动位移为设置运动位移

22、为30个像素个像素LEN=30;%设置运动角度为设置运动角度为75度度THETA=75;%建立二维仿真线性运动滤波器建立二维仿真线性运动滤波器PSFPSF=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(消除模糊后图像消除模

23、糊后图像);%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);Blurre

24、dNoisy2=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);

25、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);tit

26、le(Blurred&Noisy);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(Blurre

27、d1,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(

28、Restored2,True PSF)NP=(V*prod(size(I).2;NPOW=sum(NP(:)/prod(size(I);%噪声功率NCORR=fftshift(real(ifftn(NP);%噪声自相关函数（ACF）IP=abs(fftn(im2double(I).2;IPOW=sum(IP(:)/prod(size(I);%原始图像的功率ICORR=fftshift(real(ifftn(IP);%图像自相关函数（ACF）wnr5=deconvwnr(BlurredNoisy1,PSF1,NCORR,ICORR);wnr6=deconvwnr(BlurredNoisy2,PS

29、F2,NCORR,ICORR);figure;subplot(1,2,1);imshow(wnr5);subplot(1,2,2);imshow(wnr6)人机结合控制恢复过程以达到一些特殊的效果。如正弦干扰模式如正弦干扰模式（相关噪声）：只有虚分量，它代表一对位于频率平面上坐标分别为(u0/2,v0/2)和(u0/2,v0/2)，强度分别为 A/2和A/2的脉冲：可依靠视觉观察在频率域确定出脉冲分量的位置并在 该位置利用带阻滤波器消除；存在多个正弦分量时，在频率域里对应每个亮点的位置放1个带通滤波器H(u,v)加以滤波。5.5交互式恢复交互式恢复)sin(),(00yvxuAyx2,2 2,

30、2 2j),(0000vvuuvvuuAvuN空域相对应的结构模式：空域相对应的结构模式：从g(x,y)中减去加权的p(x,y)（其中w(x,y)称为权函数），有：点(x,y)邻域的均值和方差（最小化）),(),(),(),(yxpyxwyxgyxf),(),(),(1vuGvuHyxpF XXmYYnnymxfYXyxf),()12)(12(1),(22),(),()12)(12(1),(XXmYYnyxfnymxfYXyx),(),(),(vuGvuHvuP设w(x,y)在邻域中基本是常数能最小化 2(x,y)的w(x,y)为：222 ),(),(),(),(),(),()12)(12(1 ),(),(),(),(),(),()12)(12(1),(yxpyxwyxgnymxpyxwnymxgYXyxpyxwyxgnymxpnymxwnymxgYXyxXXmYYnXXmYYn),(),(),(),(),(),(),(22yxpyxpyxpyxgyxpyxgyxw谢谢 谢谢