智能图像处理：Python和OpenCV实现-课件 第七章 图像退化和复原

第七章图像退化和复原目录Contents图像退化与复原机理图像的运动模糊图像的逆滤波图像的维纳滤波图像质量的评价何为图像退化和复原图像的退化是指图像在形成、传输和记录过程中，由于成像系统、传输介质和设备的不完善，使图像的质量下降。针对这些问题，我们需要对退化后的图像进行复原，图像复原就是尽可能恢复退化图像的本来面目，它是沿图像退化的逆过程进行处理，也就是如果我们知道图像是经历了什么样的过程导致退化，就可以按其逆过程来复原图像。7.1图像退化和复原机理1.图像退化：图像退化就是图像因为某种原因从“正常”图像变得“不正常”，典型的表现有模糊、失真、有噪声等。图像退化模型：输入图像f(x,y)经过某个退化系统H(x,y)，再叠加上噪声n(x,y)，形成退化后的图像g(x,y)。H(x,y)+f(x,y)n(x,y)g(x,y)

复原函数g(x,y)

图像复原的步骤：确定参考图像，作为图像退化/复原模型的评估标准；设计图像退化算法，引入运动模糊和白噪声；传统算法的原理及编程实现；评价函数的设计及编程实现；过程评估和结果分析。7.2图像的运动模糊运动模糊图像是由同一图像在产生距离延迟后与原图像叠加而成。令x0(t)和y0(t)分别代表位移的x分量和y分量，那么在相机快门开启的时间T内，图像传感器上某点的总曝光量是图像在移动过程中一系列相应像素的亮度对该点作用的总和。即，如果快门开启与关闭的时间忽略不计，则有：

实现运动模糊的模型函数为motion_process(image_size,motion_angle)，它包含两个参数：图像的尺寸大小image_size以及运动的角度motion_angle。当运动位移为degree、运动角度为angle度时，则该模型函数的构建过程如下：首先是创建与图像同等大小的全0矩阵，然后找到全0矩阵的中心行数center_position，再计算出运动角度的正切（tan）值与余切（cot）值，算出运动的偏移量offset。再令α≤45°时，PSF[int(center_position+offset),int(center_position-offset)]=1α>45°时，PSF[int(center_position-offset),int(center_position+offset)]=1其中PSF为点扩散函数。【例7.1】根据运动模糊的模型函数生成任意角度的运动模糊图像。程序代码如下：importcv2importmathimportnumpyasnpimportmatplotlib.pyplotasplt#仿真运动模糊defmotion_process(image_size,motion_angle):PSF=np.zeros(image_size)center_position=(image_size[0]-1)/2

slope_tan=math.tan(motion_angle*math.pi/180)slope_cot=1/slope_tanifslope_tan<=1:foriinrange(15):offset=round(i*slope_tan)PSF[int(center_position+offset),int(center_position-offset)]=1returnPSF/PSF.sum()#对点扩散函数进行归一化亮度else:foriinrange(15):offset=round(i*slope_cot)PSF[int(center_position-offset),int(center_position+offset)]=1returnPSF/PSF.sum()#对图像进行运动模糊defmake_blurred(input,PSF,eps):input_fft=np.fft.fft2(input)#二维数组的傅里叶变换PSF_fft=np.fft.fft2(PSF)+epsblurred=np.fft.ifft2(input_fft*PSF_fft)blurred=np.abs(np.fft.fftshift(blurred))returnblurredimg=cv2.imread('d:/pics/lena.jpg')img_gray=cv2.cvtColor(img,cv2.COLOR_BGR2GRAY)plt.subplot(121),plt.axis('off')plt.title("Originimage"),plt.imshow(img_gray)#进行运动模糊处理img_h,img_w=img.shape[0:2]PSF=motion_process((img_h,img_w),60)blurred=np.abs(make_blurred(img_gray,PSF,1e-3))plt.subplot(122),plt.axis('off')plt.title("Motionblurred")plt.imshow(blurred)程序运行结果如图7-1所示原图像(b)运动模糊后的图像图7-1运动模糊图像【例7.2】根据运动模糊核（kernel）矩阵生成运动模糊图像。程序代码如下：importcv2importnumpyasnp

