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

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

2、 u vD WW H u vDD u vD W D u vD 5.4 频域滤波降低周期噪声频域滤波降低周期噪声 带阻滤波器带阻滤波器 n n阶的巴特沃思带阻滤波器阶的巴特沃思带阻滤波器 2 22 0 1 ( , ) ( , ) 1 ( , ) n H u v D u v W D u vD 高斯带阻滤波器高斯带阻滤波器 2 22 0 ( , )1 2( , ) ( , )1 Du vD D u v W H u ve 5.4 频域滤波降低周期噪声频域滤波降低周期噪声 带阻滤波器带阻滤波器 (a)理想带阻滤波器理想带阻滤波器 (b)巴特沃思带阻滤波器巴特沃思带阻滤波器 (c)高斯带阻滤波器高斯带阻滤

3、波器 5.4 频域滤波降低周期噪声频域滤波降低周期噪声 带阻滤波器带阻滤波器 (a) 被正弦噪声污染的图像被正弦噪声污染的图像 (b) 图图(a)的频谱的频谱 (c) 巴特沃思带阻滤波器巴特沃思带阻滤波器 (d) 滤波效果图滤波效果图 5.4 频域滤波降低周期噪声频域滤波降低周期噪声 ( , )1( , ) bpbr Hu vHu v 带通滤波器带通滤波器 带通滤波器执行与带阻滤波器相反的操作。带通滤波器执行与带阻滤波器相反的操作。 ( , ) ( , ): bp br Hu v Hu v 带通滤波器的传递函数可根据相应的 带阻滤波器的传递函数得到 不直接使用，损失大量不直接使用，损失大量 图

4、像细节。图像细节。 可利用带通滤波器提取噪可利用带通滤波器提取噪 声模式。声模式。 5.4 频域滤波降低周期噪声频域滤波降低周期噪声 陷波滤波器陷波滤波器 阻止阻止(或通过或通过)事先定义的中心频率邻域内的频率。事先定义的中心频率邻域内的频率。 (a) 理想陷波滤波器理想陷波滤波器 (b) 巴特沃思陷波滤波器巴特沃思陷波滤波器 (c) 高斯陷波滤波器高斯陷波滤波器 由于傅立叶变换由于傅立叶变换 是对称的是对称的,因此因此 陷波滤波器必须陷波滤波器必须 以关于原点对称以关于原点对称 的形式出现。的形式出现。 5.4 频域滤波降低周期噪声频域滤波降低周期噪声 陷波滤波器陷波滤波器 00000 ,(

5、,)(,)Du vuv半径为中心在且在对称的理想陷波滤波器的传递函数 1020 0( , )( , ) ( , ) 1 D u vDD u vD H u v 或 其他 22 1/2 100 22 1/2 200 ( , )(/2)(/2) ( , )(/2)(/2) D u vuMuvNv D u vuMuvNv 其中 5.4 频域滤波降低周期噪声频域滤波降低周期噪声 陷波滤波器陷波滤波器 2 0 12 : 1 ( , ) 1 ( , )/, n n H u v D D u vDu v 阶数为 的巴特沃思陷波带阻滤波器的传递函数为 12 2 0 ( , )/,1 2 : ( , )1 D u

