為什麼要進行高頻提升和高頻加強?
高頻濾波後的圖像,其背景平均強度 減小到接近黑色(因為高通濾波器濾除 了傅裡葉變換的零頻率成分: F(0,0)=f(x,y)=0)
解決辦法:把原始圖像加到過濾後的 結果,如拉普拉斯算子增強,這種處理 稱為高頻提升過濾
鈍化模闆(銳化或高通圖像):從一幅圖像減去 其自身模糊圖像而生成的銳化圖像構成。在頻率 域,即從圖像本身減去低通濾波(模糊)後的圖 像而得到高通濾波(銳化)的圖像。
鈍化模闆:
與
式中
是一個低通濾波器,F(u,v)是f(x,y)的傅裡葉變換。這裡
是平滑後的圖像,該圖像類似于
高頻提升過濾:
增權重重系數A
是以,當A=1,即高通過濾;當A>1,累加圖像本身
由
高頻提升過濾定義為:
系數k推導:
在原圖像上加上該模闆的權重部分
當k=1時 原圖像加上鈍化模闆,k>1時高提升濾波器;
在頻率域有
稱為高頻強調濾波器。如前面說的那樣,高頻濾波器将直流項設定為0,這樣就把濾波後圖像的平均灰階減小為0,。高頻強調濾波器不存在這一問題,因為高通項上加了1。常數k決定了影響最終結果的高頻比例。
高頻提升加強:
更為一般的公式
a >= 0, b >a
用圖像的高頻成分進行增強
當a=A-1,b=1時轉化為高頻提升過濾
當b>1,高頻得到加強,b控制高頻的貢獻
增加a的目的是使零頻率不被濾波器過濾,a控制矩原點的偏移量
例子:
使用高斯高通,D0=40,a=0.5,b=0.75
效果圖:
代碼實作:
#include "opencv2/opencv.hpp"
cv::Mat image_add_border( cv::Mat &src )
{
int w=2*src.cols;
int h=2*src.rows;
std::cout << "src: " << src.cols << "*" << src.rows << std::endl;
cv::Mat padded;
copyMakeBorder( src, padded, 0, h-src.rows, 0, w-src.cols,
cv::BORDER_CONSTANT, cv::Scalar::all(0));
padded.convertTo(padded,CV_32FC1);
std::cout << "opt: " << 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);
for(int j=0; j<src.cols; j++){
p[j] = p[j] * pow(-1, i+j);
}
}
}
//對角線交換内容
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);
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 gaussian_high_kernel( cv::Mat &scr, float D0 )
{
cv::Mat gaussian_high_pass(scr.size(),CV_32FC2);
int row_num = scr.rows;
int col_num = scr.cols;
float d0 = 2 * D0 * D0;
for(int i=0; i<row_num; i++ ){
float *p = gaussian_high_pass.ptr<float>(i);
for(int j=0; j<col_num; j++ ){
float d = pow((i - row_num/2),2) + pow((j - col_num/2),2);
p[2*j] = 0.5 + 0.75*(1 - expf(-d/d0));
p[2*j+1] = 0.5 + 0.75*(1 - expf(-d/d0));
}
}
cv::Mat temp[] = { cv::Mat::zeros(scr.size(), CV_32FC1),
cv::Mat::zeros(scr.size(), CV_32FC1) };
cv::split(gaussian_high_pass, temp);
std::string name = "高斯高通濾波器d0=" + std::to_string(D0);
cv::Mat show;
cv::normalize(temp[0], show, 1, 0, CV_MINMAX);
cv::imshow(name, show);
return gaussian_high_pass;
}
//高斯高通濾波器
cv::Mat gaussian_highpass_filter(cv::Mat &src, float D0 )
{
cv::Mat padded = image_add_border(src);
center_transform( padded );
cv::Mat gaussian_kernel = gaussian_high_kernel(padded,D0 );
cv::Mat result = frequency_filter(padded,gaussian_kernel);
return result;
}
int main(int argc, char * argv[])
{
if( argc != 3 ){
std::cerr << "Usage: " << argv[0] << " <img_name> <radius> " << std::endl;
return -1;
}
cv::Mat image = cv::imread(argv[1], cv::IMREAD_GRAYSCALE);
if( image.empty() )
return -1;
float radius = std::stof( argv[2] );
// cv::resize( image, image, cv::Size(), 0.25, 0.25);
cv::imshow("src",image);
cv::Mat src_eq;
cv::equalizeHist( image, src_eq );
cv::imshow("src_eq",src_eq);
cv::Mat gaussian_result = gaussian_highpass_filter(image, radius);
gaussian_result = gaussian_result( cv::Rect(0, 0, image.cols, image.rows) );
cv::normalize(gaussian_result, gaussian_result, 255, 0, CV_MINMAX);
gaussian_result.convertTo(gaussian_result, CV_8U);
cv::equalizeHist( gaussian_result, gaussian_result );
cv::imshow("gaussian_result",gaussian_result);
cv::waitKey(0);
return 1;
}