defmotion_blur(image,degree=12,angle=45):image=np.array(image)#生成任意角度的运动模糊kernel的矩阵，degree越大，模糊程度越高M=cv2.getRotationMatrix2D((degree/2,degree/2),angle,1)motion_blur_kernel=np.diag(np.ones(degree))#输出矩阵的对角线元素#放射变换函数motion_blur_kernel=cv2.warpAffine(motion_blur_kernel,M,(degree,degree))motion_blur_kernel=motion_blur_kernel/degreeblurred=cv2.filter2D(image,-1,motion_blur_kernel)cv2.normalize(blurred,blurred,0,255,cv2.NORM_MINMAX)#归一化函数blurred=np.array(blurred,dtype=np.uint8)returnblurred

img=cv2.imread('d:/pics/lena.jpg')img_blurred=motion_blur(img)cv2.imshow('Originimage',img)cv2.imshow('Blurredimage',img_blurred)cv2.waitKey(0)cv2.destroyAllWindows()程序运行结果如图7-2所示。(a)原图像(b)运动模糊后的图像图7-2运动模糊图像7.3图像的逆滤波

【例7.3】使用逆滤波算法对运动模糊图像进行滤波。程序代码如下：importcv2importmathimportnumpyasnpimportmatplotlib.pyplotasplt

#仿真运动模糊defmotion_process(image_size,motion_angle):PSF=np.zeros(image_size)center_position=(image_size[0]-1)/2

slope_tan=math.tan(motion_angle*math.pi/180)slope_cot=1/slope_tanifslope_tan<=1:foriinrange(15):offset=round(i*slope_tan)PSF[int(center_position+offset),int(center_position-offset)]=1returnPSF/PSF.sum()#对点扩散函数进行归一化亮度else:foriinrange(15):offset=round(i*slope_cot)PSF[int(center_position-offset),int(center_position+offset)]=1returnPSF/PSF.sum()#对图片进行运动模糊defmake_blurred(input,PSF,eps):input_fft=np.fft.fft2(input)#进行二维数组的傅里叶变换PSF_fft=np.fft.fft2(PSF)+epsblurred=np.fft.ifft2(input_fft*PSF_fft)blurred=np.abs(np.fft.fftshift(blurred))returnblurreddefinverse(input,PSF,eps):#逆滤波input_fft=np.fft.fft2(input)PSF_fft=np.fft.fft2(PSF)+eps#噪声功率，这是已知的，考虑epsilonresult=np.fft.ifft2(input_fft/PSF_fft)#计算F(u,v)的傅里叶反变换result=np.abs(np.fft.fftshift(result))returnresultif__name__=='__main__':image=cv2.imread('d:/pics/lena.jpg')image=cv2.cvtColor(image,cv2.COLOR_BGR2GRAY)

#进行运动模糊处理img_h,img_w=image.shape[0:2]PSF=motion_process((img_h,img_w),60)blurred=np.abs(make_blurred(image,PSF,1e-3))plt.subplot(221),plt.axis('off')plt.title("Motionblurred")plt.imshow(blurred)result_inv=inverse(blurred,PSF,1e-3)#逆滤波plt.subplot(222),plt.axis('off')plt.title("inversedeblurred")plt.imsave('d:\pics\lenaInv.jpg',result_inv)plt.imshow(result_inv)#添加噪声,standard_normal产生随机的函数blurred_noisy=blurred+0.1*blurred.std()*np.random.standard_normal(blurred.shape)plt.subplot(223),plt.axis('off')plt.title("motion&noisyblurred")plt.imshow(blurred_noisy)#显示添加噪声且运动模糊的图像

#对添加噪声的图像进行逆滤波result_invn=inverse(blurred_noisy,PSF,0.1+1e-3)plt.subplot(224),plt.axis('off')plt.title("Nioseinversedeblurred")plt.imsave('d:\pics\lenaInvn.jpg',result_invn)plt.imshow(result_invn)plt.show()程序运行结果如图7-3所示。图中左上图为运动模糊图像，右上图为对左上图逆滤波后的图像；左下图像为运动模糊+噪声图像，右下图为对左下图逆滤波后的图像。图7-3运动模糊逆滤波图像7.4图像的维纳滤波

