一、高斯滤波器英文介绍：https://en.wikipedia.org/wiki/Gaussian_filter
相关博客：http://www.cnblogs.com/wangguchangqing/p/6407717.html
下面是整合的代码实现：
//高斯滤波器
#include<opencv2/opencv.hpp>
#include<opencv2/imgproc/imgproc.hpp>
#include<iostream>
using namespace cv;
using namespace std;
void GaussianFilter(const Mat &src, Mat &dst, int ksize, double sigma)
{
CV_Assert(src.channels() || src.channels() == 3); // 只处理单通道或者三通道图像
const static double pi = 3.1415926;
// 根据窗口大小和sigma生成高斯滤波器模板
// 申请一个二维数组，存放生成的高斯模板矩阵
double **templateMatrix = new double*[ksize];
for (int i = 0; i < ksize; i++)
templateMatrix[i] = new double[ksize];
int origin = ksize / 2; // 以模板的中心为原点
double x2, y2;
double sum = 0;
for (int i = 0; i < ksize; i++)
{
double x2 = pow(double(i-origin), 2);
for (int j = 0; j < ksize; j++)
{
double y2 = pow(double(j-origin), 2);
// 高斯函数前的常数可以不用计算，会在归一化的过程中给消去
double g = exp(-(x2 + y2) / (2 * sigma * sigma));
sum += g;
templateMatrix[i][j] = g;
}
}
for (int i = 0; i < ksize; i++)
{
for (int j = 0; j < ksize; j++)
{
templateMatrix[i][j] /= sum;
cout << templateMatrix[i][j] << " ";
}
cout << endl;
}
// 将模板应用到图像中
int border;
border= ksize/2;
//BorderTypes没有包含进来
int borderType=BORDER_DEFAULT;
copyMakeBorder(src, dst, border, border, border, border,borderType);
int channels = dst.channels();
int rows = dst.rows - border;
int cols = dst.cols - border;
//i相当于是遍历原图像的行，所以是从border到原图像的行-border;
//j相当于是遍历原图像的列
for (int i = border; i < rows; i++)
{
for (int j = border; j < cols; j++)
{
double sum[3] = { 0 };
//以（i,j）位置的像素为中心，a表示行距离中心点的距离，b表示列距离中心点的距离
//此处实现模板有对应像素相乘
for (int a = -border; a <= border; a++)
{
for (int b = -border; b <= border; b++)
{
if (channels == 1)
{
sum[0] += templateMatrix[border + a][border + b] * dst.at<uchar>(i + a, j + b);
}
else if (channels == 3)
{
Vec3b rgb = dst.at<Vec3b>(i + a, j + b);
auto k = templateMatrix[border + a][border + b];
sum[0] += k * rgb[0];
sum[1] += k * rgb[1];
sum[2] += k * rgb[2];
}
}
}
for (int k = 0; k < channels; k++)
{
if (sum[k] < 0)
sum[k] = 0;
else if (sum[k] > 255)
sum[k] = 255;
}
if (channels == 1)
dst.at<uchar>(i, j) = static_cast<uchar>(sum[0]);
else if (channels == 3)
{
//一下注释的这种方法是不可行的
//Vec3b rgb = (static_cast<uchar>(sum[0]), static_cast<uchar>(sum[1]), static_cast<uchar>(sum[2]));
Vec3b rgb1;
rgb1[0]=static_cast<uchar>(sum[0]);
rgb1[1]=static_cast<uchar>(sum[1]);
rgb1[2]=static_cast<uchar>(sum[2]);
dst.at<Vec3b>(i, j) = rgb1;
}
}
}
// 释放模板数组
for (int i = 0; i < ksize; i++)
delete[] templateMatrix[i];
delete[] templateMatrix;
}
void separateGaussianFilter(const Mat &src, Mat &dst, int ksize, double sigma)
{
CV_Assert(src.channels()==1 || src.channels() == 3); // 只处理单通道或者三通道图像
// 生成一维的高斯滤波模板
double *matrix = new double[ksize];
double sum = 0;
int origin = ksize / 2;
for (int i = 0; i < ksize; i++)
{
// 高斯函数前的常数可以不用计算，会在归一化的过程中给消去
double g = exp(double(-(i - origin) * (i - origin) / (2 * sigma * sigma)));
sum += g;
matrix[i] = g;
}
// 归一化
for (int i = 0; i < ksize; i++)
matrix[i] /= sum;
// 将模板应用到图像中
int border = ksize / 2;
int borderType=BORDER_REFLECT;
copyMakeBorder(src, dst, border, border, border, border, borderType);
int channels = dst.channels();
int rows = dst.rows - border;
int cols = dst.cols - border;
// 水平方向
for (int i = border; i < rows; i++)
{
for (int j = border; j < cols; j++)
{
double sum[3] = { 0 };
for (int k = -border; k <= border; k++)
{
if (channels == 1)
{
sum[0] += matrix[border + k] * dst.at<uchar>(i, j + k); // 行不变，列变化；先做水平方向的卷积
}
else if (channels == 3)
{
Vec3b rgb = dst.at<Vec3b>(i, j + k);
sum[0] += matrix[border + k] * rgb[0];
sum[1] += matrix[border + k] * rgb[1];
sum[2] += matrix[border + k] * rgb[2];
}
}
for (int k = 0; k < channels; k++)
{
if (sum[k] < 0)
sum[k] = 0;
else if (sum[k] > 255)
sum[k] = 255;
}
if (channels == 1)
dst.at<uchar>(i, j) = static_cast<uchar>(sum[0]);
else if (channels == 3)
{
//Vec3b rgb = ( static_cast<uchar>(sum[0]), static_cast<uchar>(sum[1]), static_cast<uchar>(sum[2]) );
Vec3b rgb1;
rgb1[0]=static_cast<uchar>(sum[0]);
rgb1[1]=static_cast<uchar>(sum[1]);
rgb1[2]=static_cast<uchar>(sum[2]);
dst.at<Vec3b>(i, j) = rgb1;
}
}
}
// 竖直方向
for (int i = border; i < rows; i++)
{
for (int j = border; j < cols; j++)
{
double sum[3] = { 0 };
for (int k = -border; k <= border; k++)
{
if (channels == 1)
{
sum[0] += matrix[border + k] * dst.at<uchar>(i + k, j); // 列不变，行变化；竖直方向的卷积
}
else if (channels == 3)
{
Vec3b rgb = dst.at<Vec3b>(i + k, j);
sum[0] += matrix[border + k] * rgb[0];
sum[1] += matrix[border + k] * rgb[1];
sum[2] += matrix[border + k] * rgb[2];
}
}
for (int k = 0; k < channels; k++)
{
if (sum[k] < 0)
sum[k] = 0;
else if (sum[k] > 255)
sum[k] = 255;
}
if (channels == 1)
dst.at<uchar>(i, j) = static_cast<uchar>(sum[0]);
else if (channels == 3)
{
Vec3b rgb1;
rgb1[0]=static_cast<uchar>(sum[0]);
rgb1[1]=static_cast<uchar>(sum[1]);
rgb1[2]=static_cast<uchar>(sum[2]);
dst.at<Vec3b>(i, j) = rgb1;
}
}
}
delete[] matrix;
}
int main()
{
Mat img=imread("lena.jpg",1);//第二个参数大于0，,返回三通道图像，=0，返回灰度图像
Mat dstimg;
Mat dstimg1;
Mat dstimg2;
int k=5;
double sigma=0.8;
GaussianFilter(img, dstimg, k,sigma);
separateGaussianFilter(img, dstimg2, k,sigma);
GaussianBlur(img,dstimg1,Size(5,5),sigma,sigma,BORDER_DEFAULT );
imshow("src image",img);
imshow("dest image",dstimg);
imshow("dest image1",dstimg1);
imshow("dest image2",dstimg2);
waitKey();
return 0;
}运行结果:
二、双边滤波器
相关博客：
http://blog.csdn.net/abcjennifer/article/details/7616663（含matlab实现）
http://www.cnblogs.com/wangguchangqing/p/6416401.html（含c++,opencv实现）
下面，是对大神写的代码的运行：
//双边滤波器
#include <opencv2/core/core.hpp>
#include<opencv2/opencv.hpp>
#include<opencv2/imgproc/imgproc.hpp>
#include<opencv2/highgui/highgui.hpp>
#include<stdlib.h>
#include<stdio.h>
#include<vector>
#include<iostream>
using namespace cv;
using namespace std;
void myBilateralFilter(const Mat &src, Mat &dst, int ksize, double space_sigma, double color_sigma)
{
int channels = src.channels();
CV_Assert(channels == 1 || channels == 3);
double space_coeff = -0.5 / (space_sigma * space_sigma);
double color_coeff = -0.5 / (color_sigma * color_sigma);
int radius = ksize / 2;
Mat temp;
int borderType=BORDER_REFLECT;
copyMakeBorder(src, temp, radius, radius, radius, radius, borderType);
vector<double> _color_weight(channels * 256); // 存放差值的平方
vector<double> _space_weight(ksize * ksize); // 空间模板系数
vector<int> _space_ofs(ksize * ksize); // 模板窗口的坐标
double *color_weight = &_color_weight[0];
double *space_weight = &_space_weight[0];
int *space_ofs = &_space_ofs[0];
for (int i = 0; i < channels * 256; i++)
color_weight[i] = exp(i * i * color_coeff);
// 生成空间模板
int maxk = 0;
for (int i = -radius; i <= radius; i++)
{
for (int j = -radius; j <= radius; j++)
{
double r =sqrt(double (i*i +j*j));
if (r > radius)
continue;
space_weight[maxk] = exp(r * r * space_coeff); // 存放模板系数
space_ofs[maxk++] = i * temp.step + j * channels; // 存放模板的位置，和模板系数相对应
}
}
// 滤波过程
for (int i = 0; i < src.rows; i++)
{
const uchar *sptr = temp.data + (i + radius) * temp.step + radius * channels;
uchar *dptr = dst.data + i * dst.step;//局部变量显示时，此处便是错误指针了
if (channels == 1)
{
for (int j = 0; j < src.cols; j++)
{
double sum = 0, wsum = 0;
int val0 = sptr[j]; // 模板中心位置的像素
for (int k = 0; k < maxk; k++)
{
int val = sptr[j + space_ofs[k]];
double w = space_weight[k] * color_weight[abs(val - val0)]; // 模板系数 = 空间系数 * 灰度值系数
sum += val * w;
wsum += w;
}
dptr[j] = (uchar)cvRound(sum / wsum);
}
}
if (channels == 3)
{
for (int j = 0; j < src.cols*3; j += 3)
{
double sum_b = 0, sum_g = 0, sum_r = 0, wsum = 0;
int b0 = sptr[j];
int g0 = sptr[j + 1];
int r0 = sptr[j + 2];
for (int k = 0; k < maxk; k++)
{
const uchar *sptr_k = sptr + j + space_ofs[k];
int b = sptr_k[0];
int g = sptr_k[1];
int r = sptr_k[2];
double w = space_weight[k] * color_weight[abs(b - b0) + abs(g - g0) + abs(r - r0)];
sum_b += b * w;
sum_g += g * w;
sum_r += r * w;
wsum += w;
}
wsum = 1.0f / wsum;
b0 = cvRound(sum_b * wsum);
g0 = cvRound(sum_g * wsum);
r0 = cvRound(sum_r * wsum);
dptr[j] = (uchar)b0;//这里出现写入指针错误，为什么？
dptr[j + 1] = (uchar)g0;
dptr[j + 2] = (uchar)r0;
}
}
}
}
int main()
{
Mat srcimg=imread("lena.jpg");
cout<<srcimg.cols<<srcimg.rows<<endl;
Mat dstimg(srcimg.rows, srcimg.cols, CV_8UC3);
int k=20;
double space_sigma=80;
double color_sigma=80;
myBilateralFilter(srcimg,dstimg,k,space_sigma,color_sigma);
imshow("src image",srcimg);
imshow("dest image",dstimg);
waitKey();
return 0;
}
三、引导滤波器
四、Gabor滤波
英文介绍：https://en.wikipedia.org/wiki/Gabor_filter
大都根据以下公式得到：
1.维基百科上面给出的示例matlab代码，做了少许调整，分别画出滤波器组的实部和幅值图
u=5;
v=8;
m=39;
n=39;
%以上参数是自己设置的
gaborArray = cell(u,v);
fmax = 0.25;%最大频率
gama = sqrt(2);
eta = sqrt(2);%高斯分布的标准差
for i = 1:u
fu = fmax/((sqrt(2))^(i-1));%正弦波的频率，频率与波长成反比
alpha = fu/gama;
beta = fu/eta;
for j = 1:v
tetav = ((j-1)/v)*pi;%gabor滤波器的方向，v个方向，一般为8，分别为0,1/8*pi......
x = 1:m;
y = 1:n;
[X,Y]=meshgrid(x,y);
X=X-(m+1)/2;
Y=Y-(m+1)/2;
xprime = X.*cos(tetav)+Y.*sin(tetav);
yprime = -X.*sin(tetav)+Y.*cos(tetav);
gFilter = (fu^2/(pi*gama*eta)).*exp(-((alpha^2).*(xprime.^2)+(beta^2).*(yprime.^2))).*exp(1i*2*pi*fu.*xprime);%alpha和beta不明白为什么是这样的，无法理解和公式的匹配
% gFilter = exp(-((1/(2*eta^2)).*(xprime.^2)+(gama/(2*eta^2)).*(yprime.^2))).*exp(1i*2*pi*fu.*xprime);%这样子的话其实尺度没有发生什么变化
gaborArray{i,j} = gFilter;
end
end
figure('NumberTitle','Off','Name','Magnitudes of Gabor filters');
for i = 1:u
for j = 1:v
subplot(u,v,(i-1)*v+j);
mesh(abs( gaborArray{i,j}));
end
end
figure('NumberTitle','Off','Name','real part of Gabor filters');
for i = 1:u
for j = 1:v
subplot(u,v,(i-1)*v+j);
mesh(real( gaborArray{i,j}));
end
end
figure('NumberTitle','Off','Name','Magnitudes of Gabor filters');
for i = 1:u
for j = 1:v
subplot(u,v,(i-1)*v+j);
imshow(abs(gaborArray{i,j}),[]);
end
end
figure('NumberTitle','Off','Name','Real parts of Gabor filters');
for i = 1:u
for j = 1:v
subplot(u,v,(i-1)*v+j);
imshow(real(gaborArray{i,j}),[]);
end
endgabor特征的提取函数
function featureVector = gaborFeatures(img,gaborArray,d1,d2)
% GABORFEATURES extracts the Gabor features of an input image.
% It creates a column vector, consisting of the Gabor features of the input
% image. The feature vectors are normalized to zero mean and unit variance.
%
%
% Inputs:
% img : Matrix of the input image
% gaborArray : Gabor filters bank created by the function gaborFilterBank
% d1 : The factor of downsampling along rows.
% d2 : The factor of downsampling along columns.
%
% Output:
% featureVector : A column vector with length (m*n*u*v)/(d1*d2).
% This vector is the Gabor feature vector of an
% m by n image. u is the number of scales and
% v is the number of orientations in 'gaborArray'.
%
%
% Sample use:
%
% img = imread('cameraman.tif');
% gaborArray = gaborFilterBank(5,8,39,39); % Generates the Gabor filter bank
% featureVector = gaborFeatures(img,gaborArray,4,4); % Extracts Gabor feature vector, 'featureVector', from the image, 'img'.
%
%
%
% Details can be found in:
%
% M. Haghighat, S. Zonouz, M. Abdel-Mottaleb, "CloudID: Trustworthy
% cloud-based and cross-enterprise biometric identification,"
% Expert Systems with Applications, vol. 42, no. 21, pp. 7905-7916, 2015.
%
%
%
% (C) Mohammad Haghighat, University of Miami
% haghighat@ieee.org
% PLEASE CITE THE ABOVE PAPER IF YOU USE THIS CODE.
if (nargin ~= 4) % Check correct number of arguments
error('Please use the correct number of input arguments!')
end
if size(img,3) == 3 % Check if the input image is grayscale
warning('The input RGB image is converted to grayscale!')
img = rgb2gray(img);
end
img = double(img);
%% Filter the image using the Gabor filter bank
% Filter input image by each Gabor filter
[u,v] = size(gaborArray);
gaborResult = cell(u,v);
for i = 1:u
for j = 1:v
gaborResult{i,j} = imfilter(img, gaborArray{i,j});
end
end
%% Create feature vector
% Extract feature vector from input image
featureVector = [];
for i = 1:u
for j = 1:v
gaborAbs = abs(gaborResult{i,j});
gaborAbs = downsample(gaborAbs,d1);
gaborAbs = downsample(gaborAbs.',d2);
gaborAbs = gaborAbs(:);
% Normalized to zero mean and unit variance. (if not applicable, please comment this line)
gaborAbs = (gaborAbs-mean(gaborAbs))/std(gaborAbs,1);
featureVector = [featureVector; gaborAbs];
end
end
%% Show filtered images (Please comment this section if not needed!)
% Show real parts of Gabor-filtered images
figure('NumberTitle','Off','Name','Real parts of Gabor filters');
for i = 1:u
for j = 1:v
subplot(u,v,(i-1)*v+j)
imshow(real(gaborResult{i,j}),[]);
end
end
% Show magnitudes of Gabor-filtered images
figure('NumberTitle','Off','Name','Magnitudes of Gabor filters');
for i = 1:u
for j = 1:v
subplot(u,v,(i-1)*v+j)
imshow(abs(gaborResult{i,j}),[]);
end
endfunction gb=gabor_fn(sigma,theta,lambda,psi,gamma)
sigma_x = sigma;
sigma_y = sigma/gamma;
% Bounding box
nstds = 3;
xmax = max(abs(nstds*sigma_x*cos(theta)),abs(nstds*sigma_y*sin(theta)));
xmax = ceil(max(1,xmax));
ymax = max(abs(nstds*sigma_x*sin(theta)),abs(nstds*sigma_y*cos(theta)));
ymax = ceil(max(1,ymax));
xmin = -xmax; ymin = -ymax;
[x,y] = meshgrid(xmin:xmax,ymin:ymax);
% Rotation
x_theta=x*cos(theta)+y*sin(theta);
y_theta=-x*sin(theta)+y*cos(theta);
gb= exp(-.5*(x_theta.^2/sigma_x^2+y_theta.^2/sigma_y^2)).*cos(2*pi/lambda*x_theta+psi);#include <opencv2/core/core.hpp>
#include <opencv2/imgproc/imgproc.hpp>
#include <opencv2/highgui/highgui.hpp>
#include <cmath>
#include <iostream>
using namespace cv;
using namespace std;
const double PI = 3.14159265;
// ref: http://blog.csdn.net/watkinsong/article/details/7876361
Mat getMyGabor(int width, int height, int U, int V, double Kmax, double f,
double sigma, int ktype, const string& kernel_name)
{
//CV_ASSERT(width % 2 == 0 && height % 2 == 0);
//CV_ASSERT(ktype == CV_32F || ktype == CV_64F);
int half_width = width / 2;
int half_height = height / 2;
double Qu = PI*U/8;
double sqsigma = sigma*sigma;
double Kv = Kmax/pow(f,V);
double postmean = exp(-sqsigma/2);
Mat kernel_re(width, height, ktype);
Mat kernel_im(width, height, ktype);
Mat kernel_mag(width, height, ktype);
double tmp1, tmp2, tmp3;
for(int j = -half_height; j <= half_height; j++){
for(int i = -half_width; i <= half_width; i++){
tmp1 = exp(-(Kv*Kv*(j*j+i*i))/(2*sqsigma));
tmp2 = cos(Kv*cos(Qu)*i + Kv*sin(Qu)*j) - postmean;
tmp3 = sin(Kv*cos(Qu)*i + Kv*sin(Qu)*j);
if(ktype == CV_32F)
kernel_re.at<float>(j+half_height, i+half_width) =
(float)(Kv*Kv*tmp1*tmp2/sqsigma);
else
kernel_re.at<double>(j+half_height, i+half_width) =
(double)(Kv*Kv*tmp1*tmp2/sqsigma);
if(ktype == CV_32F)
kernel_im.at<float>(j+half_height, i+half_width) =
(float)(Kv*Kv*tmp1*tmp3/sqsigma);
else
kernel_im.at<double>(j+half_height, i+half_width) =
(double)(Kv*Kv*tmp1*tmp3/sqsigma);
}
}
magnitude(kernel_re, kernel_im, kernel_mag);
if(kernel_name.compare("real") == 0)
return kernel_re;
else if(kernel_name.compare("imag") == 0)
return kernel_im;
else if(kernel_name.compare("mag") == 0)
return kernel_mag;
else
printf("Invalid kernel name!\n");
}
void construct_gabor_bank()
{
double Kmax = PI/2;
double f = sqrt(2.0);
double sigma = 2*PI;
int U = 7;
int V = 4;
int GaborH = 129;
int GaborW = 129;
Mat kernel;
Mat totalMat;
for(U = 0; U < 8; U++){
Mat colMat;
for(V = 0; V < 5; V++){
kernel = getMyGabor(GaborW, GaborH, U, V,
Kmax, f, sigma, CV_64F, "real");
//show gabor kernel
normalize(kernel, kernel, 0, 1, CV_MINMAX);
printf("U%dV%d\n", U, V);
//imshow("gabor", kernel);
//waitKey(0);
if(V == 0)
colMat = kernel;
else
vconcat(colMat, kernel, colMat);
}
if(U == 0)
totalMat = colMat;
else
hconcat(totalMat, colMat, totalMat);
}
imshow("gabor bank", totalMat);
//imwrite("gabor_bank.jpg",totalMat);
waitKey(0);
}
Mat gabor_filter(Mat& img)
{
// variables for gabor filter
double Kmax = PI/2;
double f = sqrt(2.0);
double sigma = 2*PI;
int U = 7;
int V = 4;
int GaborH = 129;
int GaborW = 129;
//
Mat kernel_re, kernel_im;
Mat dst_re, dst_im, dst_mag;
// variables for filter2D
Point archor(-1,-1);
int ddepth = -1;
double delta = 0;
// filter image with gabor bank
Mat totalMat;
for(U = 0; U < 8; U++){
Mat colMat;
for(V = 0; V < 5; V++){
kernel_re = getMyGabor(GaborW, GaborH, U, V,
Kmax, f, sigma, CV_64F, "real");
kernel_im = getMyGabor(GaborW, GaborH, U, V,
Kmax, f, sigma, CV_64F, "imag");
filter2D(img, dst_re, ddepth, kernel_re);
filter2D(img, dst_im, ddepth, kernel_im);
dst_mag.create(img.rows, img.cols, CV_32FC1);
magnitude(Mat_<float>(dst_re),Mat_<float>(dst_im),
dst_mag);
//show gabor kernel
normalize(dst_mag, dst_mag, 0, 1, CV_MINMAX);
printf("U%dV%d\n", U, V);
//imshow("dst_mag", dst_mag);
//waitKey(0);
if(V == 0)
colMat = dst_mag;
else
vconcat(colMat, dst_mag, colMat);
}
if(U == 0)
totalMat = colMat;
else
hconcat(totalMat, colMat, totalMat);
}
return totalMat;
}
int main( int argc, char** argv )
{
//construct_gabor_bank();
Mat image;
image = imread("aa.jpg",0); // Read the file
if(! image.data ) // Check for invalid input
{
cout << "Could not open or find the image" << std::endl ;
return -1;
}
Mat filterd_image = gabor_filter(image);
imshow("filtered image", filterd_image);
//imwrite("filterd_image.jpg",filterd_image);
waitKey(0);
return 0;
}
3389
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
