西安交通大學(xué)數(shù)字圖像處理第六次作業(yè).docx

數(shù)字圖像與視頻處理第 六 次 作 業(yè)姓名： / 班級： 學(xué)號： 提交日期： 摘 要： 圖像的退化與復(fù)原是圖像處理的重要分支，處理目標(biāo)是改善圖像，以達到顯示或后續(xù)處理的某種需求。本文首先討論圖像退化與復(fù)原的一般化模型，在此基礎(chǔ)上對圖像的噪聲退化和運動模糊進行具體實現(xiàn)，并選擇不同參數(shù)以便于討論各參數(shù)的具體作用。在圖像復(fù)原方面，本文著重討論維納濾波和約束最小二乘方濾波，然后利用這兩種方法對退化圖像進行復(fù)原，對比兩種方法的處理效果以得出各自的優(yōu)缺點。對于以上各操作，本文在openCV3.1.0計算機視覺庫的基礎(chǔ)上利用C+編程實現(xiàn)，附錄開源，以提供參考。關(guān)鍵詞： openCV 圖像退化 圖像復(fù)原 運動模糊 維納濾波 約束最小二乘方18介紹圖像恢復(fù)技術(shù)是圖像處理領(lǐng)域一類重要的處理技術(shù)，與圖像增強等其他基本圖像處理技術(shù)類似，該技術(shù)也是以改善圖像視覺質(zhì)量為目的，所不同的是圖像恢復(fù)過程需要根據(jù)指定的圖像退化模型來完成。從該退化模型出發(fā)，對在某種情況下退化了的圖像進行恢復(fù)，以獲取未經(jīng)退化的原始圖像。換句話說，圖像恢復(fù)的處理過程實際是對退化圖像品質(zhì)的提升，并通過圖像品質(zhì)的提升來達到圖像在視覺上的改善。一般地，圖像的退化模型可以表述為 其中g(shù)(x,y)表示退化之后的圖像，h(x,y)表示退化函數(shù)，h(x,y)表示原始圖像，n(x,y)表示加性噪聲。在空域圖像的退化過程是通過與退化函數(shù)的卷積和與加性噪聲的疊加實現(xiàn)的，若對上式兩端同時進行傅里葉變換，那么退化過程的頻域表示為在此退化模型的基礎(chǔ)上對進行圖像回復(fù)，只需在輸出后添加一個復(fù)原濾波器即可，圖像的退化與復(fù)原過程可以完整地用圖1表示。圖像地退化與復(fù)原既可以在空域進行，也可以在頻域進行，但一般而言在頻域進行要相對方便，本文地后續(xù)討論均在頻域進行。圖1 圖像的退化與復(fù)原模型1 圖像的高斯噪聲退化及復(fù)原在數(shù)字圖像中，噪聲主要來源于圖像地獲取和傳輸過程1，根據(jù)噪聲分布特性的不同，可以將噪聲分為很多類，如均勻噪聲、瑞利噪聲、伽馬噪聲、高斯噪聲等，其中以高斯噪聲最為常見。高斯噪聲具有數(shù)學(xué)上的易處理性，它同時也可以粗略表征實踐中圖像噪聲的大多數(shù)情況。一維高斯隨機變量的概率分布函數(shù)可以寫為其中x表示灰度值，x和分別表示均值和方差，p(x)表示灰度x的概率。高斯分布具有單峰特征，其大致概率密度曲線如下圖所示。圖2 一維標(biāo)準(zhǔn)高斯分布的概率密度曲線1.1 Box-Muller算法盡管高斯分布性質(zhì)好且易處理，但在計算機中產(chǎn)生滿足高斯分布的隨機數(shù)卻很困難，計算機方便產(chǎn)生均勻分布的隨機數(shù)。本文借助Box-Muller算法2實現(xiàn)均勻分布向高斯分布的轉(zhuǎn)化。算法細(xì)節(jié)不作贅述，可用結(jié)果可以表述為：假設(shè)隨機變量U1、U2在區(qū)間0,1內(nèi)服從標(biāo)準(zhǔn)均勻分布，那么隨機變量Z0和Z1服從均值為零，標(biāo)準(zhǔn)差為1的標(biāo)準(zhǔn)均勻分布。上式中R和為極坐標(biāo)下的隨機變量，可以表示為1.2 退化結(jié)果本文在VS2013環(huán)境下，利用openCV計算機視覺庫實現(xiàn)圖像的高斯噪聲退化及復(fù)原。編寫函數(shù)generateGaussianNoise(double mean, double sigma)產(chǎn)生服從高斯分布的隨機數(shù)，繼而利用自編函數(shù)addGaussianNoise()(此處省略形參)為圖像添加高斯噪聲，實現(xiàn)圖像的退化過程。最后利用openCV自帶的平滑函數(shù)medianBlur()、blur()等，實現(xiàn)圖像的中值濾波、均值濾波、空域高斯濾波和頻域高斯低通濾波，完成圖像的復(fù)原。源代碼見附錄代碼1，退化及復(fù)原過程中的各圖像如圖3所示。圖3 高斯噪聲退化及復(fù)原結(jié)果，上圖從左至右、由上及下依次為(a)-(f)；(a)原始圖像；(b)添加均值為50，方差為10的高斯噪聲之后的退化圖像；(c)對退化圖像應(yīng)用模板邊長為5的中值濾波器的復(fù)原結(jié)果；(d)對退化圖像應(yīng)用邊長為5的均值濾波器的復(fù)原結(jié)果；(e)邊長為5的空域高斯平滑模板的復(fù)原結(jié)果；(f)半徑為150像素的頻域高斯低通濾波器的復(fù)原結(jié)果1.3 結(jié)果討論觀察以上實驗結(jié)果可以發(fā)現(xiàn)，添加高斯噪聲后，圖像各個區(qū)域出現(xiàn)了許多噪聲點，表現(xiàn)為區(qū)別于背景的或白或黑的細(xì)小噪點，這使得圖像的顯示質(zhì)量下降。此外，添加噪聲后的退化圖像的亮度要高于原始圖像，這是由于噪聲參數(shù)造成的，由于所添加的噪聲均值是50，而不是零，因此整體上相當(dāng)于給圖像添加了一個正灰度分量，如果噪聲均值改為0，那么圖像的整體亮度則不會改變。本文之所以選取50作為噪聲均值，是為了方便觀察退化效果。在本文涉及的四種復(fù)原方法中，中值濾波、均值濾波、高斯空域模糊都屬于空域模糊處理，高斯低通濾波器屬于頻域處理。對比四種復(fù)原手段，其效果大同小異，都在去除了部分噪聲的同時，一定程度上模糊了圖像。綜合比較去噪效果和模糊損失，效果最好的復(fù)原方法是空域高斯平滑模板，因為噪聲在空域添加，而高斯模板在空域又具備高斯函數(shù)的起伏特征，起到了保留大部分原始像素值和抑制大部分噪聲值得作用。這些復(fù)原手段嚴(yán)格來說是不理想的，理想效果應(yīng)是能近似得到原始圖像，即圖3(a)所示圖像。2 圖像的運動模糊2.1 模糊函數(shù)的推導(dǎo)圖像的運動模糊是由于在成像曝光過程中，物體和相機之間的相對運動，使得物體上同一點的光線在多個成像單元上引起響應(yīng)造成的。由于相機的成像單元在成像期間完成一個積分過程，以表示模糊后的圖像，T表示積分時間，運動模糊原理可以表示為對上式進行傅里葉變換變換，得利用傅里葉變換的移位性質(zhì)，頻域模糊函數(shù)可以表示為現(xiàn)在考慮最簡單的情況，那就是如果物體在兩個坐標(biāo)軸方向的運動是勻速直線運動，那么最終的模糊函數(shù)可以寫為此過程的建模結(jié)果與頻域濾波形式相類似，在進行處理時可以參照頻域濾波的步驟，以作為濾波器函數(shù)，對進行填充和中心化后的圖像頻譜進行直接相乘即可。需要注意的一點就是在生成時也要按照頻譜中心化的形式進行，即的大小應(yīng)是原始圖像尺寸的兩倍，且默認(rèn)中心為坐標(biāo)原點。2.2 實驗結(jié)果在openCV3.1.0計算機視覺庫基礎(chǔ)上實現(xiàn)圖像的運動模糊退化。首先編寫函數(shù)motionBlur(省略形參)，以原圖像為參考，生成一個兩倍大小的運動模糊濾波器矩陣，包括兩個部分：實部和虛部。然后按照頻域濾波步驟進行后續(xù)操作，此時濾波過程是復(fù)數(shù)運算，應(yīng)注意兩個通道的不同實虛對應(yīng)關(guān)系的相乘相加，此時需要借助函數(shù)multiply()和addWeighted()來實現(xiàn)3。具體源程序見附錄源代碼2，不同參數(shù)時的模糊結(jié)果如下圖所示。abcdef圖4 圖像的運動模糊退化，上圖從左至右、由上及下一次為(a)-(f)；（a）原始圖像；（b）參數(shù)a=0.01,b=0.01,T=1時的運動模糊退化圖像；（c）參數(shù)a=0.05,b=0.05,T=1時的運動模糊退化圖像；（d）參數(shù)a=0.1,b=0.1,T=1時的運動模糊退化圖像；（e）參數(shù)a=0，b=0.1,T=1時的運動模糊退化圖像；（f）參數(shù)a=0.05,b=0.05,T=0.5時的運動模糊退化圖像2.3 結(jié)果分析觀察上圖實驗結(jié)果可以發(fā)現(xiàn)，本文較好的實現(xiàn)了圖像的運動模糊。從不同參數(shù)的吃力結(jié)果可知，參數(shù)a，b分別決定了物體在x，y兩個方向的運動速度，反映在圖像中就是圖像在兩個方向上的模糊程度，a和b越大，圖像的模糊程度越劇烈。對比圖3的(c)和(f)可知，兩幅圖像的參數(shù)a，b完全相同，但(f)的曝光時間T要比(c)小?？梢钥闯銎毓鈺r間T越小，意味著成像單元的積分時間越短，那么相同反射強度的物體成像的亮度越暗，這就是(f)的亮度明顯小于(c)的原因。因此，直觀來說，積分時間決定了圖像整體的亮度大小。當(dāng)T小于1時，圖像整體亮度會低于原圖；當(dāng)T大于1時，圖像整體亮度會高于原圖。在進行運動模糊時，圖像的邊緣出現(xiàn)了一定寬度的“暗帶”，其寬度隨運動速度的增大而增大，當(dāng)運動速度較小時基本可以忽略，而當(dāng)運動速度較大時，這種現(xiàn)象尤其明顯。這是因為在整個積分時間內(nèi)，這些區(qū)域只有很短的非零時間，大多是時間被積量為零，所以最終的灰度很小。這種現(xiàn)象在實際攝影成像中也會出現(xiàn)，因此不需要進行更正。2.4 為模糊圖像添加噪聲在圖4(c)的基礎(chǔ)上結(jié)合1.2中方法，為運動模糊后的圖像添加高斯噪聲，此處選取運動模糊參數(shù)：a=0.05,b=0.05,T=1。由于兩個步驟均已實現(xiàn)，此處將不再附錄代碼，結(jié)果見圖5abcd圖5 為運動模糊后的圖像添加高斯噪聲；（a）運動模糊退化后的圖像；（b）添加均值為零，方差為10的高斯噪聲；（c）添加均值為10，方差為10的高斯噪聲；（d）添加均值為20，方差為10的高斯噪聲 3 維納濾波復(fù)原圖像3.1 維納濾波原理維納濾波是一種基于最小均方誤差準(zhǔn)則、對平穩(wěn)過程的最優(yōu)估計器4。這種濾波器的輸出與期望輸出之間的均方誤差為最小，因此，它是一個最佳濾波系統(tǒng)。它可用于提取被平穩(wěn)噪聲所污染的信號。維納濾波的誤差度量由下式給出帶入圖像參量及退化函數(shù)，求解最小均方誤差，得到最終的濾波器頻域表達式其中常數(shù)K表示噪聲功率譜與為退化圖像功率譜比，即3.2 濾波復(fù)原結(jié)果維納濾波實現(xiàn)原理同2中的運動模糊，基本思路都是圖像的頻域濾波，編寫函數(shù)wienerFilter()構(gòu)造維納濾波器的實部和虛部，函數(shù)代碼見附錄代碼3，復(fù)原結(jié)果如下圖所示。圖6 維納濾波復(fù)原結(jié)果，上圖從左至右依次是(a)-(c)；(a)未退化的圖像；(b)運動模糊(a=0.05,b=0.05,T=1)并添加高斯噪聲(均值位0，方差為10)的退化圖像；(c)維納濾波復(fù)原結(jié)果3.2 結(jié)果分析對比以上三幅圖像可知，維納濾波能顯著去除運動模糊，但是在未退化圖像的功率譜未知的情況下，對噪聲的去除效果并不理想，最終恢復(fù)出的圖像與未退化圖像相比仍有較大差距。因此在實際操作中，應(yīng)考慮交互式地選取參數(shù)K，以取得最佳的恢復(fù)效果。4 約束最小二乘方濾波4.1 濾波原理維納濾波是建立在最小化統(tǒng)計準(zhǔn)則的基礎(chǔ)上，因此在平均意義上它是最優(yōu)的，而約束最小二乘方濾波對每一幅圖像都能取得最優(yōu)的效果。一般地，約束最小二乘方濾波器的函數(shù)表達式為其中表示圖像的退化函數(shù)，是一個參數(shù)，為取得最佳濾波效果，其調(diào)整目標(biāo)是使得是函數(shù)的傅里葉變換，是一個時域的模板：4.2 濾波復(fù)原效果實現(xiàn)方法同3.2，同樣都是頻域濾波操作的基本步驟。對P進行傅里葉變換之前要將其填充到原始圖像的二倍大小；為簡單記，參數(shù)通過交互形式指定，實驗結(jié)果如下圖所示。圖7 約束最小二乘方濾波濾波結(jié)果，上圖從左至右依次是(a)-(c)；(a)未退化的圖像；(b)運動模糊(a=0.05,b=0.05,T=1)并添加高斯噪聲(均值位0，方差為10)的退化圖像；(c)約束最小二乘方復(fù)原結(jié)果4.3 結(jié)果分析與維納濾波不同，最小二乘方濾波對噪聲的去除效果很好，最終的復(fù)原圖像中噪聲的影響已經(jīng)非常微弱，在這一點上其性能優(yōu)于維納濾波，這通過比較圖6（c）和圖7（c）能夠明顯的體會出。但是最終的復(fù)原圖像依然存在明顯的模糊現(xiàn)象，最小二乘方濾波并不能有效恢復(fù)運動模糊（至少在本文所做實驗中是這樣）。附錄參考文獻1 Rafael C Gonzelez，Richard E Woods. 數(shù)字圖像處理. 阮秋琦，阮宇智，等/譯.北京：電子工業(yè)出版社. 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.Castleman. 數(shù)字圖像處理. 朱志剛，林學(xué)訚，石定機，等譯. 北京：/電子工業(yè)出版社. 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法生成滿足高斯分布的隨機數(shù)；para mean=均值，sigma=標(biāo)準(zhǔn)差double generateGaussianNoise(double mean, double sigma)/ 定義小值const double epsilon = std:numeric_limits:min();static double z0, z1;double u1, u2;/ 生成滿足0-1均勻分布的隨機數(shù)dou1 = rand() * (1.0 / RAND_MAX);u2 = rand() * (1.0 / RAND_MAX); while (u1 = epsilon);/ flag為真構(gòu)造高斯隨機變量Xz0 = sqrt(-2.0 * log(u1) * cos(2 * CV_PI * u2);z1 = sqrt(-2.0 * log(u1) * sin(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 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;/ 定義高斯平滑函數(shù),besarKernel:模板邊長void myGaussianBlur(const Mat& src, Mat& result, int besarKernel, double sigma)/ 高斯模板半徑int kerR = besarKernel / 2;/ 定義高斯模板因子cv:Mat kernel = Mat_(besarKernel, besarKernel);/ 歸一化參數(shù)double alpha = 1 / (2 * 3.1416 * sigma * sigma);double sum = 0;/ 模板函數(shù)生成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 = 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.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, 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 medianFilterImage, 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);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:waitKey(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.141592using namespace cv;using namespace std;/ 根據(jù)圖像生成退化函數(shù)void motionBlur(cv:Mat srcImage, cv:Mat rePart, cv:Mat imPart, double a, double 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=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: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 * 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_32FC1);shiftToMiddle(paddingImage);/ 為傅立葉變換的結(jié)果(實部和虛部)分配存儲空間cv:Mat planes = cv:Mat_(paddingImage),cv:Mat:zeros(paddingImage.size(), CV_32F) ;Mat completeI;/ 為延擴后的圖像增添一個初始化為0的通道m(xù)erge(planes, 2, completeI);/ 進行離散傅立葉變換dft(completeI, completeI);split(completeI, planes);/相乘cv:Mat rePart1(m, n, CV_32FC1);cv:Mat rePart2(m, n, CV_32FC1);cv:Mat imPart1(m, n, CV_32FC1);cv:Mat imPart2(m, n, CV_32FC1);/生成結(jié)果實部multiply(planes0, rePart, rePart1);multiply(planes1, imPart, rePart2);addWeighted(rePart1, 1, rePart2, -1, 0, planes0);/生成結(jié)果虛部multiply(planes0, imPart, imPart1);multiply(planes1, rePart, imPart2);addWeighted(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);/拷貝結(jié)果到原圖像大小畫布cv:Mat resltImage(srcImage.rows, 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);return 0;3. 維納濾波復(fù)原（openCV3.1.0 C+）（基本結(jié)構(gòu)同2，此處省略主函數(shù)，僅提供維納濾波模板的生成函數(shù)）void wienerFilter(cv:Mat srcImage, c