其中，H(u,v)为退化函数，H*(u,v)为H(u,v)的复共轭，|H(u,v)|2=H(u,v)H*(u,v)，Sη(u,v)=|N(u,v)|2=噪声的功率谱，Sf(u,v)=|F(u,v)|2=未退化图像的功率谱，H(u,v)是退化函数的傅里叶变换，G(u,v)是退化后图像的傅里叶变换。从上面的公式可以发现，如没有噪声，即Sη(u,v)=0，此时维纳滤波变为直接逆滤波，如有噪声，那么Sη(u,v)如何估计将成问题，同时Sf(u,v)的估计也成问题。在实际应用中假设退化函数已知，如果噪声为高斯白噪声，则Sη(u,v)为常数，但Sf(u,v)通常难以估计。一种近似的解决办法是用一个系数K代替Sη(u,v)/Sf(u,v)，因此公式7-5变为公式7-6所示：在假设噪声和图像不相关，其中一个或另一个有零均值，且估计中的灰度级是退化图像中灰度级的线性函数的条件下，均方误差函数的最小值在频率域有如公式7-5所示：

【例7.4】使用维纳滤波对运动模糊图像进行滤波处理。程序代码如下：importcv2importmathimportnumpyasnpimportmatplotlib.pyplotasplt

#仿真运动模糊defmotion_process(image_size,motion_angle):PSF=np.zeros(image_size)center_position=(image_size[0]-1)/2slope_tan=math.tan(motion_angle*math.pi/180)slope_cot=1/slope_tanifslope_tan<=1:foriinrange(15):offset=round(i*slope_tan)PSF[int(center_position+offset),int(center_position-offset)]=1returnPSF/PSF.sum()#对点扩散函数进行归一化亮度

else:foriinrange(15):offset=round(i*slope_cot)PSF[int(center_position-offset),int(center_position+offset)]=1returnPSF/PSF.sum()#对图像进行运动模糊defmake_blurred(input,PSF,eps):input_fft=np.fft.fft2(input)#二维数组的傅里叶变换PSF_fft=np.fft.fft2(PSF)+epsblurred=np.fft.ifft2(input_fft*PSF_fft)blurred=np.abs(np.fft.fftshift(blurred))returnblurreddefwiener(input,PSF,eps,K=0.01):#维纳滤波，K=0.01input_fft=np.fft.fft2(input)PSF_fft=np.fft.fft2(PSF)+epsPSF_fft_1=np.conj(PSF_fft)/(np.abs(PSF_fft)**2+K)result=np.fft.ifft2(input_fft*PSF_fft_1)result=np.abs(np.fft.fftshift(result))returnresultif__name__=='__main__':image=cv2.imread('d:/pics/lena.jpg')image=cv2.cvtColor(image,cv2.COLOR_BGR2GRAY)

#进行运动模糊处理img_h,img_w=image.shape[0:2]PSF=motion_process((img_h,img_w),60)blurred=np.abs(make_blurred(image,PSF,1e-3))plt.subplot(221),plt.axis('off')plt.title("Motionblurred")plt.imshow(blurred)plt.imsave('d:\pics\lenaMb.jpg',blurred)

resultwd=wiener(blurred,PSF,1e-3)#维纳滤波plt.subplot(222),plt.axis('off')plt.title("wienerdeblurred(k=0.01)")plt.imshow(resultwd)plt.imsave('d:\pics\lenaWd.jpg',resultwd)

#添加噪声,standard_normal产生随机的函数

blurred_noisy=blurred+0.1*blurred.std()*np.random.standard_normal(blurred.shape)plt.subplot(223),plt.axis('off')plt.title("motion&noisyblurred")plt.imshow(blurred_noisy)#显示添加噪声且运动模糊的图像

