频率域的拉普拉斯算子

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

//laplace滤波器
cv::Mat laplace_kernel( cv::Mat &scr )
{
    cv::Mat laplace_pass(scr.size(),CV_32FC2);
    int row_num = scr.rows;
    int col_num = scr.cols;
    for(int i=0; i<row_num; i++ ){
        float *p = laplace_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]   = -4 * CV_PI * CV_PI * d;
            p[2*j+1] = -4 * CV_PI * CV_PI * d;
        }
    }

    cv::Mat temp[] = { cv::Mat::zeros(scr.size(), CV_32FC1),
                       cv::Mat::zeros(scr.size(), CV_32FC1) };
    cv::split(laplace_pass, temp);
    std::string name = "laplace滤波器";
    cv::Mat show;
    cv::normalize(temp[0], show, 1, 0, CV_MINMAX);
    cv::imshow(name, show);
    return laplace_pass;
}

//laplace滤波器
cv::Mat laplace_filter(cv::Mat &src  )
{
    cv::Mat padded = image_add_border(src);
    center_transform( padded );
    cv::Mat lap_kernel = laplace_kernel(padded );
    cv::Mat result = frequency_filter(padded,lap_kernel);
    return result;
}

int main(int argc, char * argv[])
{
    if( argc != 1 ){
        std::cerr << "Usage: " << argv[0] << " <img_name>  " << std::endl;
        return -1;
    }

    cv::Mat image = cv::imread("458.tif", cv::IMREAD_GRAYSCALE);
    if( image.empty() )
        return -1;
    cv::resize( image, image, cv::Size(), 0.5, 0.5);
    cv::imshow("src",image);

    cv::Mat laplace_result = laplace_filter(image);
    laplace_result = laplace_result( cv::Rect(0, 0, image.cols, image.rows) );

    image.convertTo(image,CV_32FC1);
    cv::normalize(image, image, 1, 0, CV_MINMAX);
    cv::normalize(laplace_result, laplace_result, 1, 0, CV_MINMAX);
        cv::imshow("_result",laplace_result);
    cv::subtract(image,laplace_result, laplace_result );
    cv::normalize(laplace_result, laplace_result, 1, 0, CV_MINMAX);
    cv::imshow("laplace_result",laplace_result);
    cv::waitKey(0);

    return 1;
}