第五章 图像复原与重建(2)

1、数字图像处理数字图像处理 主讲人：杜宏博第五章第五章 图像复原与重建图像复原与重建5.1 图像退化图像退化/复原处理的模型复原处理的模型5.2 噪声模型噪声模型5.3 空间滤波去噪空间滤波去噪5.4 频域滤波去噪频域滤波去噪5.5 退化函数建模退化函数建模5.6 图像复原的方法图像复原的方法 直接逆滤波直接逆滤波 维纳滤波维纳滤波5.7 图像投影重建图像投影重建周期噪声的模型是二维正弦波，通过带阻、带通和陷波滤波器可以被有效去除。理想带阻滤波器的表达式理想带阻滤波器的表达式: :00001,( , )2( , )0,( , )221,( , )2WD u vDWWH u vDD u vDWD

2、u vD5.4 频域滤波降低周期噪声频域滤波降低周期噪声 带阻滤波器带阻滤波器n n阶的巴特沃思带阻滤波器阶的巴特沃思带阻滤波器22201( , )( , )1( , )nH u vD u v WD u vD高斯带阻滤波器高斯带阻滤波器2220( , )12( , )( , )1Du vDD u v WH u ve 5.4 频域滤波降低周期噪声频域滤波降低周期噪声 带阻滤波器带阻滤波器理想带阻滤波器理想带阻滤波器巴特沃思带阻滤波器巴特沃思带阻滤波器高斯带阻滤波器高斯带阻滤波器5.4 频域滤波降低周期噪声频域滤波降低周期噪声带阻滤波器带阻滤波器(a) 被正弦噪声污染的图像被正弦噪声污染的图像 (

3、b) 图图(a)的频谱的频谱(c) 巴特沃思带阻滤波器巴特沃思带阻滤波器 (d) 滤波效果图滤波效果图5.4 频域滤波降低周期噪声频域滤波降低周期噪声( , )1( , )bpbrHu vHu v 带通滤波器带通滤波器带通滤波器执行与带阻滤波器相反的操作。带通滤波器执行与带阻滤波器相反的操作。( , )( , ):bpbrHu vHu v带通滤波器的传递函数可根据相应的带阻滤波器的传递函数得到不直接使用，损失大量不直接使用，损失大量图像细节。图像细节。可利用带通滤波器提取噪可利用带通滤波器提取噪声模式。声模式。5.4 频域滤波降低周期噪声频域滤波降低周期噪声陷波滤波器陷波滤波器阻止阻止(或通过

4、或通过)事先定义的中心频率邻域内的频率。事先定义的中心频率邻域内的频率。理想陷波滤波器理想陷波滤波器巴特沃思陷波滤波器巴特沃思陷波滤波器高斯陷波滤波器高斯陷波滤波器由于傅立叶变换由于傅立叶变换是对称的是对称的,因此因此陷波滤波器必须陷波滤波器必须以关于原点对称以关于原点对称的形式出现。的形式出现。5.4 频域滤波降低周期噪声频域滤波降低周期噪声陷波滤波器陷波滤波器00000,(,)(,)Du vuv半径为中心在且在对称的理想陷波滤波器的传递函数10200( , )( , )( , )1D u vDD u vDH u v或其他22 1/210022 1/2200( , )(/2)(/2) ( ,

5、 )(/2)(/2) D u vuMuvNvD u vuMuvNv其中5.4 频域滤波降低周期噪声频域滤波降低周期噪声 陷波滤波器陷波滤波器2012:1( , )1( , )/,nnH u vDD u vDu v阶数为 的巴特沃思陷波带阻滤波器的传递函数为1220( , )/,12:( , )1D u vDu vDH u ve 高斯陷波带阻滤波器的传递函数为 还可以得到另一种陷波滤波器还可以得到另一种陷波滤波器,它能通过它能通过(而不是阻止而不是阻止)包含在陷波区的频率包含在陷波区的频率.陷波区域的形状可以是任意的陷波区域的形状可以是任意的(如矩形如矩形)。5.4 频域滤波降低周期噪声频域滤波