#对添加噪声的图像进行维纳滤波resultwdn=wiener(blurred_noisy,PSF,0.1+1e-3)plt.subplot(224),plt.axis('off')plt.title("wienerdeblurred(k=0.01)")plt.imshow(resultwdn)plt.show()程序运行结果如图7-4所示。图中左上图为运动模糊图像，右上图为对左上图维纳滤波后的图像；左下图像为运动模糊+噪声图像，右下图为对左下图维纳滤波后的图像。图7-4运动模糊维纳滤波后的图像在程序中可以修改运动模糊的角度，观察不同的角度取值情况下的模糊效果。对比7-3的逆滤波和图7-4的维纳滤波，从效果来看维纳滤波对图像去运动模糊效果较好。7.5图像质量的评价图像修复模型的好坏需要对结果品质进行评价，衡量图像修复的指标主要有以下三种：主观评价方法：由观测者对图像质量进行主观评估，大部分在实际生产中好多师傅会采用此法，虽不科学，但这种一般是最终决定产品好坏的方法。主观质量评分法是图像质量最具代表性的主观评价方法，它通过对观察者的评价判断图像质量。主观评价方法与人的主观感受相符，但费时、复杂，还会受到观测者专业背景、心理和动机等主观因素的影响，并且不能结合到其他算法中使用。客观评价方法：图像质量的客观评价方法是根据人眼的主观视觉系统建立数学模型，并通过具体的公式计算图像的质量。客观评价方法方便、快捷，容易实现并能结合到应用系统中，但它和人的主观感受有出入。通常提到的图像质量评价算法是指客观评价算法，其目标是获得与主观评价结果相一致的客观评价值。传统的图像质量客观评价方法主要包括均方误差（MSE，meansquarederror）和峰值信噪比（PSNR，peaksignaltonoiserate）。均方误差（MSE）首先计算原始图像和失真图像像素差值的均方值，然后通过均方值的大小来确定失真图像的失真程度。基于Python定义的均方误差MSE函数为：defmse(img1,img2):mse=np.mean((img1/255.-img2/255.)**2)returnmse

其中MaxValue2表示为存储最大位数2bits−1，对图像而言bits=8，MaxValue=255。PSNR指标越高，说明图像质量越好。基于Python的定义的峰值信噪比PSNR函数为：defpsnr(img1,img2):mse=np.mean((img1/1.0-img2/1.0)**2)ifmse<1.0e-10:return100return10*math.log10(255.0**2/mse)3.结构相似性评价方法（SSIM，StructuralSimilarity）。该评价方法是在2001年提出的，它是一种衡量电视、电影或者其他数字图像、视频的主观感受质量的一种方法。SSIM算法是用来测试两幅图像的相似性，其测量或者预测图像的质量是基于未压缩的或者无失真的图像作为参考的。SSIM算法在设计上考虑了人眼的视觉特性，比传统方式更符合人眼视觉感知，它从自然图像高度结构化的特征出发，通过亮度(luminance)、对比度(contrast)和结构(structure)三个方面估计感知结构信息的变化。结构信息是指像素之间有着内部的依赖性，尤其是空间上靠近的像素点，这些依赖性携带着目标对象视觉感知上的重要信息。

结构相似性指标也是一种用以衡量两张数字图像相似程度的指标。当两张图像其中一张为无失真图像，另一张为失真后的图像，二者的结构相似性可以看成是失真图像的图像品质衡量指标。相较于传统所使用的图像品质衡量指标(如MSE、PSNR)，结构相似性在图像品质的衡量上更能符合人眼对图像品质的判断。【例7.5】使用numpy库函数对维纳滤波图像求解MSE、PSNR和SSIM评价参数。程序代码如下：importcv2importmathimportnumpyasnpfromscipy.signalimportconvolve2ddefcompute_mse(img1,img2):mse=np.mean((img1/255.-img2/255.)**2)returnmsedefcompute_psnr(img1,img2):mse=np.mean((img1/1.0-img2/1.0)**2)ifmse<1.0e-10:return100return10*math.log10(255.0**2/mse)defmatlab_style_gauss2D(shape=(3,3),sigma=0.5):m,n=[(ss-1.)/2.forssinshape]y,x=np.ogrid[-m:m+1,-n:n+1]h=np.exp(-(x*x+y*y)/(2.*sigma*sigma))h[h<np.finfo(h.dtype).eps*h.max()]=0sumh=h.sum()ifsumh!=0:h/=sumhreturnhdeffilter2(x,kernel,mode='same'):returnconvolve2d(x,np.rot90(kernel,2),mode=mode)defcompute_ssim(im1,im2,k1=0.01,k2=0.03,win_size=11,L=255):ifnotim1.shape==im2.shape:raiseValueError("InputImageesmusthavethesamedimensions")iflen(im1.shape)>2:raiseValueError("Pleaseinputtheimageswith1channel")M,N=im1.shapeC1=(k1*L)**2C2=(k2*L)**2window=matlab_style_gauss2D(shape=(win_size,win_size),sigma=1.5)window=window/np.sum(np.sum(window))ifim1.dtype==np.uint8:im1=np.double(im1)ifim2.dtype==np.uint8:im2=np.double(im2)mu1=filter2(im1,window,'valid')mu2=filter2(im2,window,'valid')mu1_sq=mu1*mu1mu2_sq=mu2*mu2mu1_mu2=mu1*mu2sigma1_sq=filter2(im1*im1,window,'valid')-mu1_sqsigma2_sq=filter2(im2*im2,window,'valid')-mu2_sqsigmal2=filter2(im1*im2,window,'valid')-mu1_mu2ssim_map=((2*mu1_mu2+C1)*(2*sigmal2+C2))/((mu1_sq+mu2_sq+C1)*(sigma1_sq+sigma2_sq+C2))returnnp.mean(np.mean(ssim_map))if__name__=='__main__':img_origin=cv2.imread('d:/pics/lenaorig.jpg',0)cv2.imshow('Originimage',img_origin)img_Mb=cv2.imread('d:/pics/lenaMb.jpg',0)cv2.imshow('Moveblurimage',img_Mb)img_wd=cv2.imread('d:/pics/lenaWd.jpg',0)cv2.imshow('Winerfilterimage',img_wd)img_wdn=cv2.imread('d:/pics/lenaWdn.jpg',0)cv2.imshow('nWinernoisefilterimage',img_wdn)

