Matlab

main函數

clc clearIi=imread('Ii','bmp'); psf=imread('psf','bmp');recovery=ima_restoration(Ii,psf,0.1); figure;imagesc(recovery);axis square;

ima_restoration函數

function [ recovery ] = ima_restoration(Ii,psf,a)H=newfft2(psf); H=H./max(max(H)); Gi=newfft2(Ii); recovery=newifft2(Gi./(H+a)); recovery=abs(recovery); end

newfft2(f)函數

function u=newfft2(f) f=fftshift(f); u=fft2(f); u=fftshift(u);

newifft2(f)函數

function u=newifft2(f) f=ifftshift(f); u=ifft2(f); u=ifftshift(u);

1. 合并MATLAB程序

clc clearIi=imread('Ii','bmp'); psf=imread('psf','bmp'); ps=fftshift(psf); M=fft2(ps); H=fftshift(M); S1=max(H); S = max(S1); Hh=H ./ S; pf=fftshift(Ii); Mm=fft2(pf); Gi=fftshift(Mm); a=0.1; Hhh=Hh+a; GG=Gi./Hhh; ff=ifftshift(GG); uu=ifft2(ff); imshow(uu); recoveryy=ifftshift(uu); recovery=abs(recoveryy); figure; imagesc(recovery); axis square;

2. C++實現功能

1.main.cpp

#include"ImaRestoration.h"int main(int argc, char* argv[]) {Mat Ii = imread("Ii.bmp", 0);Mat psf = imread("psf.bmp", 0);double a = 0.1;cv::Mat color = ImaRestoration(Ii, psf, a);imshow("color", color);waitKey(0);return 0; }

2.ImaRestoration.h

ImaRestoration函數三個參數，返回一個參數
color = ImaRestoration(Ii, psf, a);

function [ recovery ] = ima_restoration(Ii,psf,a)H=newfft2(psf); H=H./max(max(H)); Gi=newfft2(Ii); recovery=newifft2(Gi./(H+a)); recovery=abs(recovery); end

ImaRestoration函數分解成：

ps=fftshift(psf); M=fft2(ps); H=fftshift(M); S1=max(H); S = max(S1); Hh=H ./ S; pf=fftshift(Ii); Mm=fft2(pf); Gi=fftshift(Mm); a=0.1; Hhh=Hh+a; GG=Gi./Hhh; ff=ifftshift(GG); uu=ifft2(ff); imshow(uu); recoveryy=ifftshift(uu); recovery=abs(recoveryy); figure; imagesc(recovery); axis square;

matlab程序中ima_restoration 主要是先對psf.bmp進行newfft2變換得到H，找到H中像素值最大的值進行歸一化處理得到H，之后對Ii.bmp進行newfft2變換得到Gi，然后，對Gi點除(H+a)后的結果進行newifft2處理，最后得到變換后的幅值圖像輸出。

使用opencv來實現：

//ImaRestoration函數 Mat ImaRestoration(cv::Mat Ii, cv::Mat psf, double a = 0.1) { Mat ps = fftshift(psf);Mat M = fft2(ps);Mat H = fftshift(M);//獲取最大值像素double minValue = -1;double maxValue = -1;cv::Vec2d comp;cv::Vec2d* HPtr = H.ptr<cv::Vec2d>(0);for (int i = 0; i < H.size().area(); i++){double val = sqrtf(HPtr[i][0] * HPtr[i][0] + HPtr[i][1] * HPtr[i][1]);if (val > maxValue){maxValue = val;comp = HPtr[i];}}Mat Hh = H.clone();int rowNumber = Hh.rows;int colNumber = Hh.cols;for (int i = 0; i < rowNumber; i++){cv::Vec2d*data = Hh.ptr<cv::Vec2d>(i);for (int j = 0; j < colNumber; j++){double a = data[j][0];double b = data[j][1];double c = comp[0];double d = comp[1];data[j][0] = (a*c + b * d) / (c*c + d * d);data[j][1] = (b*c - a * d) / (c*c + d * d);}}Mat pf = fftshift(Ii);Mat Mm = fft2(pf);Mat Gi = fftshift(Mm);Mat Hhh = Hh + a;Mat GG(Gi.size(), Gi.type());for (int i = 0; i < rowNumber; i++){cv::Vec2d*data = Hhh.ptr<cv::Vec2d>(i);cv::Vec2d*dataG = Gi.ptr<cv::Vec2d>(i);cv::Vec2d*dataGG = GG.ptr<cv::Vec2d>(i);for (int j = 0; j < colNumber; j++){double a = dataG[j][0];double b = dataG[j][1];double c = data[j][0];double d = data[j][1];dataGG[j][0] = (a*c + b * d) / (c*c + d * d);dataGG[j][1] = (b*c - a * d) / (c*c + d * d);}}Mat ff = ifftshift(GG);Mat uu = ifft2(ff);Mat recovery = ifftshift(uu);double* imptr = recovery.ptr<double>(0);int pixSZ = recovery.size().area();double maxV = std::numeric_limits<double>::min();double minV = std::numeric_limits<double>::max();int normMin = 0;int normMax = 63;for (int p = 0; p < pixSZ; p++){if (imptr[p] > maxV) maxV = imptr[p];if (imptr[p] < minV) minV = imptr[p];}double k = (double)(normMax - normMin) / (maxV - minV);double b = normMin - k * minV;cv::Mat3b color(recovery.size());cv::Vec3b* cPtr = color.ptr<cv::Vec3b>(0);for (int p = 0; p < pixSZ; p++){int index = (int)(k* (imptr[p]) + b);cPtr[p][0] = lutArr[index * 3];cPtr[p][1] = lutArr[index * 3 + 1];cPtr[p][2] = lutArr[index * 3 + 2];}return color; }