6、降低周期噪声图像退化模型：图像退化模型：5.5 退化函数建模退化函数建模),(),(*),(),(yxyxfyxhyxg),(),(),(),(vuNvuFvuHvuG退化系统一般情况下是：线性，位置不变的退化系统退化系统一般情况下是：线性，位置不变的退化系统（1 1）线性：）线性：yxfbHyxfaHyxbfyxafH,2121（2 2）位置不变性：对任意）位置不变性：对任意, yxf有有yxgyxfH,对于线性位置不变退化，图像复原其实就是一个图像反卷积对于线性位置不变退化，图像复原其实就是一个图像反卷积的过程的过程 图像观察估计法图像观察估计法给定一幅退化图像，但没有退化函数 H 的知识

7、，那么估计该函数的方法之一就是收集图像自身的信息： 寻找简单结构的子图像 寻找受噪声影响小的子图像5.5 退化函数建模退化函数建模估计退化系统模型的三种方法估计退化系统模型的三种方法构造一个估计图像，它与观察的子图像有相同大小和特性 表示观察子图像， 表示构造的子图像 和 为对应的傅立叶变换。( , )( , )?/( , )SssH u vGu vF u v),(yxgs),(yxfs),(vuGs),(vuFs假设空间不变的，由 推导出完全函数 ),(vuHs),(vuH5.5 退化函数建模退化函数建模 图像试验估计法图像试验估计法 使用与被退化图像设备相似的装置，并得到一个脉冲的冲激响应

8、，可以进行较准确的退化估计：( , )( , )G u vH u vA一个脉冲点一个脉冲点成像系统成像系统H 此处A是一个冲激的傅立叶变换，表示冲击强度，为一常数。 右图为一个放大的亮脉冲以及退化的冲激。( , )g x y退化图像退化图像 模型估计法模型估计法 建立退化模型，考虑引起退化的环境因素。建立退化模型，考虑引起退化的环境因素。22 5/6()( , )k uvH u ve 例如：例如：Hufnagel Hufnagel 等等 Stanley Stanley 的退化模型就是基于大气湍的退化模型就是基于大气湍流的物理特性而提出来的，其中流的物理特性而提出来的，其中k k为常数，与湍流特

9、性相关。为常数，与湍流特性相关。( (除了指数除了指数5/65/6，该公式与高斯低通滤波形式相同，该公式与高斯低通滤波形式相同.).)5.5 退化函数建模退化函数建模 模型估计法模型估计法5.5 退化函数建模退化函数建模大气湍流模型模拟退化模糊一幅图像：大气湍流模型模拟退化模糊一幅图像：剧烈湍流剧烈湍流(k=0.0025)(k=0.0025)中等湍流中等湍流(k=0.001)(k=0.001)轻微湍流轻微湍流(k=0.00025)(k=0.00025)可忽略可忽略的湍流的湍流22 5/6()( , )k uvH u ve如果已知系统的传递函数 ，则根据 vuFvuHvuG,vuH,vuHvuG

10、vuF,可得复原图像的谱，经傅氏逆变换即可得到复原图像在忽略噪声的影响，退化模型的傅氏变换为实际应用时存在病态的问题，即在 H(u,v) 等于零或非常小的数值点上， 将变成无穷大或非常大的数。 ),(vuF- 这就是逆滤波复原法5.6 图像复原的方法图像复原的方法-逆滤波逆滤波 vuNvuFvuHvuG,vuHvuNvuFvuHvuNvuHvuGvuF,),(,系统中存在噪声时退化模型的傅立叶变换为：写成逆滤波复原的方式：1）即使知道退化函数，也不能准确复原图像，因为噪声函数 N(u,v) 是一个随机函数，其傅里叶变换未知。2）如果退化是零或非常小的值，噪声即使数值很小，但 N(u,v)/H(

11、u,v) 之比 (上式第二项) 可能非常大，很容易错估 的值。),(vuF12 ()( , )( , )( , )( , )jux vyf x yf x yN u v Hu vedudv 5.6 图像复原的方法图像复原的方法-逆滤波逆滤波解决退化是零或非常小的值的途径：限制滤波的频率，使其接近原点值。 在频率平面离原点较远的地方，H(u,v)数值较小或为零，因此图像复原在原点周围的有限区域内进行，即将退化图像的傅立叶谱限制在没出现零点而且数值又不是太小的有限范围内，即通过将频率限制为接近原点分析，减少了遇到零值的几率。 5.6 图像复原的方法图像复原的方法-逆滤波逆滤波剧烈湍流剧烈湍流( (k

12、 k=0.0025)=0.0025)大气湍流模型模拟退化模糊一幅图像大气湍流模型模拟退化模糊一幅图像可忽略的湍流可忽略的湍流652222/)/()/(),(NvMukevuH对退化函数对退化函数H H( (u u, ,v v) )进行精确取反并进行逆滤波，结果如下图。进行精确取反并进行逆滤波，结果如下图。5.6 图像复原的方法图像复原的方法-逆滤波逆滤波全频直接全频直接逆滤波复原逆滤波复原半径为半径为4040时截止时截止H H半径为半径为7070时截止时截止H H半径为半径为8585时截止时截止H H结果表明：噪声明显影响结果表明：噪声明显影响了图像复原结果，一般直了图像复原结果，一般直接逆滤

13、波效果较差。接逆滤波效果较差。剧烈湍流图剧烈湍流图( (k k=0.0025)=0.0025)5.6 图像复原的方法图像复原的方法-逆滤波逆滤波最小均方误差复原法最小均方误差复原法 -Wiener-Wiener滤波复原滤波复原目标：目标： 寻找一个滤波器，使得复原后图像寻找一个滤波器，使得复原后图像 与原与原始图像始图像 的均方误差最小。的均方误差最小。 逆滤波没有清楚说明如何处理噪声！逆滤波没有清楚说明如何处理噪声！),(yxf)(22ffEe误差度量：误差度量： 现讨论一种滤波复原法现讨论一种滤波复原法-Wiener-Wiener滤波复原：滤波复原： 综合考虑退化函数和噪声统计特征。综合考

14、虑退化函数和噪声统计特征。E期望值。期望值。),(yxf因此维纳滤波复原又称为最小均方误差复原。因此维纳滤波复原又称为最小均方误差复原。min),(),(2yxfyxfE5.6 图像复原的方法图像复原的方法-Wiener滤波复原滤波复原),(),(/ ),(),(),(),(),(vuGvuSvuSvuHvuHvuHvuFf221f 误差函数的最小值在频域里可以通过近似图像 的傅里叶变换来计算：叶变换为复原近似图像的傅里换为退化图像的傅里叶变换为退化函数的傅里叶变其中：),(),(),(vuFvuGvuH为未退化图像的功率谱为噪声的功率谱22),(),(),(),(vuFvuSvuNvuSf)

