Opencv + FFTW3 图象高斯高低通滤波

1：LPF

公式：out(i,j)= exp(-((i-M/2)^2+(j-N/2)^2)/2/sigma^2);

代码：

        

void LPF_Fliter(vector<double>&data, int cols, int rows, double gamma){float gamma22 = 2 * gamma*gamma;float temp = 0.0;int halfRows = cvRound(rows / 2);int halfCols = cvRound(cols / 2);char debug[12] = { 0 };for (int j = 1; j <= rows; j++){for (int i = 1; i <= cols; i++){temp = (pow(j - halfRows, 2) + pow(i - halfCols, 2));data.push_back(exp(-temp / gamma22));}}}

2：HPF:

公式：out(i,j)= 1- exp(-((i-M/2)^2+(j-N/2)^2)/2/sigma^2);

代码：

 

void HPF_Fliter(vector<double>&data, int cols, int rows, double gamma){#if 1float gamma22 = 2 * gamma*gamma;float temp = 0.0;int halfRows = cvRound(rows / 2);int halfCols = cvRound(cols / 2);char debug[12] = { 0 };for (int j = 1; j <= rows; j++){for (int i = 1; i <= cols; i++){temp = (pow(j - halfRows, 2) + pow(i - halfCols, 2));data.push_back(1.0 - exp(-temp / gamma22));}}#elsefstream readFile("C:\\Users\\Tony\\Desktop\\debug.txt", ios::in);string temp;int lineCount = 0;while (getline(readFile, temp)){vector<string>dest;split(temp, "\t", dest);lineCount++;if (dest.size() > 1){for (int i = 0; i < dest.size(); i++){//printf_s("%s\tn", dest.at(i));data.push_back(atof(dest.at(i).c_str()));}}//printf_s("\n");}#endif}

//生成频谱图：

                   

//生成频谱模板void test(){int  i;fftw_complex*din, *out;fftw_plan p,backp;Mat image = imread("C:\\Users\\Tony\\Desktop\\Image\\lenaCopy.bmp", CV_LOAD_IMAGE_GRAYSCALE);din = (fftw_complex *)fftw_malloc(sizeof(fftw_complex)*image.cols*image.rows);out = (fftw_complex *)fftw_malloc(sizeof(fftw_complex)*image.cols*image.rows);Mat floatImage(image.size(), CV_32FC1);image.convertTo(floatImage, CV_32FC1);for (size_t j = 0; j < image.rows; j++){for (size_t i = 0; i < image.cols; i++){din[i + j*image.cols][0] = *floatImage.ptr<float>(j, i);din[i + j*image.cols][1] = 0;}}//forward fftp = fftw_plan_dft_2d(image.rows, image.cols, din, out, FFTW_FORWARD, FFTW_ESTIMATE);fftw_execute(p);Mat Resource(image.size(), CV_32FC1);Mat Imsource(image.size(), CV_32FC1);for (size_t j = 0; j < image.rows; j++){for (size_t i = 0; i < image.cols; i++){*Resource.ptr<float>(j, i) = out[i + j*image.cols][0];//实部*Imsource.ptr<float>(j, i) = out[i + j*image.cols][1];//虚部}}Mat Redest(image.size(), CV_32FC1);Mat Imdest(image.size(), CV_32FC1);fftshift(Resource, Redest);fftshift(Imsource, Imdest);Mat absImage(image.size(), CV_32FC1);for (size_t j = 0; j < image.rows; j++){for (size_t i = 0; i < image.cols; i++){*absImage.ptr<float>(j, i) = sqrt(pow(*Redest.ptr<float>(j, i), 2) + pow(*Imdest.ptr<float>(j, i), 2));}}absImage += Scalar::all(1);Mat logImage(image.size(), CV_32FC1);log(absImage, logImage);//SaveData(logImage);double min = 0.0, max = 0.0;minMaxLoc(logImage, &min, &max);double scale = 255/(max-min);double shift = -min*scale;Mat myImage(image.size(), CV_8UC1);convertScaleAbs(logImage, myImage, scale, shift);imwrite("C:\\Users\\Tony\\Desktop\\lenaCopy.bmp", myImage);}

高低通滤波：

void testOne(){int  i;fftw_complex*din, *out;fftw_plan p, backp;Mat image = imread("C:\\Users\\Tony\\Desktop\\blur.bmp", CV_LOAD_IMAGE_GRAYSCALE);din = (fftw_complex *)fftw_malloc(sizeof(fftw_complex)*image.cols*image.rows);out = (fftw_complex *)fftw_malloc(sizeof(fftw_complex)*image.cols*image.rows);Mat floatImage(image.size(), CV_32FC1);image.convertTo(floatImage, CV_32FC1);for (size_t j = 0; j < image.rows; j++){for (size_t i = 0; i < image.cols; i++){din[i + j*image.cols][0] = *floatImage.ptr<float>(j, i);din[i + j*image.cols][1] = 0.0;}}//forward fftp = fftw_plan_dft_2d(image.rows, image.cols, din, out, FFTW_FORWARD, FFTW_ESTIMATE);fftw_execute(p);//图像频谱中心化Mat Resource(image.size(), CV_32FC1);Mat Imsource(image.size(), CV_32FC1);//虚部Mat Redest(image.size(), CV_32FC1);Mat Imdest(image.size(), CV_32FC1);//虚部for (size_t j = 0; j < image.rows; j++){for (size_t i = 0; i < image.cols; i++){*Resource.ptr<float>(j, i) = out[i + j*image.cols][0];//实部*Imsource.ptr<float>(j, i) = out[i + j*image.cols][1];//虚部}}fftshift(Resource, Redest);fftshift(Imsource, Imdest);//虚部//图像频谱乘以不通的数据,构建滤波vector<double>data;HPF_Fliter(data, image.cols, image.rows, 10);//LPF_Fliter(data, image.cols, image.rows, 100);Mat filterRe(image.size(), CV_32FC1);Mat filterIm(image.size(), CV_32FC1);//虚部for (size_t j = 0; j < image.rows; j++){for (size_t i = 0; i < image.cols; i++){*filterRe.ptr<float>(j, i) = *Redest.ptr<float>(j, i) * data.at(j*image.cols + i);//实部*filterIm.ptr<float>(j, i) = *Imdest.ptr<float>(j, i) * data.at(j*image.cols + i);//虚部 }}//图像频谱逆中心化Mat NotCenterRedest(image.size(), CV_32FC1);Mat NotCenterImdest(image.size(), CV_32FC1);//虚部ifftShift(filterRe, NotCenterRedest);ifftShift(filterIm, NotCenterImdest);//虚部//重新构造复数，傅里叶反变换fftw_complex *spectrum = (fftw_complex *)fftw_malloc(sizeof(fftw_complex)*image.cols*image.rows);for (size_t j = 0; j < image.rows; j++){for (size_t i = 0; i < image.cols; i++){spectrum[j*image.cols + i][0] = *NotCenterRedest.ptr<float>(j, i);//实部spectrum[j*image.cols + i][1] = *NotCenterImdest.ptr<float>(j, i);//虚部}}//傅里叶反变换fftw_complex *result = (fftw_complex *)fftw_malloc(sizeof(fftw_complex)*image.cols*image.rows);backp = fftw_plan_dft_2d(image.rows, image.cols, spectrum, result, FFTW_BACKWARD, FFTW_ESTIMATE);fftw_execute(backp);int size = image.cols*image.rows;Mat IvsResource(image.size(), CV_32FC1, Scalar::all(0));//Mat IvsImsource(image.size(), CV_32FC1, Scalar::all(0));//虚部//Mat ifftImage(image.size(), CV_32FC1, Scalar::all(0));Mat myImage(image.size(), CV_8UC1);for (size_t j = 0; j < image.rows; j++){for (size_t i = 0; i < image.cols; i++){*IvsResource.ptr<float>(j, i) = result[i + j*image.cols][0] / size;//实部//*IvsImsource.ptr<float>(j, i) = result[i + j*image.cols][1] / size;//虚部//*ifftImage.ptr<float>(j, i) = sqrt(pow(result[i + j*image.cols][0], 2) + pow(result[i + j*image.cols][1], 2));*myImage.ptr<byte>(j, i) = cvRound(abs(*IvsResource.ptr<float>(j, i)));//实部}}//SaveData(IvsResource);imwrite("C:\\Users\\Tony\\Desktop\\Copy.bmp", myImage);Mat gammaImage(image.size(), CV_8UC1);GammaCorrection(myImage, gammaImage, 0.3);//释放资源fftw_destroy_plan(p);fftw_destroy_plan(backp);fftw_free(din);fftw_free(out);fftw_free(result);}