天天看點

OpenCV 頻率域實作鈍化模闆、高提升濾波和高頻強調濾波 C++

為什麼要進行高頻提升和高頻加強?

高頻濾波後的圖像,其背景平均強度 減小到接近黑色(因為高通濾波器濾除 了傅裡葉變換的零頻率成分: F(0,0)=f(x,y)=0)

解決辦法:把原始圖像加到過濾後的 結果,如拉普拉斯算子增強,這種處理 稱為高頻提升過濾

鈍化模闆(銳化或高通圖像):從一幅圖像減去 其自身模糊圖像而生成的銳化圖像構成。在頻率 域,即從圖像本身減去低通濾波(模糊)後的圖 像而得到高通濾波(銳化)的圖像。

OpenCV 頻率域實作鈍化模闆、高提升濾波和高頻強調濾波 C++

鈍化模闆:

OpenCV 頻率域實作鈍化模闆、高提升濾波和高頻強調濾波 C++

與                                            

OpenCV 頻率域實作鈍化模闆、高提升濾波和高頻強調濾波 C++

式中

OpenCV 頻率域實作鈍化模闆、高提升濾波和高頻強調濾波 C++

 是一個低通濾波器,F(u,v)是f(x,y)的傅裡葉變換。這裡

OpenCV 頻率域實作鈍化模闆、高提升濾波和高頻強調濾波 C++

是平滑後的圖像,該圖像類似于

OpenCV 頻率域實作鈍化模闆、高提升濾波和高頻強調濾波 C++

高頻提升過濾:

增權重重系數A

OpenCV 頻率域實作鈍化模闆、高提升濾波和高頻強調濾波 C++

是以,當A=1,即高通過濾;當A>1,累加圖像本身

由                                           

OpenCV 頻率域實作鈍化模闆、高提升濾波和高頻強調濾波 C++

高頻提升過濾定義為:          

OpenCV 頻率域實作鈍化模闆、高提升濾波和高頻強調濾波 C++

系數k推導:

在原圖像上加上該模闆的權重部分

OpenCV 頻率域實作鈍化模闆、高提升濾波和高頻強調濾波 C++

當k=1時 原圖像加上鈍化模闆,k>1時高提升濾波器;

在頻率域有

OpenCV 頻率域實作鈍化模闆、高提升濾波和高頻強調濾波 C++
OpenCV 頻率域實作鈍化模闆、高提升濾波和高頻強調濾波 C++

稱為高頻強調濾波器。如前面說的那樣,高頻濾波器将直流項設定為0,這樣就把濾波後圖像的平均灰階減小為0,。高頻強調濾波器不存在這一問題,因為高通項上加了1。常數k決定了影響最終結果的高頻比例。

高頻提升加強:

更為一般的公式

OpenCV 頻率域實作鈍化模闆、高提升濾波和高頻強調濾波 C++

                    a >= 0, b >a

OpenCV 頻率域實作鈍化模闆、高提升濾波和高頻強調濾波 C++

 用圖像的高頻成分進行增強

當a=A-1,b=1時轉化為高頻提升過濾

當b>1,高頻得到加強,b控制高頻的貢獻

增加a的目的是使零頻率不被濾波器過濾,a控制矩原點的偏移量

OpenCV 頻率域實作鈍化模闆、高提升濾波和高頻強調濾波 C++

例子:

使用高斯高通,D0=40,a=0.5,b=0.75

效果圖:

OpenCV 頻率域實作鈍化模闆、高提升濾波和高頻強調濾波 C++

代碼實作:

#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;

}

繼續閱讀