15、,(),(),(*vuHvuHvuH2维纳滤波器维纳滤波器5.6 图像复原的方法图像复原的方法-Wiener滤波复原滤波复原(2) (2) 未退化图像的功率谱难以知道，可用下式近似表示：未退化图像的功率谱难以知道，可用下式近似表示：KvuHvuHvuHvuHw221),(),(),(),(1) (1) 如果噪声为如果噪声为 0 0，其功率谱消失，维纳滤波就退化为逆滤波。，其功率谱消失，维纳滤波就退化为逆滤波。讨论：讨论：式中式中 K K 是根据信噪比的某种先验知识确定的常数。是根据信噪比的某种先验知识确定的常数。),(/ ),(),(),(),(),(vuSvuSvuHvuHvuHvuHfw2

16、21维纳滤波复原：维纳滤波复原：维纳滤波需要假定下述条件成立维纳滤波需要假定下述条件成立( (或近似成立或近似成立) )： 系统为线性、空间不变；系统为线性、空间不变；(1)(1) 退化图像、原始图像和噪声都是均匀随机场，噪声的均退化图像、原始图像和噪声都是均匀随机场，噪声的均值为零，且与图像不相关。值为零，且与图像不相关。5.6 图像复原的方法图像复原的方法-Wiener滤波复原滤波复原 维纳滤波复原与逆滤波复原的比较全频逆滤波半径受限逆滤波维纳滤波复原 (交互选择K) 维纳滤波的缺点： 未退化图像和噪声的功率谱必须是已知的； 功率比(信噪比)常数K 的估计一般还是没有合适的解。5.6 图像

17、复原的方法图像复原的方法-Wiener滤波复原滤波复原5.6 图像复原的方法图像复原的方法-Wiener滤波复原滤波复原维纳滤波器的维纳滤波器的matlab实现实现deconvwnr ：Deblur image using Wiener filter SyntaxJ = deconvwnr(I,PSF)J = deconvwnr(I,PSF,NSR)J = deconvwnr(I,PSF,NCORR,ICORR)其中：I是退化图像 PSF系统函数（点扩散函数） NSR信噪比 NCORR：噪声的自相关函数 ICORR：退化图像的自相关函数 J：复原图像5.6 图像复原的方法图像复原的方法-Wie

18、ner滤波复原滤波复原维纳滤波和逆滤波复原案例：维纳滤波和逆滤波复原案例：clc; clear; close clc; clear; close allall; ;f = double( imread(f = double( imread(cameraman.tifcameraman.tif););subplot(231); imshow(f,);subplot(231); imshow(f,);title(title(orginal clean imageorginal clean image); );% generate the degrade function% generate the

19、 degrade functionPSF = fspecial(PSF = fspecial(motionmotion,7,45);,7,45);subplot(232); imshow(PSF,);subplot(232); imshow(PSF,);title(title(Point spread functionPoint spread function); ); % using the PSF to degrade image% using the PSF to degrade imagegb = imfilter(f,PSF,gb = imfilter(f,PSF,circularc

20、ircular); );subplot(233); imshow(gb,);subplot(233); imshow(gb,);title(title(Blurred image caused by motionBlurred image caused by motion); ); 5.6 图像复原的方法图像复原的方法-Wiener滤波复原滤波复原维纳滤波和逆滤波复原案例：维纳滤波和逆滤波复原案例：% add noise to the degraded image% add noise to the degraded imagenoise = imnoise(zeros(size(f),noi

21、se = imnoise(zeros(size(f),gaussiangaussian,0,0.1);,0,0.1);g = gb + noise;g = gb + noise;subplot(234);imshow(g,)subplot(234);imshow(g,)title(title(Blurred image with noiseBlurred image with noise) ) % inverse filtering% inverse filteringfr1 = deconvwnr(g,PSF);fr1 = deconvwnr(g,PSF);subplot(235);imsh

22、ow(fr1,)subplot(235);imshow(fr1,)title(title(inverse filtering resultinverse filtering result) ) 5.6 图像复原的方法图像复原的方法-Wiener滤波复原滤波复原维纳滤波和逆滤波复原案例：维纳滤波和逆滤波复原案例：% wiener filtering% wiener filteringSn = abs(fft2(noise).2; % noise power spectrumnA = sum(Sn(:)/prod(size(noise); % noise average powerSf = abs

23、(fft2(f).2; % image power spectrumfA = sum(Sf(:)/prod(size(f); % image average powerR = nA/fA; % signal to noise ratiofr2 = deconvwnr(g,PSF,R);fr2 = deconvwnr(g,PSF,R);subplot(236);imshow(fr2,)subplot(236);imshow(fr2,)title(title(wiener filtering resultwiener filtering result) ) 5.6 图像复原的方法图像复原的方法-W

24、iener滤波复原滤波复原维纳滤波和逆滤波复原案例：维纳滤波和逆滤波复原案例：5.6 图像投影重建图像投影重建 概念：投影重建一般指利用物体的多个（轴向概念：投影重建一般指利用物体的多个（轴向）投影图像重建目标图像的过程。它是一类特）投影图像重建目标图像的过程。它是一类特殊的图像处理方法，输入的是一系列的投影图殊的图像处理方法，输入的是一系列的投影图，输出是重建图。，输出是重建图。 通过投影重建可以直接的看到原来被投影的物通过投影重建可以直接的看到原来被投影的物体的某种特性的空间分布，比直观观测投影图体的某种特性的空间分布，比直观观测投影图要直观的多。要直观的多。5.6 图像投影重建图像投影重

25、建 Radon变换对f(x,y)的Radon变换g(t, )定义为沿由t和定义的直线l的线积分。,( , ),( cossin)kkjg tf x ydlf x yxyt dxdy 5.6 图像投影重建图像投影重建 Radon变换(, ),( coscos)jjg tfx yxytdxdy Radon 变换揭示了函数和投影之间的关系，若函数为f (x, y)，则不同角度下的投影可写为原理：原理：“断层平面中某一点的密度值可看作这一平面内所有经过断层平面中某一点的密度值可看作这一平面内所有经过该点的射线投影之和（的平均值）该点的射线投影之和（的平均值）”5.6 图像投影重建图像投影重建 Rado

26、n变换的matlab实现radon：Radon transform SyntaxR = radon(I, theta)R,xp = radon(.)Description:R = radon(I, theta) returns the Radon transform R of the intensity image I for the angle theta degrees.The Radon transform is the projection of the image intensity along a radial line oriented at a specific angle.

27、If theta is a scalar, R is a column vector containing the Radon transform for theta degrees. If theta is a vector, R is a matrix in which each column is the Radon transform for one of the angles in theta. If you omit theta, it defaults to 0:179.R,xp = radon(.) returns a vector xp containing the radi

28、al coordinates corresponding to each row of R5.6 图像投影重建图像投影重建 Radon变换的matlab实现% generate two imagesg1 = zeros (600,600);g1(100:500,250:350)=1;g2 = phantom (Modified Shepp-Logan,600);subplot(221);imshow(g1,);sbplot(222);imshow(g2,) % radon transformtheta = 0:0.5:179.5;R1,xp1 = radon(g1,theta);R2,xp2

29、= radon(g2,theta);5.6 图像投影重建图像投影重建 Radon变换的matlab实现R1 = flipud(R1); % flip up and downR2 = flipud(R2);subplot(223);imshow(R1,XData,xp1(1 end),YData,179.5 0);axis xy;axis on;xlabel(rho);ylabel(theta);subplot(224);imshow(R2,XData,xp1(1 end),YData,179.5 0);axis xy; axis on;xlabel(rho); ylabel(theta);5.

30、6 图像投影重建图像投影重建 Radon变换的matlab实现5.6 图像投影重建图像投影重建 反投影重建法反投影重建法如何利用radon变换来重建图像f(x,y)？,( ,)( cossin,)kkkkkfx ygg xy (, ),( coscos)jjg tfx yxytdxdy 0,( , )f x yf x yd5.6 图像投影重建图像投影重建 反投影重建法反投影重建法第一步（first guess） =0 =2 =1 =31(0+1)5(2+3)1(0+1)5(2+3)把90角度的投影值加进空白图像实 例5.6 图像投影重建图像投影重建 反投影重建法反投影重建法第二步（second

31、 guess）021303331(0+1)5(2+3)1(0+1)5(2+3)+1(0+1)8(5+3)4(1+3)8(5+3)5.6 图像投影重建图像投影重建 反投影重建法反投影重建法第三步（third guess）021322443(2+1) 10(2+8)8(4+4) 12(4+8+1(0+1)8(5+3)4(1+3)8(5+35.6 图像投影重建图像投影重建 反投影重建法反投影重建法第四步（fourth guess）021332133(2+1) 10(2+8)8(4+4) 12(4+8+6(3+3)12(2+10)9(1+8)15(3+12）5.6 图像投影重建图像投影重建 反投影重建

32、法反投影重建法6129150213所有反投影的和0/36/33/39/36-6 12-69-6 15-6063902135.6 图像投影重建图像投影重建 反投影重建法反投影重建法缺点：星状伪影缺点：星状伪影000010000原始图像原始图像1/n1/n1/n1/n11/n1/n1/n1/n重构图像重构图像中心点中心点A经经n条投影线投影后，投条投影线投影后，投影值均为影值均为1：p1=p2=.=pn=1因此重建后因此重建后而其他点均为而其他点均为1/n：这类伪迹成为：这类伪迹成为星状伪影星状伪影121(.)1Anfpppn5.6 图像投影重建图像投影重建 反投影重建法反投影重建法缺点：星状伪影

33、缺点：星状伪影星状伪影星状伪影5.6 图像投影重建图像投影重建 反投影重建法反投影重建法缺点：图像模糊缺点：图像模糊图像产生模糊5.6 图像投影重建图像投影重建 傅里叶切片定理（中心切片定理）傅里叶切片定理（中心切片定理）物体空间f(x,y)Radon空间 g(R)傅立叶空间F(,)RR-1F1F2F2-1中心切片定理指出：中心切片定理指出：f(x,y)f(x,y)在某一方向上的投影函数在某一方向上的投影函数g g(R)(R)的一维傅立叶变换函数的一维傅立叶变换函数G G( () )是原函数是原函数f(x,y)f(x,y)的二的二维傅立叶变换函数维傅立叶变换函数F(F(, , ) )在在( (

34、, , ) )平面上沿同一方平面上沿同一方向且过原点的直线上的值。向且过原点的直线上的值。5.6 图像投影重建图像投影重建 傅里叶切片定理（中心切片定理）傅里叶切片定理（中心切片定理）投影函数的数学表达式：投影函数的数学表达式： dxdyRyxyxfdlyxfRg)sincos(),(),()(f(x,y)的二维傅立叶变换：的二维傅立叶变换：)()()sincos(),(),(),(122)sincos(2RgFdReRgdxdydReRyxyxfdxdyeyxfFRjRjyxj5.6 图像投影重建图像投影重建 傅里叶变换法傅里叶变换法2D IFT空间域空间域频域频域1D FT:(0 )( ,

35、 )f x y( )G插值插值( , )F u v( , )F ( )gR5.6 图像投影重建图像投影重建 滤波反投影法滤波反投影法dudeuFyxfyuxj 2,目标函数 f(x,y) 可由傅立叶函数F(u,v) 的逆变换获得，即5.6 图像投影重建图像投影重建 滤波反投影法滤波反投影法雅可比行列式sincosu频域中的笛卡尔坐标与极坐标的关系为：sincosudddduudud 5.6 图像投影重建图像投影重建 滤波反投影法滤波反投影法dudeuFyxfyuxj 2,deFddudeuFyxfyxjyuxjsincos20202sin,cos, dddudsincosu5.6 图像投影重建

36、图像投影重建 滤波反投影法滤波反投影法Let: F(cos, sin)=P(, ) dePddePddePddeFdyxfyxjyxjyxjyxjsincos200sincos200sincos2200sincos2200,sin,cos,P(w,) 为投影变换的一为投影变换的一维傅里叶变换维傅里叶变换5.6 图像投影重建图像投影重建 滤波反投影法滤波反投影法,tptp,PPdePddePddePdyxfyxjyxjyxjsincos20sincos200sincos200,5.6 图像投影重建图像投影重建 滤波反投影法滤波反投影法ddePyxftj20, P(, )表示对应于角度的单位投影的

37、傅立叶变换；里层的积分是P(, )|的逆傅立叶变换，记为g(t,)，在空间域，它表示单位投影被一频域响应为|的函数做滤波运算，故称之为滤波反投影1D Fourier transforminverse 1D Fourier transformbackprojection for all anglesfilter5.6 图像投影重建图像投影重建 滤波反投影法滤波反投影法1D FT空间域空间域频域频域1D IFT( , )f x y1( )F gR( )gR( )gR1( )F gR滤波器滤波器5.6 图像投影重建图像投影重建 滤波反投影法中滤波器的选择滤波反投影法中滤波器的选择频域中：滤波器|空间

38、域中与其对应的滤波器为： dettj2将t0代入上式计算得到(0)，即曲线|以下的面积。当时，(0)，所以上式是无法直接计算的，必须另想它法，引入限带函数（band-limiting function）5.6 图像投影重建图像投影重建 滤波反投影法中滤波器的选择滤波反投影法中滤波器的选择 2jtted滤波器是个无限频带的滤波函数, 由于其积分是发散的, 根据佩利- - 维纳准则, 这一理想滤波器是不可实现的。实际数值计算通常采用加窗的滤波函数。运用不同的窗函数可以得到不同的滤波器5.6 图像投影重建图像投影重建 滤波反投影法中滤波器的选择滤波反投影法中滤波器的选择 10H wothers 例如

39、频域滤波器例如频域滤波器5.6 图像投影重建图像投影重建 滤波反投影法中滤波器的选择滤波反投影法中滤波器的选择nRam-Lak: 矩形窗nShepp-Logan正弦窗nCosine: 余弦窗nHamming: 通用Hamming窗5.6 图像投影重建图像投影重建 滤波反投影法中滤波器的选择滤波反投影法中滤波器的选择(a)图为理想滤波器 (b)图为修正后滤波器 亦理论上滤波器亦称为Ramp滤波器，其高频分量是无限延伸的，但实际实现时必须截断处理，如图 (b)图中虚线所示，相当于在带宽之外突然衰减为零，在重建图像的边缘时会出现环状震荡条纹，称之为Gibbs现象。为有效地消除此现象，我们需对Ramp

40、滤波器稍作平滑处理，如将之与作卷积，得到Shepp-Logan滤波器； 5.6 图像投影重建图像投影重建 滤波反投影法中滤波器的选择滤波反投影法中滤波器的选择平滑了图像，损失了部分高频信息 Shepp-Logan滤波器5.6 图像投影重建图像投影重建 滤波反投影法中滤波器的选择滤波反投影法中滤波器的选择Hamming滤波器降低了高频噪声, 可得到Hamming滤波器和Hanning滤波器 5.6 图像投影重建图像投影重建反投影重建算法的反投影重建算法的matlab实现实现I = iradon(R, theta)I = iradon(P, theta, interp, filter, frequ

41、ency_scaling, output_size)I,H = iradon(.)DescriptionI = iradon(R, theta) reconstructs the image I from projection data in the two-dimensional array R. The columns of R are parallel beam projection data. iradon assumes that the center of rotation is the center point of the projections, which is defin

42、ed as ceil(size(R,1)/2).theta describes the angles (in degrees) at which the projections were taken. It can be either a vector containing the angles or a scalar specifying D_theta, the incremental angle between projections. If theta is a vector, it must contain angles with equal spacing between them

43、. If theta is a scalar specifying D_theta, the projections were taken at angles theta = m*D_theta, where m = 0,1,2,.,size(R,2)-1. If the input is the empty matrix (), D_theta defaults to 180/size(R,2).5.6 图像投影重建图像投影重建反投影重建算法的反投影重建算法的matlab实现实现I = iradon(P, theta, interp, filter, frequency_scaling, o

44、utput_size) specifies parameters to use in the inverse Radon transform. You can specify any combination of the last four arguments. iradon uses default values for any of these arguments that you erp specifies the type of interpolation to use in the back projection. The available options are

45、listed in order of increasing accuracy and computational complexity. nearest : Nearest-neighbor interpolation Linear: Linear interpolation (the default)spline: Spline interpolationcubic :Cubic interpolation from MATLAB 5. 5.6 图像投影重建图像投影重建反投影重建算法的反投影重建算法的matlab实现实现I = iradon(P, theta, interp, filter,

46、 frequency_scaling, output_size) specifies parameters to use in the inverse Radon transform. You can specify any combination of the last four arguments. iradon uses default values for any of these arguments that you erp specifies the type of interpolation to use in the back projection. The a

47、vailable options are listed in order of increasing accuracy and computational complexity. nearest : Nearest-neighbor interpolation Linear: Linear interpolation (the default)spline: Spline interpolationcubic :Cubic interpolation from MATLAB 5. 5.6 图像投影重建图像投影重建反投影重建算法的反投影重建算法的matlab实现实现I = iradon(P, t

48、heta, interp, filter, frequency_scaling, output_size) specifies parameters to use in the inverse Radon transform. You can specify any combination of the last four arguments. iradon uses default values for any of these arguments that you omit.filter specifies the filter to use for frequency domain fi

49、ltering. filter can be any of the strings that specify standard filters.Ram-LakShepp-LoganCosineHammingHannnone5.6 图像投影重建图像投影重建图像重建图像重建exampleclear;close clear;close allall; clc; clc; g = phantom (g = phantom (Modified Shepp-LoganModified Shepp-Logan,600);,600);subplot(241);imshow(g,);title (subplot

50、(241);imshow(g,);title (original imageoriginal image) )theta = 0:0.5:179.5;theta = 0:0.5:179.5;R,xp = radon(g,theta);R,xp = radon(g,theta);subplot(242);imshow(R,);title (subplot(242);imshow(R,);title (Projection imageProjection image) )f1 = iradon(R,theta,f1 = iradon(R,theta,nonenone); );subplot(243

51、);imshow(f1,);title (subplot(243);imshow(f1,);title (Reconstructed Image with back Reconstructed Image with back projectionprojection) )f2 = iradon(R,theta,f2 = iradon(R,theta,Ram-LakRam-Lak); );subplot(244);imshow(f2,);title (subplot(244);imshow(f2,);title (Reconstructed Image with Ram-Reconstructe

52、d Image with Ram-Lak filterLak filter) )5.6 图像投影重建图像投影重建图像重建图像重建examplef3 = iradon(R,theta,f3 = iradon(R,theta,Shepp-LoganShepp-Logan); );subplot(245);imshow(f3,);title (subplot(245);imshow(f3,);title (Reconstructed Image Reconstructed Image with Shepp-Logan filterwith Shepp-Logan filter) )f4 = irad

53、on(R,theta,f4 = iradon(R,theta,CosineCosine); );subplot(246);imshow(f4,);title (subplot(246);imshow(f4,);title (Reconstructed Image Reconstructed Image with Cosine filterwith Cosine filter) )f5 = iradon(R,theta,f5 = iradon(R,theta,HammingHamming); );subplot(247);imshow(f5,);title (subplot(247);imsho

54、w(f5,);title (Reconstructed Image Reconstructed Image with Hamming filterwith Hamming filter) )f6 = iradon(R,theta,f6 = iradon(R,theta,HannHann); );subplot(248);imshow(f6,);title (subplot(248);imshow(f6,);title (Reconstructed Image Reconstructed Image with Hann filterwith Hann filter) )5.6 图像投影重建图像投

55、影重建5.6 图像投影重建图像投影重建扇形扫描重建扇形扫描重建成像几何5.6 图像投影重建图像投影重建扇形扫描重建扇形扫描重建扇束情况下的重建算扇束情况下的重建算法较为复杂，但实质法较为复杂，但实质没有改变。可采用平没有改变。可采用平行束情况下的算法实行束情况下的算法实现，只需加以适当地现，只需加以适当地修正即可修正即可n 重排算法重排算法: : 把一个视图中采得的扇形数据重新组合成平行的把一个视图中采得的扇形数据重新组合成平行的 射线投影数据，然后采用平行束重建算法重建射线投影数据，然后采用平行束重建算法重建n 直接重建算法直接重建算法: : 不必数据重排，只需适当加权即可运用与平不必数据重

56、排，只需适当加权即可运用与平 行束类似的算法重建行束类似的算法重建5.6 图像投影重建图像投影重建扇形束重建- EA等角扇形束重建n当同样大小的探测器单元沿着中心为X射线焦点的弧排列时，就形成等角采样n扇形束的每一条射线可由和确定，其中是射线与中心射线（假想的通过X射线源和等中心的直线）的夹角，称为探测器角；是中心射线与y轴的夹角，称为投影角 5.6 图像投影重建图像投影重建扇形束重建- EA等角扇形束重建n 投影乘以探测器角的余弦，滤波后的样本随着到光源的距离的增长而增长n 重建公式可由用(t, ) 坐标确定(, ) 坐标上的每个样本来得到。扇形投影中的投影样本q(, )就转化为平行投影中的

57、投影样本p(t, ) 5.6 图像投影重建图像投影重建扇形投影和重建的扇形投影和重建的matlab实现实现fanbeam：Fan-beam transform Syntax:F = fanbeam(I,D)F = fanbeam(., param1, val1, param1, val2,.)F, fan_sensor_positions, fan_rotation_angles = fanbeam(.)FanRotationIncrement - Positive real scalar specifying the increment of the rotation angle of th

58、e fan-beam projections. Measured in degrees. Default value is 1.FanSensorGeometry - Text string specifying how sensors are positioned. Valid values are arc or line. In the arc geometry, sensors are spaced equally along a circular arc, This is the default value.FanSensorSpacing - Positive real scalar

59、 specifying the spacing of the fan-beam sensors. Interpretation of the value depends on the setting of FanSensorGeometry. If FanSensorGeometry is set to arc (the default), the value defines the angular spacing in degrees. Default value is 1. If FanSensorGeometry is line, the value specifies the line

60、ar spacing. Default value is 1.5.6 图像投影重建图像投影重建扇形投影和重建的扇形投影和重建的matlab实现实现fanbeam：Fan-beam transform Syntax:F = fanbeam(I,D)F = fanbeam(., param1, val1, param1, val2,.)F, fan_sensor_positions, fan_rotation_angles = fanbeam(.)FanRotationIncrement - Positive real scalar specifying the increment of the