数字图象处理第五章

1、第五章 图像的恢复与重构图像恢复与图像增强的相同和不同之处图像恢复与图像增强的相同和不同之处两者都是为了改进输入图像的视觉质量。图像增强借助人的视觉系统特性，以取得较好的视觉结果(不考虑退化原因)；而图像恢复根据相应的退化模型和知识重建或恢复原始的图像(考虑退化原因)。什么是图像退化：什么是图像退化：图像通过“系统”后的质量变坏叫做退化。退化的形式有图像模糊、图像有干扰等。图像退化的处理方法：图像退化的处理方法：如果我们对退化的类型、机制和过程都十分清楚，那么就可以利用其反过程来复原图像。规则图案变形，胶片冲洗时易发生边缘模糊， 光学系统中的孔径衍生产生退化运动模糊，或在拍摄过程中相机发生振动

2、随机噪声的叠加典型的图像恢复方法是根据图像退化的先验知识建立一个退化模型，以此模型为基础，采用滤波等手段进行处理，使得复原后的图像符合一定的准则，达到改善图像质量的目的。close all;clear all;I = imread(cameraman.tif);figure;imshow(I);title(Original Image); H = fspecial(motion,20,45);figure;mesh(H);MotionBlur = imfilter(I,H); %replicatefigure;imshow(MotionBlur);title(Motion Blurred Im

3、age);H = fspecial(disk,10); %replicatefigure;mesh(H);blurred = imfilter(I,H); figure;imshow(blurred);title(Blurred Image);运动模糊与边缘模糊close all;clear all;img=imread(noise.jpg);I=rgb2gray(img);figure;imshow(I);figure;imhist(I);J = imnoise(I, gaussian,0.05,0.01);%均值、方差均值、方差%J = imnoise(I, poisson);%均值为均值