6、vDu v D H u ve 高斯陷波带阻滤波器的传递函数为 还可以得到另一种陷波滤波器还可以得到另一种陷波滤波器,它能通过它能通过 (而不是阻止而不是阻止)包含在陷波区的频率包含在陷波区的频率. 陷波区域的形状可以是任意的陷波区域的形状可以是任意的(如矩形如矩形)。 5.4 频域滤波降低周期噪声频域滤波降低周期噪声 图像退化模型：图像退化模型： 5.5 退化函数建模退化函数建模 ),(),(*),(),(yxyxfyxhyxg ),(),(),(),(vuNvuFvuHvuG 退化系统一般情况下是：线性，位置不变的退化系统退化系统一般情况下是：线性，位置不变的退化系统 （1 1）线性：）线性

7、： yxfbHyxfaHyxbfyxafH, 2121 （2 2）位置不变性：对任意）位置不变性：对任意 , yxf 有有 yxgyxfH, 对于线性位置不变退化，图像复原其实就是一个图像反卷积对于线性位置不变退化，图像复原其实就是一个图像反卷积 的过程的过程 图像观察估计法图像观察估计法 给定一幅退化图像，但没有退化函数 H 的知识，那么估计 该函数的方法之一就是收集图像自身的信息： 寻找简单结构的子图像 寻找受噪声影响小的子图像 5.5 退化函数建模退化函数建模 估计退化系统模型的三种方法估计退化系统模型的三种方法 构造一个估计图像，它与观察的子图像有相同大小和特性 表示观察子图像， 表示

8、构造的子图像 和 为对应的傅立叶变换。 ( , )( , )?/( , ) S ss H u vGu vF u v ),(yxgs),( yxfs ),(vuGs),( vuFs 假设空间不变的，由 推导出完全函数 ),(vuH s ),(vuH 5.5 退化函数建模退化函数建模 图像试验估计法图像试验估计法 使用与被退化图像设备相似的装置，并得到一个脉冲的冲激 响应，可以进行较准确的退化估计： ( , ) ( , ) G u v H u v A 一个脉冲点一个脉冲点成像系统成像系统H 此处A是一个冲激的傅立叶变换，表示冲击强度，为一常数。 右图为一个放大的 亮脉冲以及退化的 冲激。 ( ,

9、)g x y 退化图像退化图像 模型估计法模型估计法 建立退化模型，考虑引起退化的环境因素。建立退化模型，考虑引起退化的环境因素。 22 5/6 () ( , ) k uv H u ve 例如：例如：Hufnagel Hufnagel 等等 Stanley Stanley 的退化模型就是基于大气湍的退化模型就是基于大气湍 流的物理特性而提出来的，其中流的物理特性而提出来的，其中k k为常数，与湍流特性相关。为常数，与湍流特性相关。 ( (除了指数除了指数5/65/6，该公式与高斯低通滤波形式相同，该公式与高斯低通滤波形式相同.).) 5.5 退化函数建模退化函数建模 模型估计法模型估计法 5.

10、5 退化函数建模退化函数建模 大气湍流模型模拟退化模糊一幅图像：大气湍流模型模拟退化模糊一幅图像： 剧烈湍流剧烈湍流 (k=0.0025)(k=0.0025) 中等湍流中等湍流 (k=0.001)(k=0.001) 轻微湍流轻微湍流 (k=0.00025)(k=0.00025) 可忽略可忽略 的湍流的湍流 22 5/6 () ( , ) k uv H u ve 如果已知系统的传递函数 ，则根据 vuFvuHvuG, vuH, vuHvuGvuF, 可得复原图像的谱，经傅氏逆变换即可得到复原图像 在忽略噪声的影响，退化模型的傅氏变换为 实际应用时存在病态的问题，即在 H(u,v) 等于零或非常小

11、 的数值点上， 将变成无穷大或非常大的数。 ),( vuF - 这就是逆滤波复原法 5.6 图像复原的方法图像复原的方法-逆滤波逆滤波 vuNvuFvuHvuG, vuH vuN vuF vuH vuN vuH vuG vuF , , ),( , , , , , 系统中存在噪声时退化模型的傅立叶变换为： 写成逆滤波复原的方式： 1）即使知道退化函数，也不能准确复原图像，因为噪声函数 N(u,v) 是 一个随机函数，其傅里叶变换未知。 2）如果退化是零或非常小的值，噪声即使数值很小，但 N(u,v)/H(u,v) 之比 (上式第二项) 可能非常大，很容易错估 的值。),( vuF 12 () (

12、 , )( , )( , )( , ) jux vy f x yf x yN u v Hu vedudv 5.6 图像复原的方法图像复原的方法-逆滤波逆滤波 解决退化是零或非常小的值的途径： 限制滤波的频率，使其接近原点值。 在频率平面离原点较远的地方，H(u,v)数值较小 或为零，因此图像复原在原点周围的有限区域内进 行，即将退化图像的傅立叶谱限制在没出现零点而 且数值又不是太小的有限范围内，即通过将频率限 制为接近原点分析，减少了遇到零值的几率。 5.6 图像复原的方法图像复原的方法-逆滤波逆滤波 剧烈湍流剧烈湍流( (k k=0.0025)=0.0025) 大气湍流模型模拟退化模糊一幅图

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

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

15、考虑退化函数和噪声统计特征。 E期望值。期望值。 ),(yxf 因此维纳滤波复原又称为最小均方误差复原。因此维纳滤波复原又称为最小均方误差复原。 min),(),( 2 yxfyxfE 5.6 图像复原的方法图像复原的方法-Wiener滤波复原滤波复原 ),( ),(/ ),(),( ),( ),( ),( vuG vuSvuSvuH vuH vuH vuF f 2 2 1 f 误差函数的最小值在频域里可以通过近似图像 的 傅里叶变换来计算： 叶变换为复原近似图像的傅里 换为退化图像的傅里叶变 换为退化函数的傅里叶变其中： ),( ),( ),( vuF vuG vuH 为未退化图像的功率谱

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

17、定的常数。是根据信噪比的某种先验知识确定的常数。 ),(/ ),(),( ),( ),( ),( vuSvuSvuH vuH vuH vuH f w 2 2 1 维纳滤波复原：维纳滤波复原： 维纳滤波需要假定下述条件成立维纳滤波需要假定下述条件成立( (或近似成立或近似成立) )： (1)(1) 系统为线性、空间不变；系统为线性、空间不变； (2)(2) 退化图像、原始图像和噪声都是均匀随机场，噪声的均退化图像、原始图像和噪声都是均匀随机场，噪声的均 值为零，且与图像不相关。值为零，且与图像不相关。 5.6 图像复原的方法图像复原的方法-Wiener滤波复原滤波复原 维纳滤波复原与逆滤波复原的

18、比较 全频逆滤波半径受限逆滤波维纳滤波复原 (交互选择K) 维纳滤波的缺点： 未退化图像和噪声的功率谱必须是已知的； 功率比(信噪比)常数K 的估计一般还是没有合适的解。 5.6 图像复原的方法图像复原的方法-Wiener滤波复原滤波复原 5.6 图像复原的方法图像复原的方法-Wiener滤波复原滤波复原 维纳滤波器的维纳滤波器的matlab实现实现 deconvwnr ：Deblur image using Wiener filter Syntax J = deconvwnr(I,PSF) J = deconvwnr(I,PSF,NSR) J = deconvwnr(I,PSF,NCORR,

19、ICORR) 其中：I是退化图像 PSF系统函数（点扩散函数） NSR信噪比 NCORR：噪声的自相关函数 ICORR：退化图像的自相关函数 J：复原图像 5.6 图像复原的方法图像复原的方法-Wiener滤波复原滤波复原 维纳滤波和逆滤波复原案例：维纳滤波和逆滤波复原案例： clc; clear; close clc; clear; close allall; ; f = double( imread(f = double( imread(cameraman.tifcameraman.tif);); subplot(231); imshow(f,);subplot(231); imshow(

20、f,); title(title(orginal clean imageorginal clean image); ); % generate the degrade function% generate the degrade function PSF = fspecial(PSF = fspecial(motionmotion,7,45);,7,45); subplot(232); imshow(PSF,);subplot(232); imshow(PSF,); title(title(Point spread functionPoint spread function); ); % us

21、ing the PSF to degrade image% using the PSF to degrade image gb = imfilter(f,PSF,gb = imfilter(f,PSF,circularcircular); ); subplot(233); imshow(gb,);subplot(233); imshow(gb,); title(title(Blurred image caused by motionBlurred image caused by motion); ); 5.6 图像复原的方法图像复原的方法-Wiener滤波复原滤波复原 维纳滤波和逆滤波复原案例

22、：维纳滤波和逆滤波复原案例： % add noise to the degraded image% add noise to the degraded image noise = imnoise(zeros(size(f),noise = 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 ima

23、ge with noise) ) % inverse filtering% inverse filtering fr1 = deconvwnr(g,PSF);fr1 = deconvwnr(g,PSF); subplot(235);imshow(fr1,)subplot(235);imshow(fr1,) title(title(inverse filtering resultinverse filtering result) ) 5.6 图像复原的方法图像复原的方法-Wiener滤波复原滤波复原 维纳滤波和逆滤波复原案例：维纳滤波和逆滤波复原案例： % wiener filtering% w

24、iener filtering Sn = abs(fft2(noise).2; % noise power spectrum nA = sum(Sn(:)/prod(size(noise); % noise average power Sf = abs(fft2(f).2; % image power spectrum fA = sum(Sf(:)/prod(size(f); % image average power R = nA/fA; % signal to noise ratio fr2 = deconvwnr(g,PSF,R);fr2 = deconvwnr(g,PSF,R); su

25、bplot(236);imshow(fr2,)subplot(236);imshow(fr2,) title(title(wiener filtering resultwiener filtering result) ) 5.6 图像复原的方法图像复原的方法-Wiener滤波复原滤波复原 维纳滤波和逆滤波复原案例：维纳滤波和逆滤波复原案例： 5.6 图像投影重建图像投影重建 概念：投影重建一般指利用物体的多个（轴向概念：投影重建一般指利用物体的多个（轴向 ）投影图像重建目标图像的过程。它是一类特）投影图像重建目标图像的过程。它是一类特 殊的图像处理方法，输入的是一系列的投影图殊的图像处理方法，

26、输入的是一系列的投影图 ，输出是重建图。，输出是重建图。 通过投影重建可以直接的看到原来被投影的物通过投影重建可以直接的看到原来被投影的物 体的某种特性的空间分布，比直观观测投影图体的某种特性的空间分布，比直观观测投影图 要直观的多。要直观的多。 5.6 图像投影重建图像投影重建 Radon变换 对f(x,y)的Radon变换g(t, )定义为沿由t和 定义的直线l的线积分。 ,( , ) ,( cossin) kkj g tf x ydl f x yxyt dxdy 5.6 图像投影重建图像投影重建 Radon变换 (, ),( coscos) jj g tfx yxytdxdy Radon

27、 变换揭示了函数和投影之间的关系，若函数为f (x, y)， 则不同角度下的投影可写为 原理：原理：“断层平面中某一点的密度值可看作这一平面内所有经过断层平面中某一点的密度值可看作这一平面内所有经过 该点的射线投影之和（的平均值）该点的射线投影之和（的平均值）” 5.6 图像投影重建图像投影重建 Radon变换的matlab实现 radon：Radon transform Syntax R = radon(I, theta) R,xp = radon(.) Description: R = radon(I, theta) returns the Radon transform R of the

28、 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. 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 eac

29、h 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 radial coordinates corresponding to each row of R 5.6 图像投影重建图像投影重建 Radon变换的matlab实现 % generate two images g1 = zeros (600,600); g1(100:500,250:3

30、50)=1; g2 = phantom (Modified Shepp-Logan,600); subplot(221);imshow(g1,); sbplot(222);imshow(g2,) % radon transform theta = 0:0.5:179.5; R1,xp1 = radon(g1,theta); R2,xp2 = radon(g2,theta); 5.6 图像投影重建图像投影重建 Radon变换的matlab实现 R1 = flipud(R1); % flip up and down R2 = flipud(R2); subplot(223); imshow(R1,

31、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.6 图像投影重建图像投影重建 Radon变换的matlab实现 5.6 图像投影重建图像投影重建 反投影重建法反投影重建法 如何利用radon变换来重建图像f(x,y)？ ,( ,)( cossin ,) k kkkk fx ygg

32、 xy (, ),( coscos) jj g tfx yxytdxdy 0 ,( , )f x yf x yd 5.6 图像投影重建图像投影重建 反投影重建法反投影重建法 第一步（first guess） =0 =2 =1 =3 1(0+1)5(2+3) 1(0+1)5(2+3) 把90角度的投影值加进空白图像 实 例 5.6 图像投影重建图像投影重建 反投影重建法反投影重建法 第二步（second guess） 02 13 03 33 1(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 图像投影重建图像投影重建 反投影重建法反

33、投影重建法 第三步（third guess） 02 13 22 44 3(2+1) 10(2+8) 8(4+4) 12(4+8 + 1(0+1)8(5+3) 4(1+3)8(5+3 5.6 图像投影重建图像投影重建 反投影重建法反投影重建法 第四步（fourth guess） 02 13 32 13 3(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 图像投影重建图像投影重建 反投影重建法反投影重建法 612 915 02 13 所有反投影的和 0/36/3 3/39/3 6-6 12-6 9-6 15-6 06 3

34、9 02 13 5.6 图像投影重建图像投影重建 反投影重建法反投影重建法 缺点：星状伪影缺点：星状伪影 000 010 000 原始图像原始图像 1/n1/n1/n 1/n11/n 1/n1/n1/n 重构图像重构图像 中心点中心点A经经n条投影线投影后，投条投影线投影后，投 影值均为影值均为1： p1=p2=.=pn=1 因此重建后因此重建后 而其他点均为而其他点均为1/n：这类伪迹成为：这类伪迹成为 星状伪影星状伪影 12 1 (.)1 An fppp n 5.6 图像投影重建图像投影重建 反投影重建法反投影重建法 缺点：星状伪影缺点：星状伪影 星状伪影星状伪影 5.6 图像投影重建图像

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

36、上沿同一方平面上沿同一方 向且过原点的直线上的值。向且过原点的直线上的值。 5.6 图像投影重建图像投影重建 傅里叶切片定理（中心切片定理）傅里叶切片定理（中心切片定理） 投影函数的数学表达式：投影函数的数学表达式： dxdyRyxyxf dlyxfRg )sincos(),( ),()( f(x,y)的二维傅立叶变换：的二维傅立叶变换： )( )( )sincos(),( ),(),( 1 2 2 )sincos(2 RgF dReRg dxdydReRyxyxf dxdyeyxfF Rj Rj yxj 5.6 图像投影重建图像投影重建 傅里叶变换法傅里叶变换法 2D IFT 空间域空间域

37、频域频域 1D FT :(0 ) ( , )f x y ( )G 插值插值 ( , )F u v ( , )F ( )gR 5.6 图像投影重建图像投影重建 滤波反投影法滤波反投影法 dudeuFyxf yuxj 2 , 目标函数 f(x,y) 可由傅立叶函数F(u,v) 的逆 变换获得，即 5.6 图像投影重建图像投影重建 滤波反投影法滤波反投影法 雅可比行列式 sin cos u 频域中的笛卡尔坐标与 极坐标的关系为： sin cos u dddd uu dud 5.6 图像投影重建图像投影重建 滤波反投影法滤波反投影法 dudeuFyxf yuxj 2 , deFd dudeuFyxf

38、yxj yuxj sincos2 0 2 0 2 sin,cos , dddud sin cos u 5.6 图像投影重建图像投影重建 滤波反投影法滤波反投影法 Let: F(cos, sin)=P(, ) dePd dePd dePd deFdyxf yxj yxj yxj yxj sincos2 00 sincos2 00 sincos2 2 00 sincos2 2 00 , , , sin,cos, P(w,) 为投影变换的一为投影变换的一 维傅里叶变换维傅里叶变换 5.6 图像投影重建图像投影重建 滤波反投影法滤波反投影法 ,tptp ,PP dePd dePd dePdyxf yx

39、j yxj yxj sincos2 0 sincos2 00 sincos2 00 , , , 5.6 图像投影重建图像投影重建 滤波反投影法滤波反投影法 ddePyxf tj2 0 , P(, )表示对应于角度的单位投影的傅立叶变换；里层的积分是 P(, )|的逆傅立叶变换，记为g(t,)，在空间域，它表示单位 投影被一频域响应为|的函数做滤波运算，故称之为滤波反投 影 1D Fourier transform inverse 1D Fourier transform backprojection for all angles filter 5.6 图像投影重建图像投影重建 滤波反投影法滤波

40、反投影法 1D FT 空间域空间域 频域频域 1D IFT ( , )f x y 1 ( )F gR ( )gR ( )gR 1 ( )F gR 滤波器滤波器 5.6 图像投影重建图像投影重建 滤波反投影法中滤波器的选择滤波反投影法中滤波器的选择 频域中：滤波器| 空间域中与其对应的滤波器为： det tj 2 将t0代入上式计算得到(0)，即曲线|以下的面积。 当时，(0)，所以上式是无法直接计算的，必须 另想它法，引入限带函数（band-limiting function） 5.6 图像投影重建图像投影重建 滤波反投影法中滤波器的选择滤波反投影法中滤波器的选择 2jt ted 滤波器是个无

41、限频带的滤波函数, 由于其积分是发散的, 根 据佩利- - 维纳准则, 这一理想滤波器是不可实现的。实际 数值计算通常采用加窗的滤波函数。运用不同的窗函数可以 得到不同的滤波器 5.6 图像投影重建图像投影重建 滤波反投影法中滤波器的选择滤波反投影法中滤波器的选择 1 0 H w others 例如频域滤波器例如频域滤波器 5.6 图像投影重建图像投影重建 滤波反投影法中滤波器的选择滤波反投影法中滤波器的选择 nRam-Lak: 矩形窗 nShepp-Logan正弦窗 nCosine: 余弦窗 nHamming: 通用Hamming窗 5.6 图像投影重建图像投影重建 滤波反投影法中滤波器的选

42、择滤波反投影法中滤波器的选择 (a)图为理想滤波器 (b)图为修正后滤波器 亦 理论上滤波器亦称为Ramp滤波器，其高频分量是无限延伸的，但 实际实现时必须截断处理，如图 (b)图中虚线所示，相当于在带 宽之外突然衰减为零，在重建图像的边缘时会出现环状震荡条纹， 称之为Gibbs现象。为有效地消除此现象，我们需对Ramp滤波器稍 作平滑处理，如将之与作卷积，得到Shepp-Logan滤波器； 5.6 图像投影重建图像投影重建 滤波反投影法中滤波器的选择滤波反投影法中滤波器的选择 平滑了图像，损失了部分高频信息 Shepp-Logan滤波器 5.6 图像投影重建图像投影重建 滤波反投影法中滤波器

43、的选择滤波反投影法中滤波器的选择 Hamming滤波器 降低了高频噪声, 可得到Hamming滤波 器和Hanning滤波器 5.6 图像投影重建图像投影重建 反投影重建算法的反投影重建算法的matlab实现实现 I = iradon(R, theta) I = iradon(P, theta, interp, filter, frequency_scaling, output_size) I,H = iradon(.) Description I = iradon(R, theta) reconstructs the image I from projection data in the t

44、wo- 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 defined as ceil(size(R,1)/2). theta describes the angles (in degrees) at which the projections were taken. It can be either a vecto

45、r 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. 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

46、)- 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, output_size) specifies parameters to use in the inverse Radon transform. You can specify any combination of the last four

47、arguments. iradon uses default values for any of these arguments that you omit. interp specifies the type of interpolation to use in the back projection. The available options are listed in order of increasing accuracy and computational complexity. nearest : Nearest-neighbor interpolation Linear: Li

48、near interpolation (the default) spline: Spline interpolation cubic :Cubic interpolation from MATLAB 5. 5.6 图像投影重建图像投影重建 反投影重建算法的反投影重建算法的matlab实现实现 I = iradon(P, theta, interp, filter, frequency_scaling, output_size) specifies parameters to use in the inverse Radon transform. You can specify any com

49、bination of the last four arguments. iradon uses default values for any of these arguments that you omit. interp specifies the type of interpolation to use in the back projection. The available options are listed in order of increasing accuracy and computational complexity. nearest : Nearest-neighbo

50、r interpolation Linear: Linear interpolation (the default) spline: Spline interpolation cubic :Cubic interpolation from MATLAB 5. 5.6 图像投影重建图像投影重建 反投影重建算法的反投影重建算法的matlab实现实现 I = iradon(P, theta, interp, filter, frequency_scaling, output_size) specifies parameters to use in the inverse Radon transfor

51、m. 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 filtering. filter can be any of the strings that specify standard filters. Ram-Lak Shepp-Logan Cosine Hamming Ha

52、nn none 5.6 图像投影重建图像投影重建 图像重建图像重建example clear;close clear;close allall; clc; clc; g = phantom (g = phantom (Modified Shepp-LoganModified Shepp-Logan,600);,600); subplot(241);imshow(g,);title (subplot(241);imshow(g,);title (original imageoriginal image) ) theta = 0:0.5:179.5;theta = 0:0.5:179.5; R,x

53、p = 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);imshow(f1,);title (subplot(243);imshow(f1,);title (Reconstructed Image with back Reconstructe

54、d 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-Reconstructed Image with Ram- Lak filterLak filter) ) 5.6 图像投影重建图像投影重建 图像重建图像重建example f3 = iradon(R,thet

55、a,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 = iradon(R,theta,f4 = iradon(R,theta,CosineCosine); ); subplot(246);imshow(f4,);title (subplo

56、t(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);imshow(f5,);title (Reconstructed Image Reconstructed Image with Hamming filterwith Hammin

57、g 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 图像投影重建图像投影重建 5.6 图像投影重建图像投影重建 扇形扫描重建扇形扫描重建 成像几何 5.6 图像投影重建图像投影重建 扇形扫描重建扇形扫描重建 扇束情况下的重建算扇束情

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

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

60、b实现实现 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 the fan-beam projections. Measured in degrees