Matlab
main函数
clc
clear
Ii=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
clear
Ii=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;
}