4、为10%J = imnoise(I, salt & pepper,0.2);%噪声密度噪声密度%J = imnoise(I, speckle,0.05);%方差；方差；J= I+n*I;n为一致性噪声点为一致性噪声点figure;imshow(J);figure;imhist(J);通过直方图通过直方图对噪声的再认识对噪声的再认识f(i, j)：原始图像g(i, j)：降质图像T()：成像系统的作用，则： g(x,y)=Tf(x,y)设T是线性移不变的。一幅连续的图像f(x,y)可以用函数的二维卷积表示： ddyxfyxf),(),(),(因此，令h(x,;y,) =T (x-,y-)

5、 ，则有： ddyxTfyxg),(),(),( ddyxhfyxg),;,(),(),(定义于不在原点的二维函数由于f(, )与x，y没有关系称 h(x,;y,) 为点扩散函数(PSF)或系统冲击响应。多数情况下它表现为时不变的，反映在图像中为位移不变的，则 h(x,;y,) 可以表示为h(x-,y-) ),(),(),(),(),(yxhyxfddyxhfyxg 其中*表示卷积运算。如果T()是一个h可分离系统，即： h(x,;y,)=h1 (x,) h2 (y, )则二维运算可以分解为列和行两次一维运算来代替。在加性噪声情况下，图像退化模型可以表示为： g(x,y)=f(x,y)*h(x

6、,y)+n(x,y)其中n(x, y)为噪声图像。 对于图像降质过程进行数学建模，设：f(i, j)为原始图像；y(i, j)为降质图像；h(i, j; k, l)为点扩散函数；图像为MN维。有 MkNljinlkflkjihjiy11),(),(),;,(),(假设为空间移不变h(i, j; k, l)，则：),(),(),(),(),(),(),(11jinjifjihjinlkfljkihjiyMkNl线性位移不变的图像退化模型则表示为： g(x,y)=f(x,y)*h(x,y)+n(x,y)推测结论：如果已知g(x,y)、 n(x,y)、 h(x,y)，则 f(x,y)可以计算出来。对

7、等式两端取傅立叶变换有： G(u,v)=F(u,v)H(u,v)+N(u,v) F(u,v)= (G(u,v)- N(u,v)/ H(u,v)f(x,y)= F-1F(u,v)= F-1G(u,v) / H(u,v)- N(u,v)/ H(u,v)结论的正确与否需要证明(对卷积定理的另一种理解)。g(x,y)f(x,y)n(x,y)h(x,y)+三、循环矩阵及傅立叶化一个一维离散序列通过一个系统发生失真的过程可用下图表示1210 )()()(10,.,M-,xn(x)mxhmfxgMmeee) 1() 1 () 0() 1() 1 () 0() 0()2() 1()2() 0() 1 () 1

8、() 1() 0() 1() 1 () 0(MnnnMfffhMhMhMhhhMhhhMgggeeeeeeeeeeeeeeeeeenHfgf= Hf= H-1-1g- ng- n用矩阵表示，可以写成 g(x)f(x)n(x)h(x)+如果考虑加性噪声，根据离散序列的卷积定理，有扩展为周期为扩展为周期为M的序列的序列由于离散卷积的周期性，有he(x)=he(x+M)，H可以写成) 0 () 2() 1() 2() 0 () 1 () 1 () 1() 0 (eeeeeeeeehMhMhhhhhMhhH10101-N0,1,2,.,y1-M0,1,2,.,x ),(),(),(),(MmNneee

9、eyxnnymxhnmfyxgH是一个循环阵。结论：一维离散卷积：单位取样响应构成的循环矩阵输入向量 对数字图象的二维离散函数也是如此。对图像退化模型而言，有A=5B=4M=9M=9用矩阵形式表示上式：g=Hf+ng=Hf+n g、f和n分别表示MN的函数矩阵ge(i, j)、fe(i, j)和ne(i, j)的各行前后相连而成的列矢量(堆叠矢量)。如果假设原始图像是MN维矩阵，则H是MNMN循环矩阵，且H是一个分块(MM个)循环矩阵：) 1() 1 ()0() 1() 1 ()0() 1() 1 ()0( 021201110MNnnnMNfffMNgggeeeeeeMMMeeeHHHHHHH

10、HHnHfg每一个子矩阵Hi自身也是循环矩阵NN：)0,()2,()1,()2,()0,()1 ,()1 ,()1,()0,(ihNihNihihihihihNihiheeeeeeeeeiH=+MN1 MNMN MN1 MN1 1、一维信号序列循环矩阵的对角化和傅立叶化解矩阵方程： f= Hf= H-1-1g- ng- n最简单的计算方法就是对角化，H H H H-1-1也是对角阵。对角化H H的方法求取其特征值和特征矢量。对循环矩阵而言，设：其有M个特征值和特征矢量。110 ) 1(2exp2exp1)(M-, kkMMjkMjkTwkMMjhkMjMhhkeee) 1(2exp) 1 (2

11、exp) 1()0()(由于w w(k)是由傅立叶系数构成的，因此w w(k)彼此是正交的。所以，由w w(k)构成的变换矩阵是可逆的。102)()(NuuxNjeuFxf用特征矢量组成的矩阵：WW=w w(0),w w(1)w w(M-1)生成对角矩阵D D：D=W D=W -1-1HWHW；且D(k,k)=(k)。而kiMj(k,i)wkiMjw(k,i)kMMjkMjkMHHHMMhhhMMH(k)ikMjihMMkiMMjihMMkMMjhkMjMhhkeeeMieMieeee2exp 2exp ) 1(2exp2exp1)( ) 1(000) 1 (000)0( )1() 1 (),

12、0( 2exp)(1)(2exp)(1 ) 1(2exp) 1 (2exp) 1()0()(11010为正交阵序列的傅立叶变换乘WwHWWDT1根据周期性，根据周期性，M i = - i凑一个凑一个常数常数所以所以 H=WDW H=WDW -1 -1 根据周期性，根据周期性，M i = - i对g=Hf+ng=Hf+n而言，可以写成g=WDWg=WDW-1-1f+nf+n ，有WW1 1g=DWg=DW-1-1f+ f+ WW1 1n n；其中对 WW-1-1g g列矢量的每一行G(k)而言，有序列的傅立叶变换)1() 1 (),0(2exp)(1 2exp1) 1(2exp1) 1 (1)0

13、()(10MgggkiMjigMkiMjMMgkMjMgMgkGeeeMieeee序列的傅立叶变换)1() 1 (),0(2exp)(1 2exp1) 1(2exp1) 1 (1)0()(10MfffkiMjifMkiMjMMfkMjMfMfkFeeeMieeee对WW-1-1f f列矢量的每一行F(k)而言对WW-1-1n n而言有同样的结果。所以对WW1 1g=DWg=DW-1-1f+ Wf+ W1 1n n而言，G(u)=MH(u)F(u)+N(u)G(u)=MH(u)F(u)+N(u)。上面的过程称之为循环矩阵的傅立叶化(对卷积定理的另一种理解) 。如果图像的g g、f f、n n采用

14、堆叠矢量的方法构成， g=Hf+ng=Hf+n。同一维的情况类似，不同的地方是H为块循环矩阵，以及其中的傅立叶变换是二维的，但最后结论是一样的。 G(u)=MH(u)F(u)+N(u) 一维情况下的结论G(u,v)=MNH(u,v)F(u,v)+N(u,v)二维情况下的结论在实际应用中，认为f(x,y)、h(x,y)、g(x,y)、n(x,y)的维数是相等的。因此有最后的结论：F(u,v) =G(u,v)-N(u,v)/ MNH(u,v) f(x,y) =FF(u,v) =G(u,v)-N(u,v)/ MNH(u,v) f(x,y) =F-1-1F(u,v)F(u,v)=+ 要知道一个图象降质

15、系统的H(u,v)是一件非常困难的事情。但因为f(x,y)*h(x,y)=g(x,y)，有 F(u,v)H(u,v)=G(u,v) 如果f(x,y)=(x,y)，F(x,y)=1则H(u,v)=G(u,v)所以，可以用实验的方法得到h(x,y)和H(u,v)；H(u,v)可用点光源的输出图像的傅立叶变换来近似。另外，有一些图象降质系统的H(u,v)有固定的或近似的数学模型。n一维运动模糊：通常在拍摄过程中，相机或物体移动造成的运动模糊可以用一维均匀邻域像素灰度的平均值来表示： n大气扰动模糊：这种模糊经常出现在遥感和航空摄影中，由于曝光时间过长引起的模糊可用高斯点扩散函数来表示： 式中K是一个

16、归一化常数，保证模糊的大小为单位值，2可以决定模糊的程度。其他 022 1)(LxLLxh)2exp(),(222yxKyxhn均匀不聚焦模糊 这是由于相机聚焦不准确引起的，虽然不聚焦由许多参数决定，如相机的焦距、相机光圈的大小、形状、物体和相机之间的距离等，但在研究中为了简单起见，我们用下列函数表示聚焦不准引起的模糊：n 均匀二维模糊 这是常见的一种模糊，可以用来近似聚焦不准引起的模糊：其他 0 1),(222RyxRyxh其他 02,2 1),(2LyxLLyxhclose all;clear all;img=imread(gg.jpg);img_gray=rgb2gray(img);fi

17、gure;imshow(img_gray);len = 50;theta = -45;sigma=10;psf = fspecial(average,5); %psf = fspecial(gaussian,21,21,sigma); %psf = fspecial(motion,len,theta); % 建立建立PSFfigure;freqz2(psf); Blurred = imfilter(img_gray,psf);figure; imshow(Blurred); title(Blurred Image);常见模糊演示请特别注意点扩散函数的内容逆滤波对于图像退化模型： g(x,y)=

18、f(x,y)*h(x,y)+n(x,y)两边取傅立叶变换：G(u,v)= F(u,v) H(u,v) +N(u,v) H(u,v)又称为系统的转移函数(或滤波函数)，它使图像退化。在无噪声的情况下，上式可以简化为： G(u,v)= F(u,v) H(u,v) F(u,v) =G(u,v)/H(u,v) 这种1/H(u,v)的形式称为逆滤波。再进行傅立叶逆变换就可以得到f(x,y)。 当对噪声一无所知时，使n n无约束的小。由于g=Hf+ng=Hf+n，假设通过恢复可以得到一个不错的f f的估计f f。显然，f f应满足关系g gHfHfn n，我们希望n n尽可能的小。于是问题转化为f f在什

19、么情况下n n最小对矩阵而言就是它的迹(对角线之和)的平方最小。gHffHgHff) fH)(gfHtr(gfHgfTT12 0)(2)( )(minJJ ),(),(),( ),(),(),( 11vuHvuGvuFvuMNHvuGvuFgWDfWgWWDg)(WDWgHf-1-1-1111将MN乘入Hclose all;clear all;img=imread(gg.jpg);gray_img=rgb2gray(img);figure;imshow(gray_img);img_m=double(gray_img)*double(gray_img);figure;mesh(img_m); %

20、自相关矩阵自相关矩阵img_n1 = imnoise(gray_img,salt & pepper, 0.3);%噪声的密度figure;imshow(img_n1);img_m1=double(img_n1)*double(img_n1);figure;mesh(img_m1);对角线之和的平方最小如果如果H(u,v)H(u,v)有许多值为零的点有许多值为零的点，必然使得复原的结果，必然使得复原的结果受到极大影响。受到极大影响。 或者如果或者如果H(u,v)H(u,v)不为零但是有非常小不为零但是有非常小的值，也即病态条件，也会使复原效果受到影响。的值，也即病态条件，也会使复原效果受

21、到影响。解决这个问题的方法是避开解决这个问题的方法是避开H(u,v)H(u,v)的零点。的零点。 幸好一般幸好一般的的H(u,v)H(u,v)在在低频附近低频附近的有限区域内不为零。的有限区域内不为零。因此逆滤波可以在原点附近进行，相当于在频域乘上因此逆滤波可以在原点附近进行，相当于在频域乘上一低通窗口函数一低通窗口函数W(u,v)W(u,v)。 ),(),(),(),(),(),(),(),(),(vuHvuNvuFvuHvuNvuHvuFvuHvuF ),(),(),( ),(),(),(),( vuHvuGvuFvuNvuFvuHvuG代入将为了防止随着u、v的增大H(u,v)的迅速减小

22、而增设一些条件由于截断地原因，被恢复的图象振铃较大。一种改进的方法是取202220221 )(1wv u wv u u,v/HM(u,v)其它 1 )( )(1 k k,ddu,v H u,v/HM(u,v)其中，d、k均为小于1的常数。逆滤波的应用条件：退化图像g(x,y)是信噪比较高的图像。设：),(),(),(),(vuMvuNvuGvuFclose all;clear all;img=imread(noise.jpg);img_gray=rgb2gray(img);figure;imshow(img_gray);m,n=size(img_gray);f1,f2=freqspace(m,

23、 n,meshgrid);img_gray_f=fftshift(fft2(double(img_gray);figure;mesh(f1,f2,log10(1+abs(img_gray_f);psf = fspecial(disk,5); figure;freqz2(psf);Blurred = imfilter(img_gray,psf,circular,conv);figure; imshow(Blurred); title(Blurred Image);G=fftshift(fft2(double(Blurred);figure;mesh(f1,f2,log10(1+abs(G);H=

24、freqz2(psf,n,m);%r = sqrt(f1.2 + f2.2);d=0.4;k=1H(rd) =k;H(r=d) =1./H(r=d);figure;mesh(f1,f2,abs(H);K=G.*H;figure;mesh(f1,f2,log10(1+abs(K);F=ifft2(fftshift(K);figure;imshow(uint8(abs(F);title(Restore Image);逆滤波恢复演示对降质图像对降质图像g=Hf+ng=Hf+n而言，所谓恢复可以看成是对而言，所谓恢复可以看成是对g g通过某种线性通过某种线性变换变换L L而得到原图像的估计值而得到原图

25、像的估计值f f。其中，。其中，L L是是H H的某种逆过程的某种逆过程。如。如果果f f与与f f之间的差满足某些条件之间的差满足某些条件约束条件，例如：满足给定的约束条件，例如：满足给定的均方误差，则称该图像被恢复了。均方误差，则称该图像被恢复了。TnTTfTfTTTfTnTfTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTT2LLRLHLHR)fE(fLHR)fE(f 0)E( LHR)E(ffR)E(nn R)E(ff )LnnHnfHfnHL(HffffLnfLHffLgfff LfnLHffff ffffffff)f)(ff(f)LnH(fLgf Lgf )f)(f

26、ftrE(f)f)(ffEtr(fffE 阵图像和噪声的自相关矩 2所有的红圈值有令各分项分别为而噪声是独立的，噪声是独立的，与图像不相关。与图像不相关。求其期望应为求其期望应为0。期望值越小，其迹也越小期望值越小，其迹也越小对求迹而言，式中第2项和第3项的相等，有)LLRLHLHRLHRLHRtr(RTnTTffTTff2)LLRLHLHRLH2Rtr(RTnTTfTTff2问题转化为：L L等于什么时，均方误差最小。0nTfTf22LRH2LHRH2RL由此可得最佳L L：1nTfTf)RH(HRHRL0000将R Rf f，R Rn n看成是一个循环块矩阵，对其进行对角化，有：1111

27、WBWR WAWRWRWB WRWAnfnf图像的相关度约为图像的相关度约为2030个象素个象素合理的假设合理的假设1nTfTf)RH(HRHRL可得TfnTfHR)RHL(HR两边乘以WW和WW-1-1：1Tf1n1Tf1WHWR)WWRWH(WHRWLW1T1f1n1T1f11WWHWWR)WWRWWHWWR(WHWWLW*1T1DW WHDWHW又因为当D D为对角阵时有所以*1ADB)(DADWLW当L L为分块循环矩阵时，WL WWL W-1-1也是一对角阵。其对角线上元素的值为循环序列(“主值”序列)的傅立叶变换值。同理，D D，D D* *，A A，B B都为对角阵，其对角线上的

28、元素的值为对应循环序列(“主值”序列)的傅立叶变换值。于是上式可以写成标量形式),(),(),(),(),(),(2*vuRvuRvuHvuHvuRvuLnff合理的假设合理的假设f f(x,y)=g(x,y)(x,y)=g(x,y)* *l(x,y)l(x,y)基本原理 研究发现，逆滤波复原方法对噪声极为敏感，要求信噪比较高(100以上)，通常不满足该条件。 因此希望找到一种方法，在有噪声条件下，从退化图像g(x,y)复原出f(x,y)的估计值f(x,y) ，该估计值应符合一定的准则。 用向量f f, g g, n n来表示f(x,y), g(x, y), n(x,y)，Q Q为对f f的线

29、性算子，在约束条件 |g g-HfHf|2=|n n|2下求QfQf的最小化|QfQf|2而得到f f的最佳估计f f。拉格朗日乘数法拉格朗日乘数法是一种寻找变量受一个或多个条件所限制的多元是一种寻找变量受一个或多个条件所限制的多元函数的极值的方法。这种方法将一个有函数的极值的方法。这种方法将一个有n 个变量与个变量与k 个约束条件个约束条件的最优化问题转换为一个有的最优化问题转换为一个有n + k个变量的方程组的极值问题。个变量的方程组的极值问题。简单的例子，求方程的最大值：简单的例子，求方程的最大值：f(x,y) = x2y 同时未知数满足同时未知数满足x2 + y2 = 1 因为只有一个

30、未知数的限制条件，我们只需要用一个因为只有一个未知数的限制条件，我们只需要用一个乘数乘数.g(x,y) = x2 + y2 1 (x,y,) = f(x,y) + g(x,y) = x2y + (x2 + y2 1) 将所有将所有方程的偏微分设为零，得到一个方程组，最大值是以下方程的偏微分设为零，得到一个方程组，最大值是以下方程组的解中的一个：方程组的解中的一个：2xy + 2x = 0 x2 + 2y = 0 x2 + y2 1 = 0 222)(minnfHgfQfJgHQQHHfgHfHHQQ0fHgHfQQffTTTTTT1)( 1 )()(22)(s/sJT有令gHRRHHfTnfT

31、11)(s设R Rf f和R Rn n为f f和n n的自相关矩阵，根据前述分析结果：11 WBWRWAWRnnRffRnfTnTfEE若Q QT TQ Q用R Rf f1 1R Rn n来代替，有：拉格朗日系数因此：写成标量形式为：gWWDWBWWWAWDWWWDfWWDHWDWHT111111111)( *sgWDBADDfW1111)(*s),(),(/),(),(),(),(2*vuGvuRvuRsvuHvuHvuFfn波器退化为逆滤波器。时，没有噪声；维纳滤当称为参数维纳滤波器；时，称为维纳滤波器；时，0),(),(1),(1vuRvuFsvuFsn统计量的得到是非常困难的一件事。其

32、中其中R Rf f(u,v), R(u,v), Rn n(u,v)(u,v)分别是分别是f(x,y)f(x,y)和和n(x,y)n(x,y)的功率谱。的功率谱。),(),(),(),(1),(22vuGKvuHvuHvuHvuF关于图像信噪比的求法：关于图像信噪比的求法：第一步：将原图像第一步：将原图像f归一化到归一化到0-1，求归一化图像的各个像素点，求归一化图像的各个像素点的平方和作为信噪比中的信号量。的平方和作为信噪比中的信号量。第二步：第二步： 对归一化的原图像加噪对归一化的原图像加噪结果为归一化加噪图像结果为归一化加噪图像g第三步：用归一化原图像第三步：用归一化原图像f减去加噪图像减

33、去加噪图像g，并求结果图像各个，并求结果图像各个像素点的平方和作为信噪比中的噪声量像素点的平方和作为信噪比中的噪声量第四步：求出信噪比第四步：求出信噪比 MiNjMiNjjigjifjifSNR11211210),(),(),(log10注意：与峰值注意：与峰值信噪信噪比不同比不同close all;clear all;img=imread(gg.jpg);img_gray=rgb2gray(img);figure;imshow(img_gray);m,n=size(img_gray);psf = fspecial(average,25); Blurred = imfilter(img_gra

34、y,psf,circular,conv);figure; imshow(Blurred); title(Blurred Image);F=im2double(img_gray);S=sum(F(:).2); %估算信噪比估算信噪比G=F-im2double(Blurred);N=sum(G(:).2);SNR = 10*log10(S/N);parameter_s=0.001;%参数参数srestore=deconvwnr(Blurred,psf, parameter_s*(1/SNR);figure; imshow(restore); title(Restore Image);参数维纳滤波恢

35、复示例参数维纳滤波恢复示例close all;clear all;I= im2double(imread(2-man.jpg);m,n=size(I);figure;imshow(I); noise=0.1*randn(m,n); psf=fspecial(motion,21,11);blurred=imfilter(I,psf,circular);figure;imshow(blurred); blurrednoisy=im2uint8(blurred+noise); figure;imshow(blurrednoisy); NP=abs(fftn(noise).2; %功率谱与自相关函数是

36、一对傅立叶功率谱与自相关函数是一对傅立叶NCORR=fftshift(real(ifftn(NP); IP=abs(fftn(I).2; ICORR=fftshift(real(ifftn(IP); restore=deconvwnr(blurrednoisy,psf,NCORR,ICORR);figure;imshow(restore); 噪声情况已知时噪声情况已知时维纳滤波的使用维纳滤波的使用close all;clear all;img = imread(2-man.jpg);% img = imread(Vehicle_plate.jpg);figure;imshow(img); ti

37、tle(Restored); m,n=size(img);psf = fspecial(motion,100,40); Blurred = imfilter(img,psf,circular,conv); figure; imshow(Blurred);title(Blurred Image); restore = deconvwnr(Blurred,psf); %逆滤波逆滤波figure;imshow(restore); title(Restored); 图像运动模糊恢复实例图像运动模糊恢复实例八、最小平方(最大平滑)恢复与维纳滤波恢复的不同在于Q QT TQ Q的选择。维纳滤波中Q QT

38、TQ Q用R Rf f1 1R Rn n ，根据前面的推导，其依据的准则为均方误差最小。而最小平方恢复依据的准则为222minyfxf换言之：在约束条件|g g-HfHf|2=|n n|2条件下，恢复出来的图像的梯度最小(或称最大平滑)。根据此条件选择Q Q。如果f(x,y)尺寸是AB，则图像与p(x,y)相卷的结果后的图像为MN，MA+3-1， NB+3-1。构造p(x,y)、f(x,y)的扩展周期图像pe(x,y)与fe(x,y)做卷积。1313 02020 ),(),(11 01010 ),(),(N-y M-xy xyxpyxpN-y BM-xAB-y A-xyxfyxfee* * *

39、 * * * * * * * * * * * * * * * * * *0 1 01 -4 10 1 01010),(),(),(MmNneeenymxpnmfyxg),(4), 1() 1,() 1,(), 1(),(jifjifjifjifjifjig010141010),(yxp如果图像的g g、f f采用堆叠矢量的方法构成， 有g=Cfg=Cf) 0 ,() 2,() 1,() 2 ,() 0 ,() 1 ,() 1 ,() 1,() 0 ,( 021201110jpNjpNjpjpjpjpjpNjpjpeeeeeeeeejMMMCCCCCCCCCCCMN1 MNMN MN1=对C C

40、进行对角化，有 E=WE=W-1-1CW CW 在约束条件|g g-HfHf|2=|n n|2条件下，恢复出来的图像的梯度最小(或称最大平滑)。根据此条件选择Q Q。222)(minnfHgfQfJ左乘WW-1-1，有gWWD)EWWEDW(WDgHC)CHHf1*1*1*T1TT1(ssCWWEgWDE)ED(DfW-11*1 1其中s根据对角化的讨论，有),(),(),(),(),(22*vuGvuPsvuHvuHvuFgHQQHHfTT1)( s与维纳滤波不同此值为与维纳滤波不同此值为已知，不是统计量。已知，不是统计量。剩下的问题是对s的估测。构造残差矢量r r=g g-HfHf，调节s

41、使其满足|r r|2=|n n|2a(a是一个准确度系数)。赋给s某个初始值；计算f f和|r r|2。如果满足|r r|2=|n n|2a，停止计算；如果|r r|2|n n|2+a，减少s，继续计算。点对点计算点对点计算在不同条件下拍摄的图像，一个物体的图像常会发生几何失真，出现歪斜变形的现象。从太空中宇航器拍摄的地球上的等距平行线，其图像会变为歪斜或不等距；用光学和电子扫描仪摄取的图像常会有桶形失真和枕形失真；用普通的光学摄影与测试雷达拍摄的同一地区的景物二者在几何形状上有较大的差异。以一幅图像为基准，去校正另一种方式摄入的图像，以期校正其几何失真，就叫做图像的几何失真复原或者几何失真校

42、正。1、空间变换几何基准图像的坐标系统用f(x, y)来表示需要校正的图像的坐标系统用g(x, y)表示设两个图像坐标系统之间的关系用解析式表示 x=s(x, y) y=t(x,y)通常s(x,y)和t(x,y)用多项式来表示： y),( ),(y),( ),(2121121098765426524321321ykxykxkkxkkyxtkykxkyxtykxykxkkxkkyxskykxkyxs f (x, y) g(x, y) 标准图像标准图像扭曲图像扭曲图像通常用线性失真来近似较小的几何失真 x =a0+a1x+a2y y =b0+b1x+b2y对非线性失真更精确一些可以用二次型来近似

43、x =a0+a1x+a2y+a3x2+a4xy+a5y2 y =b0+b1x+b2y+b3x2+b4xy+b5y2若基准图像为f(x,y)，失真图像为g(x ,y )，对于景物上的同一个点，假定其灰度不变，则f(x, y)=g(x, y)，可利用已知点的对应点的坐标构造方程组求取ai、bj。几何失真复原的一套方法也可以用于使图像失真的工作中：在广告制作和计算机动画中常常要使物体变形。假设的变形关系假设的变形关系有时会遇到知道两点灰度值，要计算两点之间点的灰度值问题，如加大图像尺寸等。两点之间点的灰度值问题灰度差值问题。双线性差值法g(E)=(xE-xA)g(B)-g(A)+g(A)g(F)=

44、(xF-xC)g(D)-g(C)+g(C)g(H)=(yH-yA)g(C)-g(A)+g(A)g(I) = (yI-yB)g(D)-g(B)+g(B)x3377ABCDEFyIH0图像质量的优劣既可以通过人眼主观视觉效果来判断，也可以通过客观指标来衡量。1、均方误差(MSE)：2、峰值信噪比(PSNR)：其中M 、N分别是x方向、y方向图像像素点的个数，f(i,j)和f(i,j)分别是原始图像和测试图像在(i, j)点上的取值，L是图像中灰度取值的范围，对8比特的灰度图像而言，L=255。 MiNjjifjifMNMSE112),(),(1MSELPSNR210log10 恢复图像降质图像原始

45、图像),(),(),(),(),(),(),(log10ISNR,2,210jifjiyjifjifjifjiyjifjijiISNR只是评价图像恢复算法好坏的一个客观指标，ISNR高并不一定主观视觉效果好。 在我们观察的图象中，有一部分是通过数据人为构造的图象，如：CT图象、地址断面图象、安全检查设备输出的图象等。其中一类是由投影信息重构的图象本节的重点。1917年，奥地利数学家J.Radon证明了二维或三维物体可以从许多投影来重构其内部的数据。1963年，美国科学家A.M.Cormack首先将这一理论用于医学图象重构；1972年，英国科学家G.N.Homsfield设计出第一台X射线扫描仪

46、Computer TomgraphCT。二人于1979年双双获诺贝尔医学奖。平行光按传输方向投射到一个物体后，在垂直于平行光传输方向的平面上生成影象称为投影；用某种传感器取得的投影影象的数据称为投影数据。不同物质组合体的视图自然光投影视图X射线投影视图片状X射线投影视图最后一幅图象的数据就是重构断面图象的投影数据。显然，仅凭这组数据不能得到断面图象(里面物体的方园无法判断)，我们需要更多的数据。研究表明：当强度为I0的X射线通过吸收率为(x) 的物体时，有下面的关系：badxxdxxIIeIIba)(ln0)(0(x)xabI0III00计算机采样示意图基准检测器检测器物体+I0lnI0lnI

47、+-badxxII)(ln0A/Dx对应一个吸收率不均匀物质，(x,y)不为常数。可否根据得到的投影值求得(x,y)？结论：如果(x,y)已知，即可以重构图象。方法一：解联立方程法设每个网格中的吸收系数为x1 xN ，第i条射线与第xj个象素的相交长度为aij(可以计算得出)，代表第j个像素沿第i条射线的贡献的权值。像素编号每度采样一次射线条数jxapjjijki)( ,1790,1,k)(M ,1,2,i )(如果用pi(k)表示沿射线方向在射线角度为k时的总吸收测量值，则可通过解方程组的方法得到x1 xN 。长度长度aij可以计可以计算出来，只要算出来，只要找出找出63个独立个独立方程即可

48、。方程即可。1第第 i 条射线条射线第第j个像素个像素接收器接收器23N55被照射物质内的某一确切位置的吸收函数值是固定的。构造叠代公式进行叠代计算，其收敛值应是吸收函数值。NNNNkNkMRP (i,j) (i,j)1上次求得的上次求得的k(i,j)值值某角度为某角度为时的时的投影数据值投影数据值某角度为某角度为时时射线经过的像射线经过的像元数之和元数之和同一投影线路上同一投影线路上的上次的上次k(i,j)值值之和之和2.5+(3-5)/2=1.52.5+(3-5)/2=1.52.5+(7-5)/2=3.52.5+(7-5)/2=3.5给定叠代结束条件：112.53.5

49、3.752.7321.542.752.25143213725364154 (i,j) (i,j)NkNk1给定给定初值初值2.5xysstp(s,)f(x,y)X射线方向1)、构造两个坐标系，相交成角。根据解析几何的知识有tsyxyxtscossinsincoscossinsincos根据投影数据重构图像的理论进行图像重构。根据X射线与物体吸收率的关系有：dttsfspL),(),(3)、求投影p(s,)在某一时对s的傅立叶变换：)0 ,(),( )0(2exp),( 2exp),( 2exp),(),(SFTSFdtdsTsSjtsfdtdssSjt

50、sfdssSjspSPLL其中，T=0。可以写成F(S,T)=F(S,0)角度固定后，p(s,)数据的傅立叶变换s轴上各点数据的傅立叶变换。xysstp(s,)f(x,y)X射线方向L沿L方向的微分由于吸收值与座标系统无关(仅差坐标变换系数)，有：对f(x,y)和f(s,t)的傅立叶变换为：(1) ),()cossin(),sincos(),(tsftstsfyxf(3) )(2exp),(),(2) )(2exp),(),(dtdstTsSjtsfTSFdydxyYxXjyxfYXF 将(1)代入(2)有(雅可比行列式1，dxdy=dsdt)：dtdstYXsYXjtsfdtdsJYtsXt

51、sjtsfYXF )cossin()sincos(2exp),( )cossin()sincos(2exp),(),(1tytxsysxJ雅可比行列式雅可比行列式close all;clear all;img=imread(gg.jpg);img_gray=rgb2gray(img);figure;imshow(img_gray);m,n=size(img_gray);img_fft=fftshift(fft2(img_gray);mag_img=abs(img_fft)/m/n;figure;imshow(mag_img);r_img=imrotate(img_gray,45);figure

52、;imshow(r_img);m,n=size(r_img);imgr_fft=fftshift(fft2(r_img);mag_imgr=abs(imgr_fft)/m/n;figure;imshow(mag_imgr);傅立叶变换的旋转特性傅立叶变换的旋转特性这个变换是频域的一个旋转变换，于是有TSYXYXTScossinsincoscossinsincos0)0 ,(),(),( )3()(2exp),( )cossin()sincos(2exp),(),( TSFSPTSFdtdstTsSjtsfdtdstYXsYXjtsfYXF投影定理(以前图所示)：某角度投影数据对s的傅立叶变换s

53、轴上对应点的傅立叶变换。tsyxyxtscossinsincoscossinsincos对比如果已知F(X,Y)，则dYdXyYxXjYXFyxf )(2exp),(),(sincoscossinsincos0SSYXTSYXT该式即为重构图象的生成公式。SYXSYSXJ雅可比行列式雅可比行列式dSdSySxSjSPdYdXyYxXjSFdYdXyYxXjYXFyxf 0)sincos(2exp),( )(2exp)0 ,( )(2exp),(),(公式F(S,T)= P(S,)=F(S,0) T=0 告诉我们：通过采样数据得到的傅立叶变换值，实际上是x-y座标系统旋转角度后在s-t座标系统中s轴上数据的傅立叶变换，如下图所示。如果要得到真正图象的傅立叶变换，需对下图做适当的修改。STYX4321443322111111dddddPdPdPdPC式中：C为直角坐标定位点之变换值；P1，P2，P3， P4为与C点距离最近之四点之变换值； d1，d2，d3， d4为C点与P1，P2，P3，P4四点的距离。对于直角坐标定位点在旋转座标轴上之点，可采用零级内插公式得到：)(minidPCiCXY得到频域图象后，对其进行反变换，即可得到f(x,