#原图与运动模糊后的维纳滤波图像之间的评价参数mse=compute_mse(img_origin,img_wd)print('MSE:{}'.format(mse))psnr=compute_psnr(img_origin,img_wd)print('PSNR:{}'.format(psnr))ssim=compute_ssim(img_origin,img_wd)print('SSIM:{}'.format(ssim))#原图与运动模糊+噪声后的维纳滤波图像之间的评价参数mse=compute_mse(img_origin,img_wdn)print('MSE:{}'.format(mse))psnr=compute_psnr(img_origin,img_wdn)print('PSNR:{}'.format(psnr))ssim=compute_ssim(img_origin,img_wdn)print('SSIM:{}'.format(ssim))

cv2.waitKey(0)cv2.destroyAllWindows()（1）基于运动模糊图像的维纳滤波后的参数为：MSE:0.001349162245290273PSNR:28.69935820530212SSIM:0.8315891237838533(d)带有噪声的维纳滤波图像(a)原图像(b)运动模糊图像图7-5运动模糊维纳滤波图像(c)维纳滤波图像程序运行结果如图7-5所示，其输出的图像评价参数为：（2）运动模糊+噪声图像的维纳滤波后的参数为：MSE:0.010067537870049981PSNR:19.970767280753556SSIM:0.5047339796956921【例7.6】使用numpy库函数对逆滤波图像求解MSE、PSNR和SSIM评价参数图像质量的评价。程序代码如下：importcv2importmathimportnumpyasnpfromscipy.signalimportconvolve2ddefcompute_mse(img1,img2):mse=np.mean((img1/255.-img2/255.)**2)returnmsedefcompute_psnr(img1,img2):mse=np.mean((img1/1.0-img2/1.0)**2)ifmse<1.0e-10:return100return10*math.log10(255.0**2/mse)defmatlab_style_gauss2D(shape=(3,3),sigma=0.5):m,n=[(ss-1.)/2.forssinshape]y,x=np.ogrid[-m:m+1,-n:n+1]h=np.exp(-(x*x+y*y)/(2.*sigma*sigma))h[h<np.finfo(h.dtype).eps*h.max()]=0sumh=h.sum()ifsumh!=0:h/=sumhreturnhdeffilter2(x,kernel,mode='same'):returnconvolve2d(x,np.rot90(kernel,2),mode=mode)

