《數字圖像處理》岡薩雷斯 第三版中,P179有一個指紋增強的例子;
先用opencv模拟實作:
處理這個指紋圖像 ,其污染是很明顯的。自動指紋識别的關鍵一步是增強指紋的脊線并減小污染。使用高通濾波實作;
脊線增強使用脊線包含高頻這樣一個事實來完成,而高通濾波不會改變高頻;另一方面,濾掉了低頻,而低頻分量對應圖像中慢變化的灰階,例如背景和污染。
這樣增強就可以同時減少高頻外的所有特征來達到,目前情況,高頻是感興趣的特征。
原圖像---------------------------------------------------高通濾波後圖像------------------------------------------門檻值化後圖像
先用布特沃斯高通濾波器D0=50,n=4濾波,結果失去了灰階色調,因為直流項被減小為0;最終結果是高通濾波後圖像的典型暗色調已成為主流,這樣就可以通過附加處理增強感興趣細節。
一種簡單的方法是對結果進行門檻值操作,将所有負值設為黑,所有正值設為白。
代碼實作:
#include "opencv2/opencv.hpp"
using namespace cv;
//transform to center 中心化
void center_transform( cv::Mat &src )
{
for(int i=0; i<src.rows; i++){
float *p = src.ptr<float>(i);
for(int j=0; j<src.cols; j++){
p[j] = p[j] * pow(-1, i+j);
}
}
}
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);
}
int main(int argc, char * argv[])
{
const char* filename = argc >=2 ? argv[1] : "zhiwen12.tif";
Mat image = imread(filename, IMREAD_GRAYSCALE);
if( image.empty())
return -1;
resize( image, image, Size(), 0.25, 0.25);
imshow("src",image);
image.convertTo(image,CV_32FC1);
///快速傅裡葉變換/
// int oph = getOptimalDFTSize(image.rows);
// int opw = getOptimalDFTSize(image.cols);
int oph = 2*image.rows;
int opw = 2*image.cols;
Mat padded;
copyMakeBorder(image, padded, 0, oph-image.rows, 0, opw-image.cols,
BORDER_CONSTANT, Scalar::all(0));
// center_transform(padded);
Mat pad_show;
normalize(padded, pad_show, 1, 0, CV_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);
divide(aa, oph*opw, aa);
imshow("pu",aa);
/頻域濾波///
//生成頻域濾波核
Mat butter_sharpen(padded.size(), CV_32FC2);
float D0 = 50.;
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++){
float 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 tmp[2];
// split(butter_sharpen, tmp);
// multiply(temp[0], tmp[0], temp[0]);
// multiply(temp[1], tmp[0], temp[1]);
merge(temp, 2, complexI);
multiply(complexI, butter_sharpen, complexI);
//傅裡葉反變換
idft(complexI, complexI, CV_DXT_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, 1, 0, CV_MINMAX);
show = show(Rect(0,0,image.cols, image.rows));
imshow("_filter",show);
threshold(dstSharpen[0], dstSharpen[0], 0, 255, THRESH_BINARY);
normalize(dstSharpen[0], dstSharpen[0], 0, 1, CV_MINMAX);
dstSharpen[0] = dstSharpen[0](Rect(0,0,image.cols, image.rows));
imshow("dstSharpen",dstSharpen[0]);
// imshow( "dest_finger", dest_finger );
waitKey(0);
return 1;
}