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

数字图像处理,主讲人：杜宏博,第五章图像复原与重建,5.1图像退化/复原处理的模型5.2噪声模型5.3空间滤波去噪5.4频域滤波去噪5.5退化函数建模5.6图像复原的方法直接逆滤波维纳滤波5.7图像投影重建,周期噪声的模型是二维正弦波，通过带阻、带通和陷波滤波器可以被有效去除。,理想带阻滤波器的表达式:,5.4频域滤波降低周期噪声,带阻滤波器,n阶的巴特沃思带阻滤波器,高斯带阻滤波器,5.4频域滤波降低周期噪声,带阻滤波器,理想带阻滤波器巴特沃思带阻滤波器高斯带阻滤波器,5.4频域滤波降低周期噪声,带阻滤波器,(a)被正弦噪声污染的图像(b)图(a)的频谱(c)巴特沃思带阻滤波器(d)滤波效果图,5.4频域滤波降低周期噪声,带通滤波器,带通滤波器执行与带阻滤波器相反的操作。,不直接使用，损失大量图像细节。,可利用带通滤波器提取噪声模式。,5.4频域滤波降低周期噪声,陷波滤波器,阻止(或通过)事先定义的中心频率邻域内的频率。,理想陷波滤波器巴特沃思陷波滤波器高斯陷波滤波器,由于傅立叶变换是对称的,因此陷波滤波器必须以关于原点对称的形式出现。,5.4频域滤波降低周期噪声,陷波滤波器,5.4频域滤波降低周期噪声,陷波滤波器,还可以得到另一种陷波滤波器,它能通过(而不是阻止)包含在陷波区的频率.,陷波区域的形状可以是任意的(如矩形)。,5.4频域滤波降低周期噪声,图像退化模型：,5.5退化函数建模,退化系统一般情况下是：线性，位置不变的退化系统,（1）线性：,（2）位置不变性：对任意,有,对于线性位置不变退化，图像复原其实就是一个图像反卷积的过程,图像观察估计法,给定一幅退化图像，但没有退化函数H的知识，那么估计该函数的方法之一就是收集图像自身的信息：寻找简单结构的子图像寻找受噪声影响小的子图像,5.5退化函数建模,估计退化系统模型的三种方法,构造一个估计图像，它与观察的子图像有相同大小和特性表示观察子图像，表示构造的子图像和为对应的傅立叶变换。,假设空间不变的，由推导出完全函数,5.5退化函数建模,图像试验估计法,使用与被退化图像设备相似的装置，并得到一个脉冲的冲激响应，可以进行较准确的退化估计：,一个脉冲点,成像系统H,此处A是一个冲激的傅立叶变换，表示冲击强度，为一常数。,右图为一个放大的亮脉冲以及退化的冲激。,退化图像,模型估计法,建立退化模型，考虑引起退化的环境因素。,例如：Hufnagel等Stanley的退化模型就是基于大气湍流的物理特性而提出来的，其中k为常数，与湍流特性相关。,(除了指数5/6，该公式与高斯低通滤波形式相同.),5.5退化函数建模,模型估计法,5.5退化函数建模,大气湍流模型模拟退化模糊一幅图像：,剧烈湍流(k=0.0025),中等湍流(k=0.001),轻微湍流(k=0.00025),可忽略的湍流,如果已知系统的传递函数，则根据,可得复原图像的谱，经傅氏逆变换即可得到复原图像,在忽略噪声的影响，退化模型的傅氏变换为,实际应用时存在病态的问题，即在H(u,v)等于零或非常小的数值点上，将变成无穷大或非常大的数。,-这就是逆滤波复原法,5.6图像复原的方法-逆滤波,系统中存在噪声时退化模型的傅立叶变换为：,写成逆滤波复原的方式：,1）即使知道退化函数，也不能准确复原图像，因为噪声函数N(u,v)是一个随机函数，其傅里叶变换未知。,2）如果退化是零或非常小的值，噪声即使数值很小，但N(u,v)/H(u,v)之比(上式第二项)可能非常大，很容易错估的值。,5.6图像复原的方法-逆滤波,解决退化是零或非常小的值的途径：限制滤波的频率，使其接近原点值。,在频率平面离原点较远的地方，H(u,v)数值较小或为零，因此图像复原在原点周围的有限区域内进行，即将退化图像的傅立叶谱限制在没出现零点而且数值又不是太小的有限范围内，即通过将频率限制为接近原点分析，减少了遇到零值的几率。,5.6图像复原的方法-逆滤波,剧烈湍流(k=0.0025),大气湍流模型模拟退化模糊一幅图像,可忽略的湍流,对退化函数H(u,v)进行精确取反并进行逆滤波，结果如下图。,5.6图像复原的方法-逆滤波,全频直接逆滤波复原,半径为40时截止H,半径为70时截止H,半径为85时截止H,结果表明：噪声明显影响了图像复原结果，一般直接逆滤波效果较差。,剧烈湍流图(k=0.0025),5.6图像复原的方法-逆滤波,最小均方误差复原法-Wiener滤波复原,目标：寻找一个滤波器，使得复原后图像与原始图像的均方误差最小。,逆滤波没有清楚说明如何处理噪声！,误差度量：,现讨论一种滤波复原法-Wiener滤波复原：综合考虑退化函数和噪声统计特征。,期望值。,因此维纳滤波复原又称为最小均方误差复原。,5.6图像复原的方法-Wiener滤波复原,误差函数的最小值在频域里可以通过近似图像的傅里叶变换来计算：,维纳滤波器,5.6图像复原的方法-Wiener滤波复原,(2)未退化图像的功率谱难以知道，可用下式近似表示：,(1)如果噪声为0，其功率谱消失，维纳滤波就退化为逆滤波。,讨论：,式中K是根据信噪比的某种先验知识确定的常数。,维纳滤波复原：,维纳滤波需要假定下述条件成立(或近似成立)：系统为线性、空间不变；退化图像、原始图像和噪声都是均匀随机场，噪声的均值为零，且与图像不相关。,5.6图像复原的方法-Wiener滤波复原,维纳滤波复原与逆滤波复原的比较,全频逆滤波,半径受限逆滤波,维纳滤波复原(交互选择K),维纳滤波的缺点：未退化图像和噪声的功率谱必须是已知的；功率比(信噪比)常数K的估计一般还是没有合适的解。,5.6图像复原的方法-Wiener滤波复原,5.6图像复原的方法-Wiener滤波复原,维纳滤波器的matlab实现,deconvwnr：DeblurimageusingWienerfilterSyntaxJ=deconvwnr(I,PSF)J=deconvwnr(I,PSF,NSR)J=deconvwnr(I,PSF,NCORR,ICORR),其中：I是退化图像PSF系统函数（点扩散函数）NSR信噪比NCORR：噪声的自相关函数ICORR：退化图像的自相关函数J：复原图像,5.6图像复原的方法-Wiener滤波复原,维纳滤波和逆滤波复原案例：,clc;clear;closeall;f=double(imread(cameraman.tif);subplot(231);imshow(f,);title(orginalcleanimage);%generatethedegradefunctionPSF=fspecial(motion,7,45);subplot(232);imshow(PSF,);title(Pointspreadfunction);%usingthePSFtodegradeimagegb=imfilter(f,PSF,circular);subplot(233);imshow(gb,);title(Blurredimagecausedbymotion);,5.6图像复原的方法-Wiener滤波复原,维纳滤波和逆滤波复原案例：,%addnoisetothedegradedimagenoise=imnoise(zeros(size(f),gaussian,0,0.1);g=gb+noise;subplot(234);imshow(g,)title(Blurredimagewithnoise)%inversefilteringfr1=deconvwnr(g,PSF);subplot(235);imshow(fr1,)title(inversefilteringresult),5.6图像复原的方法-Wiener滤波复原,维纳滤波和逆滤波复原案例：,%wienerfilteringSn=abs(fft2(noise).2;%noisepowerspectrumnA=sum(Sn(:)/prod(size(noise);%noiseaveragepowerSf=abs(fft2(f).2;%imagepowerspectrumfA=sum(Sf(:)/prod(size(f);%imageaveragepowerR=nA/fA;%signaltonoiseratiofr2=deconvwnr(g,PSF,R);subplot(236);imshow(fr2,)title(wienerfilteringresult),5.6图像复原的方法-Wiener滤波复原,维纳滤波和逆滤波复原案例：,5.6图像投影重建,概念：投影重建一般指利用物体的多个（轴向）投影图像重建目标图像的过程。它是一类特殊的图像处理方法，输入的是一系列的投影图，输出是重建图。通过投影重建可以直接的看到原来被投影的物体的某种特性的空间分布，比直观观测投影图要直观的多。,5.6图像投影重建,Radon变换,对f(x,y)的Radon变换g(t,)定义为沿由t和定义的直线l的线积分。,5.6图像投影重建,Radon变换,Radon变换揭示了函数和投影之间的关系，若函数为f(x,y)，则不同角度下的投影可写为,原理：“断层平面中某一点的密度值可看作这一平面内所有经过该点的射线投影之和（的平均值）”,5.6图像投影重建,Radon变换的matlab实现,radon：RadontransformSyntaxR=radon(I,theta)R,xp=radon(.)Description:R=radon(I,theta)returnstheRadontransformRoftheintensityimageIfortheanglethetadegrees.TheRadontransformistheprojectionoftheimageintensityalongaradiallineorientedataspecificangle.Ifthetaisascalar,RisacolumnvectorcontainingtheRadontransformforthetadegrees.Ifthetaisavector,RisamatrixinwhicheachcolumnistheRadontransformforoneoftheanglesintheta.Ifyouomittheta,itdefaultsto0:179.R,xp=radon(.)returnsavectorxpcontainingtheradialcoordinatescorrespondingtoeachrowofR,5.6图像投影重建,Radon变换的matlab实现,%generatetwoimagesg1=zeros(600,600);g1(100:500,250:350)=1;g2=phantom(ModifiedShepp-Logan,600);subplot(221);imshow(g1,);sbplot(222);imshow(g2,)%radontransformtheta=0:0.5:179.5;R1,xp1=radon(g1,theta);R2,xp2=radon(g2,theta);,5.6图像投影重建,Radon变换的matlab实现,R1=flipud(R1);%flipupanddownR2=flipud(R2);subplot(223);imshow(R1,XData,xp1(1end),YData,179.50);axisxy;axison;xlabel(rho);ylabel(theta);subplot(224);imshow(R2,XData,xp1(1end),YData,179.50);axisxy;axison;xlabel(rho);ylabel(theta);,5.6图像投影重建,Radon变换的matlab实现,5.6图像投影重建,反投影重建法,如何利用radon变换来重建图像f(x,y)？,5.6图像投影重建,反投影重建法,第一步（firstguess）,把90角度的投影值加进空白图像,实例,5.6图像投影重建,反投影重建法,第二步（secondguess）,+,5.6图像投影重建,反投影重建法,第三步（thirdguess）,+,5.6图像投影重建,反投影重建法,第四步（fourthguess）,+,5.6图像投影重建,反投影重建法,所有反投影的和,5.6图像投影重建,反投影重建法,缺点：星状伪影,原始图像,重构图像,中心点A经n条投影线投影后，投影值均为1：p1=p2=.=pn=1因此重建后而其他点均为1/n：这类伪迹成为星状伪影,5.6图像投影重建,反投影重建法,缺点：星状伪影,星状伪影,5.6图像投影重建,反投影重建法,缺点：图像模糊,图像产生模糊,5.6图像投影重建,傅里叶切片定理（中心切片定理）,物体空间f(x,y),Radon空间g(R),傅立叶空间F(,),中心切片定理指出：f(x,y)在某一方向上的投影函数g(R)的一维傅立叶变换函数G()是原函数f(x,y)的二维傅立叶变换函数F(,)在(,)平面上沿同一方向且过原点的直线上的值。,5.6图像投影重建,傅里叶切片定理（中心切片定理）,投影函数的数学表达式：,f(x,y)的二维傅立叶变换：,5.6图像投影重建,傅里叶变换法,5.6图像投影重建,滤波反投影法,目标函数f(x,y)可由傅立叶函数F(u,v)的逆变换获得，即,5.6图像投影重建,滤波反投影法,雅可比行列式,频域中的笛卡尔坐标与极坐标的关系为：,5.6图像投影重建,滤波反投影法,5.6图像投影重建,滤波反投影法,Let:F(cos,sin)=P(,),P(w,)为投影变换的一维傅里叶变换,5.6图像投影重建,滤波反投影法,5.6图像投影重建,滤波反投影法,P(,)表示对应于角度的单位投影的傅立叶变换；里层的积分是P(,)|的逆傅立叶变换，记为g(t,)，在空间域，它表示单位投影被一频域响应为|的函数做滤波运算，故称之为滤波反投影,1DFouriertransform,inverse1DFouriertransform,backprojectionforallangles,filter,5.6图像投影重建,滤波反投影法,滤波器,5.6图像投影重建,滤波反投影法中滤波器的选择,频域中：滤波器|空间域中与其对应的滤波器为：,将t0代入上式计算得到(0)，即曲线|以下的面积。当时，(0)，所以上式是无法直接计算的，必须另想它法，引入限带函数（band-limitingfunction）,5.6图像投影重建,滤波反投影法中滤波器的选择,滤波器是个无限频带的滤波函数,由于其积分是发散的,根据佩利-维纳准则,这一理想滤波器是不可实现的。实际数值计算通常采用加窗的滤波函数。运用不同的窗函数可以得到不同的滤波器,5.6图像投影重建,滤波反投影法中滤波器的选择,例如频域滤波器,5.6图像投影重建,滤波反投影法中滤波器的选择,Ram-Lak:矩形窗Shepp-Logan正弦窗Cosine:余弦窗Hamming:通用Hamming窗,5.6图像投影重建,滤波反投影法中滤波器的选择,(a)图为理想滤波器(b)图为修正后滤波器,亦,理论上滤波器亦称为Ramp滤波器，其高频分量是无限延伸的，但实际实现时必须截断处理，如图(b)图中虚线所示，相当于在带宽之外突然衰减为零，在重建图像的边缘时会出现环状震荡条纹，称之为Gibbs现象。为有效地消除此现象，我们需对Ramp滤波器稍作平滑处理，如将之与作卷积，得到Shepp-Logan滤波器；,5.6图像投影重建,滤波反投影法中滤波器的选择,平滑了图像，损失了部分高频信息,Shepp-Logan滤波器,5.6图像投影重建,滤波反投影法中滤波器的选择,Hamming滤波器,降低了高频噪声,可得到Hamming滤波器和Hanning滤波器,5.6图像投影重建,反投影重建算法的matlab实现,I=iradon(R,theta)I=iradon(P,theta,interp,filter,frequency_scaling,output_size)I,H=iradon(.)DescriptionI=iradon(R,theta)reconstructstheimageIfromprojectiondatainthetwo-dimensionalarrayR.ThecolumnsofRareparallelbeamprojectiondata.iradonassumesthatthecenterofrotationisthecenterpointoftheprojections,whichisdefinedasceil(size(R,1)/2).thetadescribestheangles(indegrees)atwhichtheprojectionsweretaken.ItcanbeeitheravectorcontainingtheanglesorascalarspecifyingD_theta,theincrementalanglebetweenprojections.Ifthetaisavector,itmustcontainangleswithequalspacingbetweenthem.IfthetaisascalarspecifyingD_theta,theprojectionsweretakenatanglestheta=m*D_theta,wherem=0,1,2,.,size(R,2)-1.Iftheinputistheemptymatrix(),D_thetadefaultsto180/size(R,2).,5.6图像投影重建,反投影重建算法的matlab实现,I=iradon(P,theta,interp,filter,frequency_scaling,output_size)specifiesparameterstouseintheinverseRadontransform.Yerpspecifiesthetypeofinterpolationtouseinthebackprojection.Theavailableoptionsarelistedinorderofincreasingaccuracyandcomputationalcomplexity.,nearest:Nearest-neighborinterpolationLinear:Linearinterpolation(thedefault)spline:Splineinterpolationcubic:CubicinterpolationfromMATLAB5.,5.6图像投影重建,反投影重建算法的matlab实现,I=iradon(P,theta,interp,filter,frequency_scaling,output_size)specifiesparameterstouseintheinverseRadontransform.Yerpspecifiesthetypeofinterpolationtouseinthebackprojection.Theavailableoptionsarelistedinorderofincreasingaccuracyandcomputationalcomplexity.,nearest:Nearest-neighborinterpolationLinear:Linearinterpolation(thedefault)spline:Splineinterpolationcubic:CubicinterpolationfromMATLAB5.,5.6图像投影重建,反投影重建算法的matlab实现,I=iradon(P,theta,interp,filter,frequency_scaling,output_size)specifiesparameterstouseintheinverseRadontransform.Youcanspecifyanycombinationofthelastfourarguments.iradonusesdefaultvaluesforanyoftheseargumentsthatyouomit.filterspecifiesthefiltertouseforfrequencydomainfiltering.filtercanbeanyofthestringsthatspecifystandardfilters.,Ram-LakShepp-LoganCosineHammingHannnone,5.6图像投影重建,图像重建example,clear;closeall;clc;g=phantom(ModifiedShepp-Logan,600);subplot(241);imshow(g,);title(originalimage)theta=0:0.5:179.5;R,xp=radon(g,theta);subplot(242);imshow(R,);title(Projectionimage)f1=iradon(R,theta,none);subplot(243);imshow(f1,);title(ReconstructedImagewithbackprojection)f2=iradon(R,theta,Ram-Lak);subplot(244);imshow(f2,);title(ReconstructedImagewithRam-Lakfilter),5.6图像投影重建,图像重建example,f3=iradon(R,theta,Shepp-Logan);subplot(245);imshow(f3,);title(ReconstructedImagewithShepp-Loganfilter)f4=iradon(R,theta,Cosine);subplot(246);imshow(f4,);title(ReconstructedImagewithCosinefilter)f5=iradon(R,theta,Hamming);subplot(247);imshow(f5,);title(ReconstructedImagewithHammingfilter)f6=iradon(R,theta,Hann);subplot(248);imshow(f6,);title(ReconstructedImagewithHannfilter),5.6图像投影重建,5.6图像投影重建,扇形扫描重建,成像几何,5.6图像投影重建,扇形扫描重建,扇束情况下的重建算法较为复杂，但实质没有改变。可采用平行束情况下的算法实现，只需加以适当地修正即可,重排算法:把一个视图中采得的扇形数据重新组合成平行的射线投影数据，然后采用平行束重建算法重建直接重建算法:不必数据重排，只需适当加权即可运用与平行束类似的算法重建,5.6图像投影重建,扇形束重建-EA等角扇形束重建,当同样大小的探测器单元沿着中心为X射线焦点的弧排列时，就形成等角采样扇形束的每一条射线可由和确定，其中是射线与中心射线（假想的通过X射线源和等中心的直线）的夹角，称为探测器角；是中心射线与y轴的夹角，称为投影角,5.6图像投影重建,扇形束重建-EA等角扇形束重建,投影乘以探测器角的余弦，滤波后的样本随着到光源的距离的增长而增长,重建公式可由用(t,)坐标确定(,)坐标上的每个样本来得到。扇形投影中的投影样本q(,)就转化为平行投影中的投影样本p(t,),5.6图像投影重建,扇形投影和重建的matlab实现,fanbeam：Fan-beamtransformSyntax:F=fanbeam(I,D)F=fanbeam(.,param1,val1,param1,val2,.)F,fan_sensor_positions,fan_rotation_angles=fanbeam(.),FanRotationIncrement-Positiverealscalarspecifyingtheincrementoftherotationangleofthefan-beamprojections.Measuredindegrees.Defaultvalueis1.FanSensorGeometry-Textstringspecifyinghowsensorsarepositioned.Validvaluesarearcorline.Inthearcgeometry,sensorsarespacedequallyalongacirculararc,Thisisthedefaultvalue.FanSensorSpacing-Positiverealscalarspecifyingthespacingofthefan-beamsensors.InterpretationofthevaluedependsonthesettingofFanSensorGeometry.IfFanSensorGeometryissettoarc(thedefault),thevaluedefinestheangularspacingindegrees.Defaultvalueis1.IfFanSensorGeometryisline,thevaluespecifiesthelinearspacing.Defaultvalueis1.,5.6图像投影重建,扇形投影和重建的matlab实现,fanbeam：Fan-beamtransformSyntax:F=fanbeam(I,D)F=fanbeam(.,param1,val1,param1,val2,.)F,fan_sensor_positions,fan_rotation_angles=fanbeam(.),FanRotationIncrement-Positiverealscalarspecifyingtheincrementoftherotationangleofthefan-beamprojections.Measuredindegrees.Defaultvalueis1.FanSensorGeometry-Textstringspecifyinghowsensorsarepositioned.Validvaluesarearcorline.Inthearcgeometry,sensorsarespacedequallyalongacirculararc,Thisisthedefaultvalue.,5.6图像投影重建,扇形投影和重建的matlab实现,ifanbeam：Inversefan-beamtransformSyntax：I=ifanbeam(F,D)I=ifanbeam(.,param1,val1,param2,val2,.)I,H=ifanbeam(.),fan2para:Convertfan-beamprojectionstoparallel-beamSyntax:P=fan2para(F,D)P=fan2para(.,param1,val1,param2,val2,.)P,parallel_locations,parallel_rotation_angles=fan2para(.),5.6图像投影重建,扇形投影和重建的matlab实现,ifanbeam：Inversefan-beamtransformSyntax：I=ifanbeam(F,D)I=ifanbeam(.,param1,val1,param2,val2,.)I,H=ifanbeam(.),fan2para:Convertfan-beamprojectionstoparallel-beamSyntax:P=fan2para(F,D)P=fan2para(.,param1,val1,param2,val2,.)P,parallel_locations,parallel_rotation_angles=fan2para(.),para2fan:Convertparallel-beamprojectionstofan-beamSyntax:F=para2fan(P,D)I=para2fan(.,param1,val1,param2,val2,.)F,fan_positions,fan_rotation_angles=fan2para(.),5.6图像投影重建,扇形投影和重建的matlab实现,ifanbeam：Inversefan-beamtransformSyntax：I=ifanbeam(F,D)I=ifanbeam(.,param1,val1,param2,val2,.)I,H=ifanbeam(.),fan2para:Convertfan-beamprojectionstoparallel-beamSyntax:P=fan2para(F,D)P=fan2para(.,param1,val1,param2,val2,.)P,parallel_locations,parallel_rotation_angles=fan2para(.),para2fan:Convertparallel-beamprojectionstofan-beamSyntax:F=para2fan(P,D)I=para2fan(.,param1,val1,param2,val2,.)F,fan_positions,fan_rotation_angles=fan2para(.),5.6图像投影重建,扇形