循環位移函數circshift

Mat circshift(Mat out, const Point delta) {Size sz = out.size();assert(sz.height > 0 && sz.width > 0);if ((sz.height == 1 && sz.width == 1) || (delta.x == 0 && delta.y == 0))return Mat();int x = delta.x;int y = delta.y;if (x > 0) x = x % sz.width;if (y > 0) y = y % sz.height;if (x < 0) x = x % sz.width + sz.width;if (y < 0) y = y % sz.height + sz.height;vector<Mat> planes;split(out, planes);for (size_t i = 0; i < planes.size(); i++){// 豎直方向移動Mat tmp0, tmp1, tmp2, tmp3;Mat q0(planes[i], Rect(0, 0, sz.width, sz.height - y));Mat q1(planes[i], Rect(0, sz.height - y, sz.width, y));q0.copyTo(tmp0);q1.copyTo(tmp1);tmp0.copyTo(planes[i](Rect(0, y, sz.width, sz.height - y)));tmp1.copyTo(planes[i](Rect(0, 0, sz.width, y)));// 水平方向移動Mat q2(planes[i], Rect(0, 0, sz.width - x, sz.height));Mat q3(planes[i], Rect(sz.width - x, 0, x, sz.height));q2.copyTo(tmp2);q3.copyTo(tmp3);tmp2.copyTo(planes[i](Rect(x, 0, sz.width - x, sz.height)));tmp3.copyTo(planes[i](Rect(0, 0, x, sz.height)));}merge(planes, out);return out; }

fftshift函數

Mat fftshift(Mat Src) {Size sz = Src.size();Point pt(0, 0);pt.x = (int)floor(sz.width / 2.0);pt.y = (int)floor(sz.height / 2.0);Mat out = circshift(Src, pt);return out; }

fft2函數

Mat fft2(const Mat src) {int mat_type = src.type();Mat Fourier;assert(mat_type < 15); //不支持的數據格式if (mat_type < 7){Mat planes[] = { Mat_<double>(src), Mat::zeros(src.size(), CV_64F) };merge(planes, 2, Fourier);dft(Fourier, Fourier);}else // 7<mat_type<15{Mat tmp;dft(src, tmp);vector<Mat> planes;split(tmp, planes);magnitude(planes[0], planes[1], planes[0]); //將復數轉化為幅值Fourier = planes[0];}return Fourier; }

ifft2函數

Mat ifft2(const Mat src) {Mat Fourier;int mat_type = src.type();assert(mat_type < 15); //不支持的數據格式if (mat_type < 7){Mat planes[] = { Mat_<double>(src), Mat::zeros(src.size(), CV_64F) };merge(planes, 2, Fourier);dft(Fourier, Fourier, DFT_INVERSE + DFT_SCALE, 0);}else // 7<mat_type<15{Mat tem;dft(src, tem, DFT_INVERSE + DFT_SCALE, 0);vector<Mat> planes;split(tem, planes);magnitude(planes[0], planes[1], planes[0]); //將復數轉化為幅值Fourier = planes[0];}return Fourier; }

ifftshift函數

Mat ifftshift(Mat Scr) {Size sz = Scr.size();Point pt(0, 0);pt.x = (int)ceil(sz.width / 2.0);pt.y = (int)ceil(sz.height / 2.0);Mat out = circshift(Scr, pt);return out; }

完整的ImaRestoration程序：ImaRestoration.h

#ifndef IMA_RESTORATION_H #define IMA_RESTORATION_H #include<opencv2/opencv.hpp> #include <opencv2/core/core.hpp> #include <opencv2/highgui/highgui.hpp> #include <opencv2/imgproc/imgproc.hpp> #include <math.h> #include <stdio.h> using namespace cv; using namespace std;int lutArr[192] = {135, 42, 53,147, 48, 54,160, 55, 54,173, 61, 53,186, 67, 50,199, 74, 44,212, 83, 32,221, 92, 15,225, 99, 3,225, 104, 2,224, 109, 4,222, 113, 8,220, 117, 13,218, 121, 16,216, 125, 18,214, 129, 20,212, 133, 20,211, 137, 19,210, 142, 16,210, 147, 12,209, 152, 9,207, 156, 7,205, 160, 6,202, 164, 6,198, 167, 6,194, 169, 7,190, 172, 10,185, 174, 15,180, 177, 21,175, 179, 29,169, 181, 37,164, 183, 46,158, 185, 56,152, 187, 66,146, 188, 77,140, 189, 89,134, 190, 101,128, 191, 113,123, 191, 124,119, 191, 135,115, 191, 146,111, 191, 156,107, 190, 165,103, 190, 174,100, 189, 183,96, 188, 192,93, 188, 200,89, 187, 209,86, 186, 217,82, 185, 225,78, 185, 233,74, 185, 241,68, 187, 248,61, 190, 253,55, 195, 255,50, 200, 254,46, 206, 252,42, 211, 250,38, 216, 247,33, 222, 245,29, 228, 245,24, 235, 245,19, 243, 246,14, 251, 249 };Mat circshift(Mat out, const Point delta) {Size sz = out.size();assert(sz.height > 0 && sz.width > 0);if ((sz.height == 1 && sz.width == 1) || (delta.x == 0 && delta.y == 0))return Mat();int x = delta.x;int y = delta.y;if (x > 0) x = x % sz.width;if (y > 0) y = y % sz.height;if (x < 0) x = x % sz.width + sz.width;if (y < 0) y = y % sz.height + sz.height;vector<Mat> planes;split(out, planes);for (size_t i = 0; i < planes.size(); i++){// 豎直方向移動Mat tmp0, tmp1, tmp2, tmp3;Mat q0(planes[i], Rect(0, 0, sz.width, sz.height - y));Mat q1(planes[i], Rect(0, sz.height - y, sz.width, y));q0.copyTo(tmp0);q1.copyTo(tmp1);tmp0.copyTo(planes[i](Rect(0, y, sz.width, sz.height - y)));tmp1.copyTo(planes[i](Rect(0, 0, sz.width, y)));// 水平方向移動Mat q2(planes[i], Rect(0, 0, sz.width - x, sz.height));Mat q3(planes[i], Rect(sz.width - x, 0, x, sz.height));q2.copyTo(tmp2);q3.copyTo(tmp3);tmp2.copyTo(planes[i](Rect(x, 0, sz.width - x, sz.height)));tmp3.copyTo(planes[i](Rect(0, 0, x, sz.height)));}merge(planes, out);return out; }Mat fftshift(Mat Src) {Size sz = Src.size();Point pt(0, 0);pt.x = (int)floor(sz.width / 2.0);pt.y = (int)floor(sz.height / 2.0);Mat out = circshift(Src, pt);return out; }Mat fft2(const Mat src) {int mat_type = src.type();Mat Fourier;assert(mat_type < 15); //不支持的數據格式if (mat_type < 7){Mat planes[] = { Mat_<double>(src), Mat::zeros(src.size(), CV_64F) };merge(planes, 2, Fourier);dft(Fourier, Fourier);}else // 7<mat_type<15{Mat tmp;dft(src, tmp);vector<Mat> planes;split(tmp, planes);magnitude(planes[0], planes[1], planes[0]); //將復數轉化為幅值Fourier = planes[0];}return Fourier; }Mat ifft2(const Mat src) {Mat Fourier;int mat_type = src.type();assert(mat_type < 15); //不支持的數據格式if (mat_type < 7){Mat planes[] = { Mat_<double>(src), Mat::zeros(src.size(), CV_64F) };merge(planes, 2, Fourier);dft(Fourier, Fourier, DFT_INVERSE + DFT_SCALE, 0);}else // 7<mat_type<15{Mat tem;dft(src, tem, DFT_INVERSE + DFT_SCALE, 0);vector<Mat> planes;split(tem, planes);magnitude(planes[0], planes[1], planes[0]); //將復數轉化為幅值Fourier = planes[0];}return Fourier; }Mat ifftshift(Mat Scr) {Size sz = Scr.size();Point pt(0, 0);pt.x = (int)ceil(sz.width / 2.0);pt.y = (int)ceil(sz.height / 2.0);Mat out = circshift(Scr, pt);return out; }Mat ImaRestoration(cv::Mat Ii, cv::Mat psf, double a = 0.1) {Mat ps = fftshift(psf);Mat M = fft2(ps);Mat H = fftshift(M);//獲取最大值像素double minValue = -1;double maxValue = -1;cv::Vec2d comp;cv::Vec2d* HPtr = H.ptr<cv::Vec2d>(0);for (int i = 0; i < H.size().area(); i++){double val = sqrtf(HPtr[i][0] * HPtr[i][0] + HPtr[i][1] * HPtr[i][1]);if (val > maxValue){maxValue = val;comp = HPtr[i];}}Mat Hh = H.clone();int rowNumber = Hh.rows;int colNumber = Hh.cols;for (int i = 0; i < rowNumber; i++){cv::Vec2d*data = Hh.ptr<cv::Vec2d>(i);for (int j = 0; j < colNumber; j++){double a = data[j][0];double b = data[j][1];double c = comp[0];double d = comp[1];data[j][0] = (a*c + b * d) / (c*c + d * d);data[j][1] = (b*c - a * d) / (c*c + d * d);}}Mat pf = fftshift(Ii);Mat Mm = fft2(pf);Mat Gi = fftshift(Mm);Mat Hhh = Hh + a;Mat GG(Gi.size(), Gi.type());for (int i = 0; i < rowNumber; i++){cv::Vec2d*data = Hhh.ptr<cv::Vec2d>(i);cv::Vec2d*dataG = Gi.ptr<cv::Vec2d>(i);cv::Vec2d*dataGG = GG.ptr<cv::Vec2d>(i);for (int j = 0; j < colNumber; j++){double a = dataG[j][0];double b = dataG[j][1];double c = data[j][0];double d = data[j][1];dataGG[j][0] = (a*c + b * d) / (c*c + d * d);dataGG[j][1] = (b*c - a * d) / (c*c + d * d);}}Mat ff = ifftshift(GG);Mat uu = ifft2(ff);Mat recovery = ifftshift(uu);double* imptr = recovery.ptr<double>(0);int pixSZ = recovery.size().area();double maxV = std::numeric_limits<double>::min();double minV = std::numeric_limits<double>::max();int normMin = 0;int normMax = 63;for (int p = 0; p < pixSZ; p++){if (imptr[p] > maxV) maxV = imptr[p];if (imptr[p] < minV) minV = imptr[p];}double k = (double)(normMax - normMin) / (maxV - minV);double b = normMin - k * minV;cv::Mat3b color(recovery.size());cv::Vec3b* cPtr = color.ptr<cv::Vec3b>(0);for (int p = 0; p < pixSZ; p++){int index = (int)(k* (imptr[p]) + b);cPtr[p][0] = lutArr[index * 3];cPtr[p][1] = lutArr[index * 3 + 1];cPtr[p][2] = lutArr[index * 3 + 2];}return color; } #endif

圖像復原——維納濾波
main.cpp

#include"ImaRestoration.h"int main(int argc, char* argv[]) {Mat Ii = imread("Ii.bmp", 0);Mat psf = imread("psf.bmp", 0);double a = 0.1;cv::Mat color = ImaRestoration(Ii, psf, a);imshow("color", color);waitKey(0);return 0; }

總結

以上是生活随笔為你收集整理的图像复原——维纳滤波的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。