---
### 要求:
1.去噪
2.提取脊线、汗孔处理
---
### 思路过程:
**1.滤波**
先判断是什么类型的噪声,再进行滤波处理。
- 灰度直方图可以判断出椒盐噪声,指数噪声,均匀噪声但需要一块平坦区域,由于指纹的图片没有平坦区域,因此效果不好。
- 傅里叶变换可以判断出周期性噪声。先试试看。做完傅里叶变换后发现我看不出来。
- 直接用滤波试一试:均值,高斯,中值,双边,感觉中值的效果好一点。
**2.提取脊线**
就是对图像里剩余的一些小白点,小黑点进行消除,采用了腐蚀膨胀,再对连通域进行处理的方法
---
### 遇到的错误,经验
1
**问题:** 读入图片报错
```
OpenCV Error: Assertion failed
```
**解决办法:** Mat srcImage = imread(" ../../im/27.bmp "); 冒号里接空格导致,删去即可。
2
**问题:** 双边滤波too few arguments to function ‘void cv::bilateralFilter
解决办法:==未解决哭了!==
3
**问题:**
```
Process finished with exit code 139 (interrupted by signal 11: SIGSEGV)
```
**解决办法:** 由于Mat A,只声明,没有初始化导致
---
### 代码文件
代码放在有道云里了,但是作为刚刚开始学编程的人,代码写的惨不忍睹: [http://note.youdao.com/noteshare?id=f5e38dbf56f27f36af93d23ba6fea77c](http://note.youdao.com/noteshare?id=f5e38dbf56f27f36af93d23ba6fea77c)
---
### 效果图
图片看不了。。。点这个链接[http://note.youdao.com/noteshare?id=2f34834905c45b74d2d63fbe8b09002b](http://note.youdao.com/noteshare?id=2f34834905c45b74d2d63fbe8b09002b)
---
### 具体内容:数字图像处理第四、五章涉及的内容基本都涉及了一下: 傅里叶变换,频率域滤波,空间滤波。
**1.灰度直方图**
```
//灰度直方图
//opencv有计算灰度直方图的函数calcHist,但用起来不太方便,封装一下,直接用封装好的就可以了。
class Histogram1D
{
private:
int histSize[1]; // 项的数量
float hranges[2]; // 统计像素的最大值和最小值
const float* ranges[1];
int channels[1]; // 仅计算一个通道
public:
Histogram1D()
{
// 准备1D直方图的参数
histSize[0] = 256;
hranges[0] = 0.0f;
hranges[1] = 255.0f;
ranges[0] = hranges;
channels[0] = 0;
}
MatND getHistogram(const Mat &image)
{
MatND hist;
// 计算直方图
calcHist(&image ,// 要计算图像的
1, // 只计算一幅图像的直方图
channels, // 通道数量
Mat(), // 不使用掩码
hist, // 存放直方图
1, // 1D直方图
histSize, // 统计的灰度的个数
ranges); // 灰度值的范围
return hist;
}
Mat getHistogramImage(const Mat &image)
{
MatND hist = getHistogram(image);
// 最大值,最小值
double maxVal = 0.0f;
double minVal = 0.0f;
minMaxLoc(hist, &minVal, &maxVal);
//显示直方图的图像
Mat histImg(histSize[0], histSize[0], CV_8U, Scalar(255));
// 设置最高点为nbins的90%
int hpt = static_cast<int>(0.9 * histSize[0]);
//每个条目绘制一条垂直线
for (int h = 0; h < histSize[0]; h++)
{
float binVal = hist.at<float>(h);
int intensity = static_cast<int>(binVal * hpt / maxVal);
// 两点之间绘制一条直线
line(histImg, Point(h, histSize[0]), Point(h, histSize[0] - intensity), Scalar::all(0));
}
return histImg;
}
};
//这么调用
Histogram1D hist;
Mat histImg;
histImg = hist.getHistogramImage(srcImage);
imshow("Histogram", histImg);
waitKey(0);
```
**2.傅里叶变换**
```
//傅里叶变换
//将输入图像拓展到最佳的尺寸,边界使用0进行补充
int m = getOptimalDFTSize(srcImage.rows);
int n = getOptimalDFTSize(srcImage.cols);
//将添加的像素初始化为0
Mat padded;
copyMakeBorder(srcImage, padded, m - srcImage.rows, 0, n - srcImage.cols, 0, BORDER_CONSTANT, Scalar::all(0));
//为傅里叶变换的结果分配存储空间(分为实部和虚部)
//将planes数组合并成为一个多通道的数组complexI
Mat planes[] = { Mat_<float>(padded), Mat::zeros(padded.size(), CV_32F) };
Mat complexI;
merge(planes, 2, complexI);
//就地进行离散傅里叶变换
dft(complexI, complexI);
//将负数转换为幅值,即为:log(1+sqrt(Re(DFT(I))^2+Im(DFT(I))^2))
split(complexI, planes); //将多通道数组分为两个单通道数组
//planes[0] = Re(DFT(I)),planes[1] = Im(DFT(I))
magnitude(planes[0], planes[1], planes[0]);
Mat MagnitudeImage = planes[0];
//进行对数尺度放缩
MagnitudeImage += Scalar::all(1);
log(MagnitudeImage, MagnitudeImage); //求自然对数
//剪切和重新分布幅度图的象限
//若有奇数行和奇数列,进行频谱剪裁
MagnitudeImage = MagnitudeImage(Rect(0, 0, MagnitudeImage.cols &-2, MagnitudeImage.rows &-2));
//这句中的 & -2 是位操作,因为-2的二进制为除最后一位为0外其余均为1,因此与之后的结果为将最低位的1与掉
//实现将原来的奇数变为偶数
//重新排列傅里叶图像中的象限,使得原点位原图像中心
int cx = MagnitudeImage.cols / 2;
int cy = MagnitudeImage.rows / 2;
Mat q0(MagnitudeImage, Rect(0, 0, cx, cy)); //ROI区域的左上方
Mat q1(MagnitudeImage, Rect(cx, 0, cx, cy)); //ROI区域的右上方
Mat q2(MagnitudeImage, Rect(0, cy, cx, cy)); //ROI区域的左下方
Mat q3(MagnitudeImage, Rect(cx, cy, cx, cy)); //ROI区域的右下方
//交换象限,左上与右下交换
Mat tmp;
q0.copyTo(tmp);
q3.copyTo(q0);
tmp.copyTo(q3);
//交换象限,右上与左下进行交换
q1.copyTo(tmp);
q2.copyTo(q1);
tmp.copyTo(q2);
//归一化,用0~1之间的浮点数将矩阵变换为可视图的格式
normalize(MagnitudeImage, MagnitudeImage, 0, 1, CV_MINMAX);
//显示效果图
imshow("27fourier", MagnitudeImage);
imwrite("../../result/27fourier.bmp", MagnitudeImage);
waitKey();
```
**3.陷波滤波**(这里没用上,因为在频谱图上并看不出哪里是噪声点。)
```
#include <iostream>
#include <opencv2/opencv.hpp>
using namespace std;
using namespace cv;
//扩充边界
Mat image_add_border( Mat &src )
{
int w=2*src.cols;
int h=2*src.rows;
std::cout << "src: " << src.cols << "*" << src.rows << std::endl;
cv::Mat padded; //panding填充的意思,一般都用来扩充边界
copyMakeBorder( src, padded, 0, h-src.rows, 0, w-src.cols,cv::BORDER_CONSTANT, cv::Scalar::all(0));
//copyMakerBorder用来扩充边界的函数,中间的四个参数表示上下左右的扩充大小。
//BORDER_CONSTANT以一个常量像素值(由参数 value给定)填充扩充的边界值,这里用Scalar::all(0)黑色来赋值
padded.convertTo(padded,CV_32FC1);
std::cout << "dst: " << padded.cols << "*" << padded.rows << std::endl;
return padded;
}
//transform to center 中心化
//????只是给坐标和为奇数的加了一个负号,怎么就中心化了
void center_transform( cv::Mat &src )
{
for(int i=0; i<src.rows; i++){
float *p = src.ptr<float>(i); //p是指向每一行的指针
for(int j=0; j<src.cols; j++){
p[j] = p[j] * pow(-1, i+j);//pow(a,b),求a的b次幂
}
}
}
//对角线交换内容
void zero_to_center(cv::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;//以下的操作是移动图像 (零频移到中心)
cv::Mat part1_r(freq_plane, cv::Rect(0,0,cx,cy)); //元素坐标表示为(cx,cy)
cv::Mat part2_r(freq_plane, cv::Rect(cx,0,cx,cy));
cv::Mat part3_r(freq_plane, cv::Rect(0,cy,cx,cy));
cv::Mat part4_r(freq_plane, cv::Rect(cx,cy,cx,cy));
cv::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);
}
void show_spectrum( cv::Mat &complexI )
{
cv::Mat temp[] = {cv::Mat::zeros(complexI.size(),CV_32FC1),cv::Mat::zeros(complexI.size(),CV_32FC1)};
//显示频谱图
cv::split(complexI, temp);
cv::Mat aa;
cv::magnitude(temp[0], temp[1], aa);
// zero_to_center(aa);
cv::divide(aa, aa.cols*aa.rows, aa);
double min, max;
cv::Point min_loc, max_loc;
cv::minMaxLoc(aa, &min, &max, &min_loc, &max_loc );
std::cout << "min: " << min << " max: " << max << std::endl;
std::cout << "min_loc: " << min_loc << " max_loc: " << max_loc << std::endl;
cv::circle( aa, min_loc, 20, cv::Scalar::all(1), 3);
cv::circle( aa, max_loc, 20, cv::Scalar::all(1), 3);
cv::imshow("src_img_spectrum",aa);
}
//频率域滤波
cv::Mat frequency_filter(cv::Mat &padded,cv::Mat &blur)
{
cv::Mat plane[]={padded, cv::Mat::zeros(padded.size(), CV_32FC1)};
cv::Mat complexIm;
cv::merge(plane,2,complexIm);
cv::dft(complexIm,complexIm);//fourior transform
show_spectrum(complexIm);
cv::multiply(complexIm, blur, complexIm);
cv::idft(complexIm, complexIm, CV_DXT_INVERSE); //idft
cv::Mat dst_plane[2];
cv::split(complexIm, dst_plane);
center_transform(dst_plane[0]);
// center_transform(dst_plane[1]);
cv::magnitude(dst_plane[0],dst_plane[1],dst_plane[0]); //求幅值(模)
// center_transform(dst_plane[0]); //center transform
return dst_plane[0];
}
//陷波滤波器
cv::Mat notch_kernel( cv::Mat &scr, std::vector<cv::Point> ¬ch_pot, float D0 )
{
cv::Mat notch_pass(scr.size(),CV_32FC2);
int row_num = scr.rows;
int col_num = scr.cols;
int n = 4;
for(int i=0; i<row_num; i++ ){
float *p = notch_pass.ptr<float>(i);
for(int j=0; j<col_num; j++ ){
float h_nr = 1.0;
for( unsigned int k = 0; k < notch_pot.size(); k++ ){
int u_k = notch_pot.at(k).y;
int v_k = notch_pot.at(k).x;
float dk = sqrt( pow((i-u_k),2) +pow((j-v_k),2) );
float d_k = sqrt( pow((i+u_k),2) +pow((j+v_k),2) );
float dst_dk = 1.0/(1.0 + pow(D0/dk, 2*n));
float dst_d_k = 1.0/(1.0 + pow(D0/d_k, 2*n));
h_nr = dst_dk * dst_d_k * h_nr;
// std::cout << "notch_pot: " << notch_pot.at(k) << std::endl;
}
p[2*j] = h_nr;
p[2*j+1] = h_nr;
}
}
cv::Mat temp[] = { cv::Mat::zeros(scr.size(), CV_32FC1),cv::Mat::zeros(scr.size(), CV_32FC1) };
cv::split(notch_pass, temp);
std::string name = "notch滤波器d0=" + std::to_string(D0);
cv::Mat show;
cv::normalize(temp[0], show, 1, 0, CV_MINMAX);
cv::imshow(name, show);
return notch_pass;
}
std::string name_win("Notch_filter");
cv::Rect g_rectangle;
bool g_bDrawingBox = false;//是否进行绘制
cv::RNG g_rng(12345);
std::vector<cv::Point> notch_point;
void on_MouseHandle(int event, int x, int y, int flags, void*param);
void DrawRectangle(cv::Mat& img, cv::Rect box);
void on_MouseHandle(int event, int x, int y, int falgs, void* param)
{
Mat& image = *(cv::Mat*)param;
switch (event){ //鼠标移动消息
case cv::EVENT_MOUSEMOVE:{
if (g_bDrawingBox){ //标识符为真,则记录下长和宽到Rect型变量中
g_rectangle.width = x - g_rectangle.x;
g_rectangle.height = y - g_rectangle.y;
}
}
break;
case cv::EVENT_LBUTTONDOWN:{
g_bDrawingBox = true;
g_rectangle = cv::Rect(x, y, 0, 0);//记录起点
std::cout << "start point( " << g_rectangle.x << "," << g_rectangle.y << ")" << std::endl;
}
break;
case cv::EVENT_LBUTTONUP: {
g_bDrawingBox = false;
bool w_less_0 = false, h_less_0 = false;
if (g_rectangle.width < 0){ //对宽高小于0的处理
g_rectangle.x += g_rectangle.width;
g_rectangle.width *= -1;
w_less_0 = true;
}
if (g_rectangle.height < 0){
g_rectangle.y += g_rectangle.height;
g_rectangle.height *= -1;
h_less_0 = true;
}
std::cout << "end Rect[ " << g_rectangle.x << "," << g_rectangle.y << "," << g_rectangle.width<< ","<< g_rectangle.height << "]" <<std::endl;
if( g_rectangle.height > 0 && g_rectangle.width > 0 ){
Mat imageROI = image(g_rectangle).clone();
double min, max;
cv::Point min_loc, max_loc;
cv::minMaxLoc( imageROI, &min, &max, &min_loc, &max_loc);
cv::circle(imageROI, max_loc, 10, 1);
max_loc.x += g_rectangle.x;
max_loc.y += g_rectangle.y;
notch_point.push_back(max_loc);
cv::circle(image, max_loc, 10, 1);
// cv::imshow( "ROI", imageROI );
cv::imshow( "src", image );
}
}
break;
}
}
cv::Mat notchFilter_BTW(int rows,int cols,std::vector<cv::Point> np,
float* D,int n=1,int cvtype=CV_32FC2)
{
cv::Mat filt(rows,cols,cvtype,cv::Scalar::all(0));
int cx=cols/2,cy=rows/2;
int numNotch=np.size();
float* D2=D;
for(int i=0;i<numNotch;i++)
{
D2[i]=D[i]*D[i];
}
int uk[numNotch],vk[numNotch];//在画面上的实际陷波坐标点
int u_k[numNotch],v_k[numNotch];//陷波共轭点
float Dk[numNotch],D_k[numNotch];//陷波半径r
float Hk[numNotch],H_k[numNotch];
for(int i=0;i<rows;i++)
{
for(int j=0;j<cols;j++)
{
int u=cx-j,v=cy-i;//中心坐标
for(int s=0;s<numNotch;s++)
{
uk[s]=u+np[s].x,vk[s]=v+np[s].y;
Dk[s]=uk[s]*uk[s]+vk[s]*vk[s];//距离中心半径的平方
Hk[s]=1-1/(1+pow(Dk[s]/D2[s],n));
u_k[s]=u-np[s].x,v_k[s]=v-np[s].y;
D_k[s]=u_k[s]*u_k[s]+v_k[s]*v_k[s];
H_k[s]=1-1/(1+pow(D_k[s]/D2[s],n));
}
//.data返回的是uchar*型指针,所以要强制转换成浮点数型
float* p=(float*)(filt.data+i*filt.step[0]+j*filt.step[1]);
for(int c=0;c<filt.channels();c++)
{
p[c]=Hk[0]*H_k[0];
for(int k=1;k<numNotch;k++)
{
p[c]*=Hk[k]*H_k[k];
}
}
}
}
return filt;
}
int main(int argc, char * argv[])
{
// if( argc != 2){
// std::cerr << "Usage: " << argv[0] << " <img_name> " << std::endl;
// return -1;
// }
//Mat srcImage = cv::imread(argv[1], cv::IMREAD_GRAYSCALE);
Mat srcImage = cv::imread("../../im/30.bmp", 0);
if( srcImage.empty() )
return -1;
imshow( "src_img", srcImage );
Mat padded = image_add_border(srcImage);
center_transform( padded );
Mat plane[]={padded, cv::Mat::zeros(padded.size(), CV_32FC1)};
Mat complexIm;
merge(plane,2,complexIm);
cv::dft(complexIm,complexIm);//fourior transform
Mat temp[] = {cv::Mat::zeros(complexIm.size(),CV_32FC1),cv::Mat::zeros(complexIm.size(),CV_32FC1)};
//显示频谱图
split(complexIm, temp);
Mat aa;
magnitude(temp[0], temp[1], aa);
divide(aa, aa.cols*aa.rows, aa);
```
**4.空间滤波**
```
//中值滤波
Mat dstImage_median;
medianBlur(srcImage, dstImage_median, 3);
imshow("median", dstImage_median);
//waitKey(0);
//高斯滤波
Mat dstImage_Gaussian;
GaussianBlur(dstImage_median, dstImage_Gaussian, Size(3,3),0 , 0);//高斯滤波多一点参数要写
imshow("Gauss", dstImage_Gaussian);
//waitKey(0);
//均值滤波
Mat dstImage_mean;
blur( dstImage_median, dstImage_mean, Size(3,3));
imshow("mean", dstImage_mean);
//waitKey(0);
```
**5.腐蚀膨胀**
```
//腐蚀膨胀前要二值化
Mat dstImage_binary;
threshold(dstImage_mean, dstImage_binary, 60, 255, CV_THRESH_BINARY);
imshow("binary", dstImage_binary);
//waitKey(0);
//腐蚀
Mat element = getStructuringElement(MORPH_RECT, Size(3,3));
Mat dstImage_erode;
erode (dstImage_binary, dstImage_erode, element);
imshow("erode",dstImage_erode);
//waitKey();
//膨胀
Mat dstImage_dilate;
dilate (dstImage_erode, dstImage_dilate, element);
imshow("dilate",dstImage_dilate);
//waitKey();
```
5.去除小联通域
```
void RemoveSmallRegion(Mat &Src, Mat &Dst, int AreaLimit, int CheckMode, int NeihborMode)
{
int RemoveCount = 0;
//新建一幅标签图像初始化为0像素点,为了记录每个像素点检验状态的标签,0代表未检查,1代表正在检查,2代表检查不合格(需要反转颜色),3代表检查合格或不需检查
//初始化的图像全部为0,未检查
Mat PointLabel = Mat::zeros(Src.size(), CV_8UC1);
if (CheckMode == 1)//去除小连通区域的白色点
{
//cout << "去除小连通域.";
for (int i = 0; i < Src.rows; i++)
{
for (int j = 0; j < Src.cols; j++)
{
if (Src.at<uchar>(i, j) < 10)
{
PointLabel.at<uchar>(i, j) = 3;//将背景黑色点标记为合格,像素为3
}
}
}
}
else//去除孔洞,黑色点像素
{
//cout << "去除孔洞";
for (int i = 0; i < Src.rows; i++)
{
for (int j = 0; j < Src.cols; j++)
{
if (Src.at<uchar>(i, j) > 10)
{
PointLabel.at<uchar>(i, j) = 3;//如果原图是白色区域,标记为合格,像素为3
}
}
}
}
vector<Point2i>NeihborPos;//将邻域压进容器
NeihborPos.push_back(Point2i(-1, 0));
NeihborPos.push_back(Point2i(1, 0));
NeihborPos.push_back(Point2i(0, -1));
NeihborPos.push_back(Point2i(0, 1));
if (NeihborMode == 1)
{
//cout << "Neighbor mode: 8邻域." << endl;
NeihborPos.push_back(Point2i(-1, -1));
NeihborPos.push_back(Point2i(-1, 1));
NeihborPos.push_back(Point2i(1, -1));
NeihborPos.push_back(Point2i(1, 1));
}
else int a = 0;//cout << "Neighbor mode: 4邻域." << endl;
int NeihborCount = 4 + 4 * NeihborMode;
int CurrX = 0, CurrY = 0;
//开始检测
for (int i = 0; i < Src.rows; i++)
{
for (int j = 0; j < Src.cols; j++)
{
if (PointLabel.at<uchar>(i, j) == 0)//标签图像像素点为0,表示还未检查的不合格点
{ //开始检查
vector<Point2i>GrowBuffer;//记录检查像素点的个数
GrowBuffer.push_back(Point2i(j, i));
PointLabel.at<uchar>(i, j) = 1;//标记为正在检查
int CheckResult = 0;
for (int z = 0; z < GrowBuffer.size(); z++)
{
for (int q = 0; q < NeihborCount; q++)
{
CurrX = GrowBuffer.at(z).x + NeihborPos.at(q).x;
CurrY = GrowBuffer.at(z).y + NeihborPos.at(q).y;
if (CurrX >= 0 && CurrX<Src.cols&&CurrY >= 0 && CurrY<Src.rows) //防止越界
{
if (PointLabel.at<uchar>(CurrY, CurrX) == 0)
{
GrowBuffer.push_back(Point2i(CurrX, CurrY)); //邻域点加入buffer
PointLabel.at<uchar>(CurrY, CurrX) = 1; //更新邻域点的检查标签,避免重复检查
}
}
}
}
if (GrowBuffer.size()>AreaLimit) //判断结果(是否超出限定的大小),1为未超出,2为超出
CheckResult = 2;
else
{
CheckResult = 1;
RemoveCount++;//记录有多少区域被去除
}
for (int z = 0; z < GrowBuffer.size(); z++)
{
CurrX = GrowBuffer.at(z).x;
CurrY = GrowBuffer.at(z).y;
PointLabel.at<uchar>(CurrY, CurrX) += CheckResult;//标记不合格的像素点,像素值为2
}
//********结束该点处的检查**********
}
}
}
CheckMode = 255 * (1 - CheckMode);
//开始反转面积过小的区域
for (int i = 0; i < Src.rows; ++i)
{
for (int j = 0; j < Src.cols; ++j)
{
if (PointLabel.at<uchar>(i, j) == 2)
{
Dst.at<uchar>(i, j) = CheckMode;
}
else if (PointLabel.at<uchar>(i, j) == 3)
{
Dst.at<uchar>(i, j) = Src.at<uchar>(i, j);
}
}
}
//cout << RemoveCount << " objects removed." << endl;
}
//===============去除小连通域===============================
```
---
**碎碎念:**
好啦,差不多就这些啦,完结撒花~
其实还有蛮多地方没弄懂的,有点函数也是直接找过来的,并不清楚,等我过两天(真的是过两天吗???)搞懂了再来更新,我一定会回来的!