《西安交通大学数字图像处理第六次作业(共18页).docx》由会员分享，可在线阅读，更多相关《西安交通大学数字图像处理第六次作业(共18页).docx（18页珍藏版）》请在taowenge.com淘文阁网|工程机械CAD图纸|机械工程制图|CAD装配图下载|SolidWorks_CaTia_CAD_UG_PROE_设计图分享下载上搜索。
1、精选优质文档-倾情为你奉上数字图像与视频处理第 六 次 作 业姓名： / 班级： 学号： 提交日期： 摘 要： 图像的退化与复原是图像处理的重要分支，处理目标是改善图像，以达到显示或后续处理的某种需求。本文首先讨论图像退化与复原的一般化模型，在此基础上对图像的噪声退化和运动模糊进行具体实现，并选择不同参数以便于讨论各参数的具体作用。在图像复原方面，本文着重讨论维纳滤波和约束最小二乘方滤波，然后利用这两种方法对退化图像进行复原，对比两种方法的处理效果以得出各自的优缺点。对于以上各操作，本文在openCV3.1.0计算机视觉库的基础上利用C+编程实现，附录开源，以提供参考。关键词： openCV
2、图像退化 图像复原 运动模糊 维纳滤波 约束最小二乘方专心-专注-专业介绍图像恢复技术是图像处理领域一类重要的处理技术，与图像增强等其他基本图像处理技术类似，该技术也是以改善图像视觉质量为目的，所不同的是图像恢复过程需要根据指定的图像退化模型来完成。从该退化模型出发，对在某种情况下退化了的图像进行恢复，以获取未经退化的原始图像。换句话说，图像恢复的处理过程实际是对退化图像品质的提升，并通过图像品质的提升来达到图像在视觉上的改善。一般地，图像的退化模型可以表述为 其中g(x,y)表示退化之后的图像，h(x,y)表示退化函数，h(x,y)表示原始图像，n(x,y)表示加性噪声。在空域图像的退化过程
3、是通过与退化函数的卷积和与加性噪声的叠加实现的，若对上式两端同时进行傅里叶变换，那么退化过程的频域表示为在此退化模型的基础上对进行图像回复，只需在输出后添加一个复原滤波器即可，图像的退化与复原过程可以完整地用图1表示。图像地退化与复原既可以在空域进行，也可以在频域进行，但一般而言在频域进行要相对方便，本文地后续讨论均在频域进行。图1 图像的退化与复原模型1 图像的高斯噪声退化及复原在数字图像中，噪声主要来源于图像地获取和传输过程1，根据噪声分布特性的不同，可以将噪声分为很多类，如均匀噪声、瑞利噪声、伽马噪声、高斯噪声等，其中以高斯噪声最为常见。高斯噪声具有数学上的易处理性，它同时也可以粗略表征
4、实践中图像噪声的大多数情况。一维高斯随机变量的概率分布函数可以写为其中x表示灰度值，x和分别表示均值和方差，p(x)表示灰度x的概率。高斯分布具有单峰特征，其大致概率密度曲线如下图所示。图2 一维标准高斯分布的概率密度曲线1.1 Box-Muller算法尽管高斯分布性质好且易处理，但在计算机中产生满足高斯分布的随机数却很困难，计算机方便产生均匀分布的随机数。本文借助Box-Muller算法2实现均匀分布向高斯分布的转化。算法细节不作赘述，可用结果可以表述为：假设随机变量U1、U2在区间0,1内服从标准均匀分布，那么随机变量Z0和Z1服从均值为零，标准差为1的标准均匀分布。上式中R和为极坐标下的
5、随机变量，可以表示为1.2 退化结果本文在VS2013环境下，利用openCV计算机视觉库实现图像的高斯噪声退化及复原。编写函数generateGaussianNoise(double mean, double sigma)产生服从高斯分布的随机数，继而利用自编函数addGaussianNoise()(此处省略形参)为图像添加高斯噪声，实现图像的退化过程。最后利用openCV自带的平滑函数medianBlur()、blur()等，实现图像的中值滤波、均值滤波、空域高斯滤波和频域高斯低通滤波，完成图像的复原。源代码见附录代码1，退化及复原过程中的各图像如图3所示。图3 高斯噪声退化及复原结果，上
6、图从左至右、由上及下依次为(a)-(f)；(a)原始图像；(b)添加均值为50，方差为10的高斯噪声之后的退化图像；(c)对退化图像应用模板边长为5的中值滤波器的复原结果；(d)对退化图像应用边长为5的均值滤波器的复原结果；(e)边长为5的空域高斯平滑模板的复原结果；(f)半径为150像素的频域高斯低通滤波器的复原结果1.3 结果讨论观察以上实验结果可以发现，添加高斯噪声后，图像各个区域出现了许多噪声点，表现为区别于背景的或白或黑的细小噪点，这使得图像的显示质量下降。此外，添加噪声后的退化图像的亮度要高于原始图像，这是由于噪声参数造成的，由于所添加的噪声均值是50，而不是零，因此整体上相当于给
7、图像添加了一个正灰度分量，如果噪声均值改为0，那么图像的整体亮度则不会改变。本文之所以选取50作为噪声均值，是为了方便观察退化效果。在本文涉及的四种复原方法中，中值滤波、均值滤波、高斯空域模糊都属于空域模糊处理，高斯低通滤波器属于频域处理。对比四种复原手段，其效果大同小异，都在去除了部分噪声的同时，一定程度上模糊了图像。综合比较去噪效果和模糊损失，效果最好的复原方法是空域高斯平滑模板，因为噪声在空域添加，而高斯模板在空域又具备高斯函数的起伏特征，起到了保留大部分原始像素值和抑制大部分噪声值得作用。这些复原手段严格来说是不理想的，理想效果应是能近似得到原始图像，即图3(a)所示图像。2 图像的运
8、动模糊2.1 模糊函数的推导图像的运动模糊是由于在成像曝光过程中，物体和相机之间的相对运动，使得物体上同一点的光线在多个成像单元上引起响应造成的。由于相机的成像单元在成像期间完成一个积分过程，以表示模糊后的图像，T表示积分时间，运动模糊原理可以表示为对上式进行傅里叶变换变换，得利用傅里叶变换的移位性质，频域模糊函数可以表示为现在考虑最简单的情况，那就是如果物体在两个坐标轴方向的运动是匀速直线运动，那么最终的模糊函数可以写为此过程的建模结果与频域滤波形式相类似，在进行处理时可以参照频域滤波的步骤，以作为滤波器函数，对进行填充和中心化后的图像频谱进行直接相乘即可。需要注意的一点就是在生成时也要按照
9、频谱中心化的形式进行，即的大小应是原始图像尺寸的两倍，且默认中心为坐标原点。2.2 实验结果在openCV3.1.0计算机视觉库基础上实现图像的运动模糊退化。首先编写函数motionBlur(省略形参)，以原图像为参考，生成一个两倍大小的运动模糊滤波器矩阵，包括两个部分：实部和虚部。然后按照频域滤波步骤进行后续操作，此时滤波过程是复数运算，应注意两个通道的不同实虚对应关系的相乘相加，此时需要借助函数multiply()和addWeighted()来实现3。具体源程序见附录源代码2，不同参数时的模糊结果如下图所示。abcdef图4 图像的运动模糊退化，上图从左至右、由上及下一次为(a)-(f)；
10、（a）原始图像；（b）参数a=0.01,b=0.01,T=1时的运动模糊退化图像；（c）参数a=0.05,b=0.05,T=1时的运动模糊退化图像；（d）参数a=0.1,b=0.1,T=1时的运动模糊退化图像；（e）参数a=0，b=0.1,T=1时的运动模糊退化图像；（f）参数a=0.05,b=0.05,T=0.5时的运动模糊退化图像2.3 结果分析观察上图实验结果可以发现，本文较好的实现了图像的运动模糊。从不同参数的吃力结果可知，参数a，b分别决定了物体在x，y两个方向的运动速度，反映在图像中就是图像在两个方向上的模糊程度，a和b越大，图像的模糊程度越剧烈。对比图3的(c)和(f)可知，两幅
11、图像的参数a，b完全相同，但(f)的曝光时间T要比(c)小。可以看出曝光时间T越小，意味着成像单元的积分时间越短，那么相同反射强度的物体成像的亮度越暗，这就是(f)的亮度明显小于(c)的原因。因此，直观来说，积分时间决定了图像整体的亮度大小。当T小于1时，图像整体亮度会低于原图；当T大于1时，图像整体亮度会高于原图。在进行运动模糊时，图像的边缘出现了一定宽度的“暗带”，其宽度随运动速度的增大而增大，当运动速度较小时基本可以忽略，而当运动速度较大时，这种现象尤其明显。这是因为在整个积分时间内，这些区域只有很短的非零时间，大多是时间被积量为零，所以最终的灰度很小。这种现象在实际摄影成像中也会出现，
12、因此不需要进行更正。2.4 为模糊图像添加噪声在图4(c)的基础上结合1.2中方法，为运动模糊后的图像添加高斯噪声，此处选取运动模糊参数：a=0.05,b=0.05,T=1。由于两个步骤均已实现，此处将不再附录代码，结果见图5abcd图5 为运动模糊后的图像添加高斯噪声；（a）运动模糊退化后的图像；（b）添加均值为零，方差为10的高斯噪声；（c）添加均值为10，方差为10的高斯噪声；（d）添加均值为20，方差为10的高斯噪声 3 维纳滤波复原图像3.1 维纳滤波原理维纳滤波是一种基于最小均方误差准则、对平稳过程的最优估计器4。这种滤波器的输出与期望输出之间的均方误差为最小，因此，它是一个最佳滤
13、波系统。它可用于提取被平稳噪声所污染的信号。维纳滤波的误差度量由下式给出带入图像参量及退化函数，求解最小均方误差，得到最终的滤波器频域表达式其中常数K表示噪声功率谱与为退化图像功率谱比，即3.2 滤波复原结果维纳滤波实现原理同2中的运动模糊，基本思路都是图像的频域滤波，编写函数wienerFilter()构造维纳滤波器的实部和虚部，函数代码见附录代码3，复原结果如下图所示。图6 维纳滤波复原结果，上图从左至右依次是(a)-(c)；(a)未退化的图像；(b)运动模糊(a=0.05,b=0.05,T=1)并添加高斯噪声(均值位0，方差为10)的退化图像；(c)维纳滤波复原结果3.2 结果分析对比以
14、上三幅图像可知，维纳滤波能显著去除运动模糊，但是在未退化图像的功率谱未知的情况下，对噪声的去除效果并不理想，最终恢复出的图像与未退化图像相比仍有较大差距。因此在实际操作中，应考虑交互式地选取参数K，以取得最佳的恢复效果。4 约束最小二乘方滤波4.1 滤波原理维纳滤波是建立在最小化统计准则的基础上，因此在平均意义上它是最优的，而约束最小二乘方滤波对每一幅图像都能取得最优的效果。一般地，约束最小二乘方滤波器的函数表达式为其中表示图像的退化函数，是一个参数，为取得最佳滤波效果，其调整目标是使得是函数的傅里叶变换，是一个时域的模板：4.2 滤波复原效果实现方法同3.2，同样都是频域滤波操作的基本步骤。
15、对P进行傅里叶变换之前要将其填充到原始图像的二倍大小；为简单记，参数通过交互形式指定，实验结果如下图所示。图7 约束最小二乘方滤波滤波结果，上图从左至右依次是(a)-(c)；(a)未退化的图像；(b)运动模糊(a=0.05,b=0.05,T=1)并添加高斯噪声(均值位0，方差为10)的退化图像；(c)约束最小二乘方复原结果4.3 结果分析与维纳滤波不同，最小二乘方滤波对噪声的去除效果很好，最终的复原图像中噪声的影响已经非常微弱，在这一点上其性能优于维纳滤波，这通过比较图6（c）和图7（c）能够明显的体会出。但是最终的复原图像依然存在明显的模糊现象，最小二乘方滤波并不能有效恢复运动模糊（至少在本
16、文所做实验中是这样）。附录参考文献1 Rafael C Gonzelez，Richard E Woods. 数字图像处理. 阮秋琦，阮宇智，等/译.北京：电子工业出版社. 2016.2 Box G E P, Muller M E. A Note on the Generation of Random Normal DeviatesJ. /Annals of Mathematical Statistics, 1958, 29(2):610-611.3 Robert Langaniere. OpenCV计算机视觉编程攻略. 相初银译. 北京：人民邮电/出版社，2015.4 Kenneth R.Ca
17、stleman. 数字图像处理. 朱志刚，林学訚，石定机，等译. 北京：/电子工业出版社. 2001.源代码1. 图像的高斯噪声退化（openCV3.1.0 C+）#include #include #include opencv2/imgproc/imgproc.hpp#include #include #include #include using namespace cv;using namespace std;/Box-Muller法生成满足高斯分布的随机数；para mean=均值，sigma=标准差double generateGaussianNoise(double mean,
18、double sigma)/ 定义小值const double epsilon = std:numeric_limits:min();static double z0, z1;double u1, u2;/ 生成满足0-1均匀分布的随机数dou1 = rand() * (1.0 / RAND_MAX);u2 = rand() * (1.0 / RAND_MAX); while (u1 = epsilon);/ flag为真构造高斯随机变量Xz0 = sqrt(-2.0 * log(u1) * cos(2 * CV_PI * u2);z1 = sqrt(-2.0 * log(u1) * sin(
19、2 * CV_PI * u2);return z0 * sigma + mean;/ 为图像添加高斯噪声cv:Mat addGaussianNoise(cv:Mat& srcImage,double mean, double sigma)cv:Mat resultImage = srcImage.clone();int channels = resultImage.channels();int nRows = resultImage.rows;int nCols = resultImage.cols * channels;for (int i = 0; i nRows; +i)for (int
20、 j = 0; j nCols; +j)/ 添加高斯噪声int val = resultImage.ptr(i)j + generateGaussianNoise(mean, sigma);if (val 255)val = 255;resultImage.ptr(i)j = (uchar)val;return resultImage;/ 定义高斯平滑函数,besarKernel:模板边长void myGaussianBlur(const Mat& src, Mat& result, int besarKernel, double sigma)/ 高斯模板半径int kerR = besarK
21、ernel / 2;/ 定义高斯模板因子cv:Mat kernel = Mat_(besarKernel, besarKernel);/ 归一化参数double alpha = 1 / (2 * 3.1416 * sigma * sigma);double sum = 0;/ 模板函数生成for (int i = -kerR; i = kerR; i+)for (int j = -kerR; j = kerR; j+)kernel.at(i + kerR, j + kerR)= exp(-(j * j) + (i * i) /(sigma * sigma * 2) * alpha;sum =
22、sum + kernel.at(i + kerR, j + kerR);for (int i = -kerR; i = kerR; i+)for (int j = -kerR; j = kerR; j+)kernel.at(i + kerR, j + kerR)= kernel.at(i + kerR, j + kerR) / sum;/输出高斯模板for (int k = -kerR; k = kerR; +k)for (int n = -kerR; n = kerR; +n)cout kernel.at(kerR + k, kerR + n) ;cout endl;result = src
23、.clone();double pix;/ 遍历图像（i,j）为像素索引for (int i = kerR; i src.rows - kerR; +i)for (int j = kerR; j src.cols - kerR; +j)/ 窗搜索完成高斯模板平滑，（k,n)为模板索引pix = 0;for (int k = -kerR; k = kerR; +k)for (int n = -kerR; n = kerR; +n)pix = pix + src.at(i + k, j + n) *kernel.at(kerR + k, kerR + n);result.at(i - kerR,
24、j - kerR) = pix;int main()cv:Mat srcImage =cv:imread(lena.bmp, 1);if (!srcImage.data)return -1;cv:imshow(srcImage, srcImage);cv:Mat resultImage = addGaussianNoise(srcImage,50,sqrt(10);cout sqrt(10) endl;cv:imshow(noisedIamge, resultImage);cv:imwrite(lenaAddGaussianNoise.bmp, resultImage);cv:Mat medi
25、anFilterImage, meanFilterImage, gaussianFilterImage, geoMeanFilter;/中值滤波去噪medianBlur(resultImage, medianFilterImage, 5);cv:imshow(medianFilterImage, medianFilterImage);cv:imwrite(lenaMedianFilterImage_5.bmp, medianFilterImage);/均值滤波去噪blur(resultImage, meanFilterImage, cv:Size(5, 5), cv:Point(-1, -1)
26、;cv:imshow(meanFilterImage, meanFilterImage);cv:imwrite(lenaMeanFilterImage_5.bmp, meanFilterImage);/空域高斯滤波去噪GaussianBlur(resultImage, gaussianFilterImage, cv:Size(5, 5), 0, 0);cv:imshow(gaussianFilterImage, gaussianFilterImage);cv:imwrite(gaussianNoise_50_sqrt(10).bmp, gaussianFilterImage);cv:waitK
27、ey(0);return 0;2. 图像的运动模糊 （openCV3.1.0 C+）#include opencv2/core/core.hpp#include opencv2/imgproc/imgproc.hpp#include opencv2/highgui/highgui.hpp#include #include #define pi 3.using namespace cv;using namespace std;/ 根据图像生成退化函数void motionBlur(cv:Mat srcImage, cv:Mat rePart, cv:Mat imPart, double a, d
28、ouble b, double t)int m = 2* srcImage.rows;int n = 2* srcImage.cols;for (int i = 0; im; i+)for (int j = 0; jn; j+)int ii = i - m / 2;int jj = j - n / 2;double cof = pi*(ii*a + jj*b)+1e-10;rePart.at(i, j) = (t / cof)*sin(cof)*cos(cof);imPart.at(i, j) = (-1)*(t / cof)*sin(cof)*sin(cof);/填充M*N图像为P*Q,P=
29、2*M,Q=2*Nvoid imagePadding(cv:Mat srcImage, cv:Mat paddingImage)int m = srcImage.rows;int n = srcImage.cols;int p = 2 * m;int q = 2 * n;for (int i = 0; ip; i+)for (int j = 0; jq; j+)if (i m & j n)paddingImage.at(i, j) = srcImage.at(i, j);else paddingImage.at(i, j) = 0;/乘-1幂进行中心化void shiftToMiddle(cv
30、:Mat srcImage)srcImage.convertTo(srcImage, CV_32FC1);for (int i = 0; isrcImage.rows; i+) /中心化float *p = srcImage.ptr(i);for (int j = 0; jsrcImage.cols; j+)pj = pj * pow(-1, i + j);int main()cv:Mat srcImage = cv:imread(lena.bmp, 0);if (srcImage.empty()return -1;cv:imshow(srcImage, srcImage);int m = 2
31、 * srcImage.rows;int n = 2 * srcImage.cols;/生成频域运动模糊滤波器cv:Mat rePart(m, n, CV_32FC1);cv:Mat imPart(m, n, CV_32FC1);motionBlur(srcImage, rePart, imPart, 0.05, 0.05, 0.5);/填充为2倍大小cv:Mat paddingImage(m, n, CV_32FC1);imagePadding(srcImage, paddingImage);/乘-1幂变换到中心paddingImage.convertTo(paddingImage, CV_
32、32FC1);shiftToMiddle(paddingImage);/ 为傅立叶变换的结果(实部和虚部)分配存储空间cv:Mat planes = cv:Mat_(paddingImage),cv:Mat:zeros(paddingImage.size(), CV_32F) ;Mat completeI;/ 为延扩后的图像增添一个初始化为0的通道merge(planes, 2, completeI);/ 进行离散傅立叶变换dft(completeI, completeI);split(completeI, planes);/相乘cv:Mat rePart1(m, n, CV_32FC1);c
33、v:Mat rePart2(m, n, CV_32FC1);cv:Mat imPart1(m, n, CV_32FC1);cv:Mat imPart2(m, n, CV_32FC1);/生成结果实部multiply(planes0, rePart, rePart1);multiply(planes1, imPart, rePart2);addWeighted(rePart1, 1, rePart2, -1, 0, planes0);/生成结果虚部multiply(planes0, imPart, imPart1);multiply(planes1, rePart, imPart2);addWe
34、ighted(imPart1, 1, imPart2, 1, 0, planes1);merge(planes, 2, completeI);cv:Mat idftImage,idftImageConverted;/进行逆变换&去中心化cv:idft(completeI, idftImage, cv:DFT_SCALE | cv:DFT_REAL_OUTPUT);shiftToMiddle(idftImage);idftImage.convertTo(idftImageConverted, CV_8U);/拷贝结果到原图像大小画布cv:Mat resltImage(srcImage.rows,
35、 srcImage.cols, CV_32FC1);for (int i = 0; isrcImage.rows; i+)for (int j = 0; jsrcImage.cols; j+)resltImage.at(i, j) = idftImageConverted.at(i, j);resltImage.convertTo(resltImage, CV_8UC1);imwrite(motionBluredImage_a=0.05_b=0.05_T=0.5.bmp, resltImage);imshow(resltImage, resltImage);cv:waitKey(0);retu
36、rn 0;3. 维纳滤波复原（openCV3.1.0 C+）（基本结构同2，此处省略主函数，仅提供维纳滤波模板的生成函数）void wienerFilter(cv:Mat srcImage, cv:Mat rePart, cv:Mat imPart, double a, double b, double t,double k)int m = 2* srcImage.rows;int n = 2* srcImage.cols;for (int i = 0; im; i+)for (int j = 0; jn; j+)int ii = i - m / 2;int jj = j - n / 2;do
37、uble cof = pi*(ii*a + jj*b)+1e-10;double denominator = t* t* sin(cof) * sin(cof) + k* cof* cof;rePart.at(i, j) = (t* cof* sin(cof) * cos(cof) / (denominator);imPart.at(i, j) = (t * cof * sin(cof) * sin(cof) / (denominator);4. 约束最小乘方滤波 （openCV3.1.0 C+）（基本结构同2，此处省略主函数，仅提供维纳滤波模板的生成函数）void leastSquareFi
38、lter(cv:Mat srcImage, cv:Mat Pre, cv:Mat Pim, cv:Mat rePart, cv:Mat imPart, double a, double b, double t, double gama)int m = 2 * srcImage.rows;int n = 2 * srcImage.cols;for (int i = 0; im; i+)for (int j = 0; jn; j+)int ii = i - m / 2;int jj = j - n / 2;double cof = pi*(ii*a + jj*b) + 1e-10;double d
39、enominator1 = t* sin(cof) + gama*cof*Pre.at(i, j);double denominator2 = gama*cof*Pim.at(i, j);double denominator = denominator1*denominator1 + denominator2*denominator2;double re = t*sin(cof)*cos(cof);double im = t*sin(cof)*sin(cof);rePart.at(i, j) = (re*denominator1 + im*denominator2) / (denominator);imPart.at(i, j) = (im*denominator1 - re*denominator2) / (denominator);