圖像復原的方法很多，常用的比較經典的是反向濾波法和約束還原法。博主在做反向濾波實驗的過程中，發現圖像復原的關鍵是退化模型的建立，可以誇張地說：要有好的復原效果就得根據各個圖像的退化特點建立相關的退化模型，並在退化模型的基礎上做相關的濾波或者說對待處理的像素做相應的處理，從而盡可能地復原圖像。再說一遍，復原方法的關鍵是退化模型。可以想到的是，由於造成圖像退化的原因五花八門。簡單的有加性退化、減性、乘性、除性退化等；復雜的有非線性退化等。從這一點看來似乎沒有通用的復原方法，這樣一來似乎只能使用深度學習等智能算法做一些通用復原算法的研究了。博主推薦在使用神經網絡等算法進行研究之前，先修《數據分析與數據挖掘》和《復雜性思維》，或許能從中找到通用復原算法的鑰匙！
實驗內容
利用逆濾波和其他逆卷積算法對運動模糊或散焦模糊圖像進行圖像復原，並給出實現結果。
【背景知識】
1. 圖像退化模型
圖形復原處理的關鍵是建立退化模型，原圖像f(x,y)是通過一個系統H及加入一個外來加性噪聲n(x,y)而退化成一幅圖像g(x,y)的。
這樣圖像的退化過程的數學表達式就可以寫為g(x,y)=H[f(x,y)]+n(x,y)。容易想到圖像退化及復原的過程如下：
如果系統H滿足下面兩個式子，那么系統就是線性和空間位置不變的系統。在圖像復原處理中，非線性和空間變化的系統的模型雖然更具普遍性和准確性，但它卻給處理工作帶來巨大的困難，它常常沒有解或很難用計算機來處理。實際的成像系統在一定條件下往往可以近似地視為線性和空間不變的系統，因此在圖像復原處理中，往往用線性和空間不變的系統模型加以近似地視為線性和空間不變的系統。
2.連續的退化模型
線性系統H可由其沖激響應來表征，當系統H空間位置不變時，則 h(x-α,y-β)=H[δ(x-α,y-β)]。系統H對輸入f(x,y)的響應就是系統輸入信號f(x,y)與系統沖激響應的卷積。考慮加性噪聲n(x,y)時。可寫為：
其中n(x,y)與圖像中位置無關。
3.離散的退化模型
如果給出 A×B 大小的數字圖像，以及 C×D 大小的點擴散函數，可首先擴展大小為 M×N 的周期延拓圖像。
為避免折疊，要求 M≥A+C-1，N≥B+D-1。這樣一來，fe(x,y)和he(x,y)分別成為二維周期函數，它們在x和y方向上的周期分別為M和N。由此得到二維退化模型為一個二維卷積形式：
式中，x=0,1,2,3，…，M-1；y=0,1,2，…，N-1；ge(x,y)也為周期函數，其周期與fe(x,y)和he(x,y)完全一樣。上式也可用矩陣表示為 g=Hf+n。G(u,v)、F(u,v)和N(u,v)分別是g(x,y)、f(x,y)和n(x,y)的二維傅里葉變換。於是有G(u,v)=MNH(u,v)F(u,v)+N(u,v)。這樣就將求f(x,y)的過程轉換為求解F(u,v)的過程，簡化了計算過程，同時上式也是進行圖像復原的基礎。
圖像的復原方法：
反向濾波法
約束還原法
運動模糊圖像的復原
模糊模型
水平勻速直線運動引起模糊的復原
復原的關鍵在退化模型的估計，復原的基礎建立在退化模型之上。
逆濾波代碼如下：
#include <opencv2/highgui.hpp> #include <opencv2/imgproc.hpp> #include <iostream> #define GRAY_THRESH 10 using namespace std; using namespace cv; void idfft(Mat complexI, Mat image, int oph, int opw, string str); void zero_to_center(Mat& freq_plane); int main() { char imgname[] = "1.png"; Mat img = imread(imgname, 0); if (img.empty()) return 0; imshow("srcimg", img); string filename = "1.png"; Mat image = imread(filename, IMREAD_GRAYSCALE); // 灰度圖像 if (image.empty()) return 0; imshow("src", image); image.convertTo(image, CV_32FC1); // 轉換到32位浮點型單通道矩陣 ///快速傅里葉變換/ int oph = getOptimalDFTSize(image.rows); int opw = getOptimalDFTSize(image.cols); Mat padded; copyMakeBorder(image, padded, 0, oph - image.rows, 0, opw - image.cols, BORDER_CONSTANT, Scalar::all(0)); // 增加邊框 Mat pad_show; normalize(padded, pad_show, 0, 1, NORM_MINMAX); //歸一化 imshow("padded", pad_show); // 顯示添加了邊框的原圖 Mat temp[] = { padded, Mat::zeros(padded.size(),CV_32FC1) }; // 二通道 Mat complexI; merge(temp, 2, complexI); // 合並多個陣列以形成單個多通道陣列 dft(complexI, complexI); // 傅里葉變換 //顯示頻譜圖 split(complexI, temp); zero_to_center(temp[0]); zero_to_center(temp[1]); Mat aa; magnitude(temp[0], temp[1], aa); // 計算計算二維矢量的幅值。dst(I)=sqrt(x(I)^2+y(I)^2) divide(aa, oph * opw, aa); // 除法 imshow("pu", aa); merge(temp, 2, complexI); /*----------------------------------退化-------------------------------------------------------*/ // 生成退化函數 const float ee = 1e-3; Mat Degenerate_function(padded.size(), CV_32FC2); // 32位浮點型雙通道 for (int i = 0; i < Degenerate_function.rows; i++) { float* dq = Degenerate_function.ptr<float>(i); for (int j = 0; j < Degenerate_function.cols; j++) { double hmi = -pow(pow(i - oph / 2, 2) + pow(j - opw / 2, 2), 5.0 / 6); float mh = exp(hmi) < ee ? ee : exp(hmi); dq[2 * j] = mh; dq[2 * j + 1] = mh; } } multiply(complexI, Degenerate_function, complexI); // 計算兩個數組的每元素縮放乘積。 Mat ccomplexI = complexI.clone(); idfft(ccomplexI, image, oph, opw, "退化"); /*----------------------------------運算部分---------------------------------------------------*/ // 頻域濾波 //生成頻域濾波核 巴特沃斯低通濾波器 （注：不使用濾波器效果更佳） Mat butter_sharpen(padded.size(), CV_32FC2); // 32位浮點型雙通道 double D0 = 50.0; int n = 4; for (int i = 0; i < butter_sharpen.rows; i++) { float* q = butter_sharpen.ptr<float>(i); for (int j = 0; j < butter_sharpen.cols; j++) { double d = sqrt(pow(i - oph / 2, 2) + pow(j - opw / 2, 2)); q[2 * j] = 1.0 / (1 + pow(D0 / d, 2 * n)); q[2 * j + 1] = 1.0 / (1 + pow(D0 / d, 2 * n)); } } // 生成退化函數 Mat Degenerate_function1(padded.size(), CV_32FC2); // 32位浮點型雙通道 float MM = oph * opw; for (int i = 0; i < Degenerate_function1.rows; i++) { float* dq = Degenerate_function1.ptr<float>(i); for (int j = 0; j < Degenerate_function1.cols; j++) { double hmi = -pow(pow(i - oph / 2, 2) + pow(j - opw / 2, 2), 5.0 / 6); float mh = exp(hmi) * MM < ee ? ee : exp(hmi) * MM; dq[2 * j] = mh; dq[2 * j + 1] = mh; } } multiply(complexI, butter_sharpen, complexI); // 計算兩個數組的每元素縮放乘積。 divide(complexI, Degenerate_function1, complexI); // 除法 idfft(complexI, image, oph, opw, "復原"); waitKey(); return 0; } void idfft(Mat complexI, Mat image, int oph, int opw, string str) { string s1 = "_filter", s2 = "dstSharpen"; //傅里葉反變換 idft(complexI, complexI, DFT_INVERSE); Mat dstSharpen[2]; split(complexI, dstSharpen); // 將多通道陣列划分為多個單通道陣列。 // magnitude(dstSharpen[0],dstSharpen[1],dstSharpen[0]); //求幅值(模) for (int i = 0; i < oph; i++) { float* q = dstSharpen[0].ptr<float>(i); for (int j = 0; j < opw; j++) { q[j] = q[j] * pow(-1, i + j); } } Mat show; // divide(dstSharpen[0], dstSharpen[0].cols*dstSharpen[0].rows, show); normalize(dstSharpen[0], show, 0, 1, NORM_MINMAX); //.... show = show(Rect(0, 0, image.cols, image.rows)); imshow(str+s1, show); threshold(dstSharpen[0], dstSharpen[0], 0, 255, THRESH_BINARY); // 閾值化 normalize(dstSharpen[0], dstSharpen[0], 0, 1, NORM_MINMAX); // 映射到0~1 dstSharpen[0] = dstSharpen[0](Rect(0, 0, image.cols, image.rows)); imshow(str+s2, dstSharpen[0]); } void zero_to_center(Mat& freq_plane) { // freq_plane = freq_plane(Rect(0, 0, freq_plane.cols & -2, freq_plane.rows & -2)); //這里為什么&上-2具體查看opencv文檔 //其實是為了把行和列變成偶數 -2的二進制是11111111.......10 最后一位是0 int cx = freq_plane.cols / 2; int cy = freq_plane.rows / 2; //以下的操作是移動圖像 (零頻移到中心) 與函數center_transform()作用相同，只是一個先處理，一個dft后再變換 Mat part1_r(freq_plane, Rect(0, 0, cx, cy)); //元素坐標表示為(cx,cy) Mat part2_r(freq_plane, Rect(cx, 0, cx, cy)); Mat part3_r(freq_plane, Rect(0, cy, cx, cy)); Mat part4_r(freq_plane, Rect(cx, cy, cx, cy)); Mat tmp; part1_r.copyTo(tmp); //左上與右下交換位置(實部) part4_r.copyTo(part1_r); tmp.copyTo(part4_r); part2_r.copyTo(tmp); //右上與左下交換位置(實部) part3_r.copyTo(part2_r); tmp.copyTo(part3_r); }
#include <iostream> #include <opencv2/opencv.hpp> #include <cmath> using namespace cv; using namespace std; void calcPSF(Mat& outputImg, Size filterSize, int len, double theta); // 維納濾波 void fftshift(const Mat& inputImg, Mat& outputImg); // 頻譜平移 void filter2DFreq(const Mat& inputImg, Mat& outputImg, const Mat& H); // 卷積 void calcWnrFilter(const Mat& input_h_PSF, Mat& output_G, double nsr); // 傅里葉變換 void edgetaper(const Mat& inputImg, Mat& outputImg, double gamma = 5.0, double beta = 0.2); // 減少圖像恢復時的振鈴效果 int LEN = 50; int THETA = 360; int snr = 8000; Mat imgIn; Rect roi; Mat himage; int main() { string strFileName = "1.png"; Mat imgIn1 = imread(strFileName, IMREAD_GRAYSCALE); imshow("原始圖像", imgIn1); string strInFileName = "2.png"; Mat imgIn = imread(strInFileName, IMREAD_GRAYSCALE); imshow("退化圖像", imgIn); Rect roi = Rect(0, 0, imgIn.cols & -2, imgIn.rows & -2); // ***&-2指保留最大的2次冪 ***&11111110 int D0 = 30; // 截止頻率 int n = 1; // 巴特沃斯濾波系數的分母冪次 float k = 0.0025; // 退化函數的冪系數 0.001 0.00025 int rlen = imgIn.rows & -2, clen = imgIn.cols & -2; float hf = 1.0 / (rlen * rlen); Mat HimgIn(roi.size(), CV_32F, Scalar(0)); // 巴特沃斯低通濾波器 Mat Himage(roi.size(), CV_32F, Scalar(0)); // 點擴散函數 for (int i = 0; i < rlen; i++) { float* himg = HimgIn.ptr<float>(i); float* hima = Himage.ptr<float>(i); for (int j = 0; j < clen; j++) { // 巴特沃斯 float D = sqrt((i - rlen / 2) * (i - rlen / 2) + (j - clen / 2) * (j - clen / 2)); float hfen = pow((1 + D / D0), 2 * n); hfen >= 0 ? hfen : 1e-3; // hfen >= 255 ? 255 : hfen; himg[j] = 1.0 / hfen; // 退化 float hmi = -k * pow(((i - rlen / 2) * (i - rlen / 2) + (j - clen / 2) * (j - clen / 2)), 5.0 / 6); float hhh = 1.0 / exp(hmi) * hf; hima[j] = hhh; } } himage = Himage; cout << rlen << "\t" << clen << endl; Mat outimg, fhimg; Mat Gimg = imgIn.clone(); calcWnrFilter(HimgIn, fhimg, 1.0 / double(snr)); Gimg.convertTo(Gimg, CV_32F); edgetaper(Gimg, Gimg); filter2DFreq(Gimg(roi), outimg, fhimg); outimg.convertTo(outimg, CV_8U); normalize(outimg, outimg, 0, 255, NORM_MINMAX); imshow("result", outimg); // normalize(fGimg, fGimg, 0, 255, NORM_MINMAX); // imshow("頻譜", fGimg); // cout << pow((1 + 1), 5.0 / 6) << endl; Mat imgOut; // cout << (3 & -2) << endl; // 維納濾波開始 /*Mat Hw, h; calcPSF(h, roi.size(), LEN, THETA); calcWnrFilter(h, Hw, 1.0 / double(snr)); imgIn.convertTo(imgIn, CV_32F); edgetaper(imgIn, imgIn); filter2DFreq(imgIn(roi), imgOut, Hw); imgOut.convertTo(imgOut, CV_8U); normalize(imgOut, imgOut, 0, 255, NORM_MINMAX); imshow("result", imgOut);*/ waitKey(); return 0; } //! [calcPSF] void calcPSF(Mat& outputImg, Size filterSize, int len, double theta) { Mat h(filterSize, CV_32F, Scalar(0)); Point point(filterSize.width / 2, filterSize.height / 2); ellipse(h, point, Size(0, cvRound(float(len) / 2.0)), 90.0 - theta, 0, 360, Scalar(255), FILLED); Scalar summa = sum(h); outputImg = h / summa[0]; } //! [calcPSF] //! [fftshift] void fftshift(const Mat& inputImg, Mat& outputImg) { outputImg = inputImg.clone(); int cx = outputImg.cols / 2; int cy = outputImg.rows / 2; Mat q0(outputImg, Rect(0, 0, cx, cy)); Mat q1(outputImg, Rect(cx, 0, cx, cy)); Mat q2(outputImg, Rect(0, cy, cx, cy)); Mat q3(outputImg, Rect(cx, cy, cx, cy)); Mat tmp; q0.copyTo(tmp); q3.copyTo(q0); tmp.copyTo(q3); q1.copyTo(tmp); q2.copyTo(q1); tmp.copyTo(q2); } //! [fftshift] //! [filter2DFreq] void filter2DFreq(const Mat& inputImg, Mat& outputImg, const Mat& H) { Mat planes[2] = { Mat_<float>(inputImg.clone()), Mat::zeros(inputImg.size(), CV_32F) }; Mat complexI; merge(planes, 2, complexI); dft(complexI, complexI, DFT_SCALE); // 二維的 Mat planesH[2] = { Mat_<float>(H.clone()), Mat::zeros(H.size(), CV_32F) }; Mat complexH; merge(planesH, 2, complexH); Mat complexIH; mulSpectrums(complexI, complexH, complexIH, 0); Mat planesHH[2] = { Mat_<float>(himage.clone()),Mat::zeros(himage.size(),CV_32F) }; Mat complexHH; merge(planesHH, 2, complexHH); Mat complexIHH; mulSpectrums(complexIH, complexHH, complexIHH, 0); complexIHH += Scalar::all(1); // switch to logarithmic scale log(complexIHH, complexIHH); log(complexIHH, complexIHH); idft(complexIHH, complexIHH); split(complexIHH, planes); outputImg = planes[0]; /*idft(complexIH, complexIH); split(complexIH, planes); outputImg = planes[0];*/ } //! [filter2DFreq] //! [calcWnrFilter] void calcWnrFilter(const Mat& input_h_PSF, Mat& output_G, double nsr) { Mat h_PSF_shifted; fftshift(input_h_PSF, h_PSF_shifted); Mat planes[2] = { Mat_<float>(h_PSF_shifted.clone()), Mat::zeros(h_PSF_shifted.size(), CV_32F) }; Mat complexI; merge(planes, 2, complexI); dft(complexI, complexI); split(complexI, planes); Mat denom; pow(abs(planes[0]), 2, denom); denom += nsr; divide(planes[0], denom, output_G); } //! [calcWnrFilter] //! [edgetaper] void edgetaper(const Mat& inputImg, Mat& outputImg, double gamma, double beta) { int Nx = inputImg.cols; int Ny = inputImg.rows; Mat w1(1, Nx, CV_32F, Scalar(0)); Mat w2(Ny, 1, CV_32F, Scalar(0)); float* p1 = w1.ptr<float>(0); float* p2 = w2.ptr<float>(0); float dx = float(2.0 * CV_PI / Nx); float x = float(-CV_PI); for (int i = 0; i < Nx; i++) { p1[i] = float(0.5 * (tanh((x + gamma / 2) / beta) - tanh((x - gamma / 2) / beta))); x += dx; } float dy = float(2.0 * CV_PI / Ny); float y = float(-CV_PI); for (int i = 0; i < Ny; i++) { p2[i] = float(0.5 * (tanh((y + gamma / 2) / beta) - tanh((y - gamma / 2) / beta))); y += dy; } Mat w = w2 * w1; multiply(inputImg, w, outputImg); } //! [edgetaper]
-----------------------------------continue-------------------------------------------