defcompute_ssim(im1,im2,k1=0.01,k2=0.03,win_size=11,L=255):ifnotim1.shape==im2.shape:raiseValueError("InputImageesmusthavethesamedimensions")iflen(im1.shape)>2:raiseValueError("Pleaseinputtheimageswith1channel")M,N=im1.shapeC1=(k1*L)**2C2=(k2*L)**2window=matlab_style_gauss2D(shape=(win_size,win_size),sigma=1.5)window=window/np.sum(np.sum(window))ifim1.dtype==np.uint8:im1=np.double(im1)ifim2.dtype==np.uint8:im2=np.double(im2)mu1=filter2(im1,window,'valid')mu2=filter2(im2,window,'valid')mu1_sq=mu1*mu1mu2_sq=mu2*mu2mu1_mu2=mu1*mu2sigma1_sq=filter2(im1*im1,window,'valid')-mu1_sqsigma2_sq=filter2(im2*im2,window,'valid')-mu2_sqsigmal2=filter2(im1*im2,window,'valid')-mu1_mu2ssim_map=((2*mu1_mu2+C1)*(2*sigmal2+C2))/((mu1_sq+mu2_sq+C1)*(sigma1_sq+sigma2_sq+C2))returnnp.mean(np.mean(ssim_map))if__name__=='__main__':img_origin=cv2.imread('d:/pics/lenaorig.jpg',0)cv2.imshow('Originimage',img_origin)img_Mb=cv2.imread('d:/pics/lenaMb.jpg',0)cv2.imshow('Moveblurimage',img_Mb)

img_inv=cv2.imread('d:/pics/lenainv.jpg',0)cv2.imshow('Inversefilterimage',img_inv)img_invn=cv2.imread('d:/pics/lenainvn.jpg',0)cv2.imshow('Noiseinversefilterimage',img_invn)#原图与运动模糊后的逆滤波图像之间的评价参数mse=compute_mse(img_origin,img_inv)print('MSE:{}'.format(mse))psnr=compute_psnr(img_origin,img_inv)print('PSNR:{}'.format(psnr))ssim=compute_ssim(img_origin,img_inv)print('SSIM:{}'.format(ssim))#原图与运动模糊+噪声后的逆滤波图像之间的评价参数mse=compute_mse(img_origin,img_invn)print('MSE:{}'.format(mse))psnr=compute_psnr(img_origin,img_invn)print('PSNR:{}'.format(psnr))ssim=compute_ssim(img_origin,img_invn)print('SSIM:{}'.format(ssim))cv2.waitKey(0)cv2.destroyAllWindows()（1）基于运动模糊图像的逆滤波后的参数为：MSE:1.2172241445597847e-06PSNR:59.146294416695355SSIM:0.9995992013825972(a)原图像(b)运动模糊图像(c)逆滤波图像(d)带有噪声的逆滤波图像图7-6运动模糊逆滤波图像（2）运动模糊+噪声图像的逆滤波后的参数为：MSE:0.04287537024221453PSNR:13.677921166644737SSIM:0.18683314611701046程序运行结果如图7-6所示，其输出的图像评价参数为：【例7.7】使用skimage库函数对维纳滤波图像求解MSE、PSNR和SSIM评价参数。程序代码如下：importcv2fromskimage.measureimportcompare_msefromskimage.measureimportcompare_psnrfromskimage.measureimportcompare_ssim

img_origin=cv2.imread('d:/pics/lenaorig.jpg',0)cv2.imshow('Originimage',img_origin)img_Mb=cv2.imread('d:/pics/lenaMb.jpg',0)cv2.imshow('Moveblurimage',img_Mb)img_wd=cv2.imread('d:/pics/lenaWd.jpg',0)cv2.imshow('Winerfilterimage',img_wd)img_wdn=cv2.imread('d:/pics/lenaWdn.jpg',0)cv2.imshow('nWinernoisefilterimage',img_wdn)#原图与运动模糊后的维纳滤波图像之间的评价参数mse=compare_mse(img_origin,img_wd)print('MSE:{}'.format(mse))psnr=compare_psnr(img_origin,img_wd)print('PSNR:{}'.format(psnr))ssim=compare_ssim(img_origin,img_wd)print('SSIM:{}'.format(ssim))#原图与运动模糊+噪声后的维纳滤波图像之间的评价参数mse=compare_mse(img_origin,img_wdn)print('MSE:{}'.format(mse))psnr=compare_psnr(img_origin,img_wdn)print('PSNR:{}'.format(psnr))ssim=compare_ssim(img_origin,img_wdn)print('SSIM:{}'.format(ssim))cv2.waitKey(0)cv2.destroyAllWindows()(d)带有噪声的维纳滤波图像(a)原图像(b)运动模糊图像图7-7运动模糊维纳滤波图像(c)维纳滤波图像