#include "opencv2/core/core.hpp"
#include "opencv2/imgproc/imgproc.hpp"
#include "opencv2/highgui/highgui.hpp"
#include <iostream>
#include "highgui.h"
#include "cv.h"
using namespace cv;
using namespace std;
IplImage* img = cvLoadImage("C:\\Users\\CDZHYF\\Desktop\\101.png",0);
double HEIGHT = img->height;
double WIDTH = img->width;
void Calculate_grad_square(double *pData, double *pGrad_data, short isX)
{
int i, j;
IplImage *src = cvLoadImage("C:\\Users\\CDZHYF\\Desktop\\101.png",0);
int HEIGHT_IN = src ->height;
int WIDTH_IN = src -> width;
if (isX == 1)
{
for (i = 0; i < HEIGHT_IN; i++)
{
for (j = 0; j < WIDTH_IN; j++)
{
if (i == HEIGHT_IN - 1)
{
pGrad_data[i * WIDTH_IN + j] = 0.0;
}
else
{
pGrad_data[i * WIDTH_IN + j] = pow((pData[(i + 1) * WIDTH_IN + j] - pData[i * WIDTH_IN + j]), 2);
}
}
}
}
else
{
for (i = 0; i < HEIGHT_IN; i++)
{
for (j = 0; j < WIDTH_IN; j++)
{
if (i == WIDTH_IN - 1)
{
pGrad_data[i * WIDTH_IN + j] = 0.0;
}
else
{
pGrad_data[i * WIDTH_IN + j] = pow((pData[i * WIDTH_IN + j + 1] - pData[i * WIDTH_IN + j]), 2);
}
}
}
}
}
void PrintMat(cv::Mat m)
{
int w = m.cols;
int h = m.rows;
for (int y=0; y<h; y++)
{
for (int x=0; x<w; x++)
{
if (m.type() == CV_64FC1)
{
printf("%f ", m.at<double>(y,x));
}
else if (m.type() == CV_8UC1)
{
printf("%d ", m.at<uchar>(y,x));
}
else if (m.type() == CV_32FC1)
{
printf("%f ", m.at<float>(y,x));
}
}
printf("\n");
}
}
void FftTest()
{
Mat I = imread("C:\\Users\\CDZHYF\\Desktop\\101.png", CV_LOAD_IMAGE_GRAYSCALE);
if( I.empty())
{
printf("img load error \n");
system("pause");
return;
}
Mat padded; //expand input image to optimal size
int m = getOptimalDFTSize( I.rows );
int n = getOptimalDFTSize( I.cols ); // on the border add zero values
copyMakeBorder(I, padded, 0, m - I.rows, 0, n - I.cols, BORDER_CONSTANT, Scalar::all(0));
Mat planes[] = {Mat_<float>(padded), Mat::zeros(padded.size(), CV_32F)};
Mat complexI;
merge(planes, 2, complexI); // Add to the expanded another plane with zeros
dft(complexI, complexI); // this way the result may fit in the source matrix
// compute the magnitude and switch to logarithmic scale
// => log(1 + sqrt(Re(DFT(I))^2 + Im(DFT(I))^2))
split(complexI, planes); // planes[0] = Re(DFT(I), planes[1] = Im(DFT(I))
magnitude(planes[0], planes[1], planes[0]);// planes[0] = magnitude
Mat magI = planes[0];
magI += Scalar::all(1); // switch to logarithmic scale
log(magI, magI);
// crop the spectrum, if it has an odd number of rows or columns
magI = magI(Rect(0, 0, magI.cols & -2, magI.rows & -2));
// rearrange the quadrants of Fourier image so that the origin is at the image center
int cx = magI.cols/2;
int cy = magI.rows/2;
Mat q0(magI, Rect(0, 0, cx, cy)); // Top-Left - Create a ROI per quadrant
Mat q1(magI, Rect(cx, 0, cx, cy)); // Top-Right
Mat q2(magI, Rect(0, cy, cx, cy)); // Bottom-Left
Mat q3(magI, Rect(cx, cy, cx, cy)); // Bottom-Right
Mat tmp; // swap quadrants (Top-Left with Bottom-Right)
q0.copyTo(tmp);
q3.copyTo(q0);
tmp.copyTo(q3);
q1.copyTo(tmp); // swap quadrant (Top-Right with Bottom-Left)
q2.copyTo(q1);
tmp.copyTo(q2);
normalize(magI, magI, 0, 1, CV_MINMAX); // Transform the matrix with float values into a
// viewable image form (float between values 0 and 1).
imshow("Input Image" , I ); // Show the result
imshow("spectrum magnitude", magI);
waitKey(0);
}
//傅里叶正变换
void Fft2(IplImage *src, IplImage *dst)
{
IplImage *image_Re = 0, *image_Im = 0, *Fourier = 0;
image_Re = cvCreateImage(cvGetSize(src), IPL_DEPTH_64F, 1); //实部
image_Im = cvCreateImage(cvGetSize(src), IPL_DEPTH_64F, 1); //虚部
IplImage* im = cvLoadImage("C:\\Users\\CDZHYF\\Desktop\\101.png",0);
int HEIGHT_IN = im ->height;
int WIDTH_IN = im -> width;
int i, j, k;
double *pDouble_in = (double*) malloc(WIDTH_IN * HEIGHT_IN * sizeof(double));
double *pFx = (double*) malloc(WIDTH_IN * HEIGHT_IN * sizeof(double));
double *pFy = (double*) malloc(WIDTH_IN * HEIGHT_IN * sizeof(double));
uchar* pData = (uchar*)src->imageData;
for (i = 0; i < HEIGHT_IN; i++)
{
for (j = 0; j< WIDTH_IN; j++)
{
pDouble_in[i * WIDTH_IN + j] = pData[i * src->widthStep + j];
}
}
//printf("YES\n");
Calculate_grad_square(pDouble_in, pFx, 1);
Calculate_grad_square(pDouble_in, pFy, 0);
//2 channels (image_Re, image_Im)
Fourier = cvCreateImage(cvGetSize(src), IPL_DEPTH_64F, 2);
// Real part conversion from u8 to 64f (double)
cvConvertScale(src, image_Re, 1, 0);
// Imaginary part (zeros)
cvZero(image_Im);
// Join real and imaginary parts and stock them in Fourier image
cvMerge(image_Re, image_Im, 0, 0, Fourier);
// Application of the forward Fourier transform
cvDFT(Fourier, dst, CV_DXT_FORWARD);
//for (k = 0; k < 1 ;k++)
//{
// for (i = 0; i < HEIGHT_IN; i++)
// {
// for (j = 0; j< WIDTH_IN; j++)
// {
// if ((pFx[i * WIDTH_IN + j] + pFy[i * WIDTH_IN + j]) == 0.0)
// {
// //printf("%f , %f\n",pFx[i * WIDTH_IN + j] , pFy[i * WIDTH_IN + j]);
// //dst->imageData[i * dst->widthStep + j] = 0.0;
// //pUf_imag[i * WIDTH_IN + j] = 0.0;
// }
// else
// {
// dst->imageData[i * dst->widthStep + j] = pFx[i * WIDTH_IN + j] / (pFx[i * WIDTH_IN + j] + pFy[i * WIDTH_IN + j]) * (dst->imageData[i * dst->widthStep + j]);
// //pUf_imag[i * WIDTH_IN + j] = (pFx[i * WIDTH_IN + j] / (pFx[i * WIDTH_IN + j] + pFy[i * WIDTH_IN + j])) * pIf_imag[i * WIDTH_IN + j];
// //printf("(%d,%d)= %f +i %f\n", i, j, pUf_real[i * WIDTH_IN + j], pUf_imag[i * WIDTH_IN + j]);
// //printf("%f\n", pow((pFx[i * WIDTH_IN + j] / (pFx[i * WIDTH_IN + j] + pFy[i * WIDTH_IN + j])),1));
// }
// }
// }
//}
//cvEqualizeHist(dst,dst);
free(pDouble_in);
free(pFx);
free(pFy);
cvReleaseImage(&image_Re);
cvReleaseImage(&image_Im);
cvReleaseImage(&Fourier);
}
void CalcFftFrequencySpectrum(IplImage *src, IplImage *dst, IplImage* pTmp)
{
IplImage *image_Re = 0, *image_Im = 0;
int nRow, nCol, i, j, cy, cx;
double scale, shift, tmp13, tmp24;
double sum = 0;
double hei, wid;
double sumwinr = 0.0;
double sumwini = 0.0;
double* meanx = (double*)malloc(HEIGHT * sizeof(double));
double *meany = (double*)malloc(WIDTH * sizeof(double));
memset(meanx, 0.0, HEIGHT * sizeof(double));
memset(meany, 0.0, WIDTH * sizeof(double));
int wi, wj;
printf("YES\n");
//IplImage* pTmp = 0;
//pTmp = cvCreateImage(cvGetSize(src), IPL_DEPTH_64F, 2);
image_Re = cvCreateImage(cvGetSize(src), IPL_DEPTH_64F, 1);
//Imaginary part
image_Im = cvCreateImage(cvGetSize(src), IPL_DEPTH_64F, 1);
cvSplit( src, image_Re, image_Im, 0, 0 );
//具体原理见冈萨雷斯数字图像处理p123
// Compute the magnitude of the spectrum Mag = sqrt(Re^2 + Im^2)
//计算傅里叶谱
cvPow( image_Re, image_Re, 2.0);
cvPow( image_Im, image_Im, 2.0);
cvAdd( image_Re, image_Im, image_Re);
cvPow( image_Re, image_Re, 0.5 );
//对数变换以增强灰度级细节(这种变换使以窄带低灰度输入图像值映射
//一宽带输出值,具体可见冈萨雷斯数字图像处理p62)
// Compute log(1 + Mag);
double minVal = 0, maxVal = 0;
cvAddS( image_Re, cvScalar(1.0), image_Re ); // 1 + Mag
cvLog( image_Re, image_Re ); // log(1 + Mag)
//Rearrange the quadrants of Fourier image so that the origin is at the image center
nRow = src->height;
nCol = src->width;
cy = nRow/2; // image center
cx = nCol/2;
//CV_IMAGE_ELEM为OpenCV定义的宏,用来读取图像的像素值,这一部分就是进行中心变换
for( j = 0; j < cy; j++ )
{
for( i = 0; i < cx; i++ )
{
//中心化,将整体份成四块进行对角交换
tmp13 = CV_IMAGE_ELEM( image_Re, double, j, i);
CV_IMAGE_ELEM( image_Re, double, j, i) = CV_IMAGE_ELEM(image_Re, double, j+cy, i+cx);
CV_IMAGE_ELEM( image_Re, double, j+cy, i+cx) = tmp13;
tmp24 = CV_IMAGE_ELEM( image_Re, double, j, i+cx);
CV_IMAGE_ELEM( image_Re, double, j, i+cx) = CV_IMAGE_ELEM( image_Re, double, j+cy, i);
CV_IMAGE_ELEM( image_Re, double, j+cy, i) = tmp24;
}
}
//归一化处理将矩阵的元素值归一为[0,255]
//[(f(x,y)-minVal)/(maxVal-minVal)]*255
minVal = 0, maxVal = 0;
// Localize minimum and maximum values
cvMinMaxLoc( image_Re, &minVal, &maxVal );
// Normalize image (0 - 255) to be observed as an u8 image
scale = 255/(maxVal - minVal);
shift = -minVal * scale;
cvConvertScale(image_Re, dst, scale, shift);
printf("YES\n");
for (i = 0; i < nRow; i++)
{
for (j = 0; j < nCol; j++)
{
sum = sum + CV_IMAGE_ELEM(dst, double, i, j);
}
}
printf("YES\n");
hei = sum / nRow;
wid = sum / nCol;
for (i = 0; i < nRow; i++)
{
for (j = 0; j < nCol; j++)
{
meanx[i] = meanx[i] + CV_IMAGE_ELEM(dst, double, i, j);
meany[j] = meany[j] + CV_IMAGE_ELEM(dst, double, i, j);
}
}
printf("YES\n");
uchar* imagre = (uchar*)image_Re -> imageData;
uchar* imagim = (uchar*)image_Im -> imageData;
for (i = 2; i < nRow - 2; i++)
{
for (j = 2; j < nCol - 2; j++)
{
if (meany[j] < wid)
{
if (meanx[i] > hei)
{
if (i < nRow * 5 / 12 || i > nRow * 7 / 12)
{
for (wi = -2; wi <= 2; wi++)
{
for (wj = -2; wj <= 2; wj++)
{
sumwinr += imagre[(i + wi) * nCol + j + wj];
sumwini += imagim[(i + wi) * nCol + j + wj];
}
}
imagre[i * nCol + j] = sumwinr / 25;
imagim[i * nCol + j] = sumwini / 25;
}
}
}
}
}
for (i = 0; i < nRow; i++)
{
for (j = 0; j < nCol; j++)
{
image_Re -> imageData[i * nCol + j] = imagre[i * nCol + j];
image_Im -> imageData[i * nCol + j] = imagim[i * nCol + j];
}
}
for( j = 0; j < cy; j++ )
{
for( i = 0; i < cx; i++ )
{
//中心化,将整体份成四块进行对角交换
tmp13 = CV_IMAGE_ELEM( image_Re, double, j, i);
CV_IMAGE_ELEM( image_Re, double, j, i) = CV_IMAGE_ELEM(image_Re, double, j+cy, i+cx);
CV_IMAGE_ELEM( image_Re, double, j+cy, i+cx) = tmp13;
tmp24 = CV_IMAGE_ELEM( image_Re, double, j, i+cx);
CV_IMAGE_ELEM( image_Re, double, j, i+cx) = CV_IMAGE_ELEM( image_Re, double, j+cy, i);
CV_IMAGE_ELEM( image_Re, double, j+cy, i) = tmp24;
}
}
cvMerge(image_Re, image_Im, 0, 0, pTmp);
free(meanx);
meanx = NULL;
free(meany);
meany = NULL;
cvReleaseImage(&image_Re);
cvReleaseImage(&image_Im);
}
void Fft2Test()
{
//IplImage* tt = cvLoadImage("C:\\Users\\CDZHYF\\Desktop\\300.png",1);
//IplImage* ttt = cvCreateImage(cvGetSize(tt), tt->depth, 1);
IplImage *src = cvLoadImage("C:\\Users\\CDZHYF\\Desktop\\101.png",0);
IplImage *Fourier= cvCreateImage(cvGetSize(src),IPL_DEPTH_64F,2);
IplImage *dst = cvCreateImage(cvGetSize(src),IPL_DEPTH_64F,2);
IplImage *ImageRe= cvCreateImage(cvGetSize(src),IPL_DEPTH_64F,1);
IplImage *ImageIm= cvCreateImage(cvGetSize(src),IPL_DEPTH_64F,1);
IplImage *Image= cvCreateImage(cvGetSize(src),src->depth,src->nChannels);
IplImage *ImageDst= cvCreateImage(cvGetSize(src),src->depth,src->nChannels);
double m,M;
double scale;
double shift;
//cvCvtColor(tt, ttt, CV_BGR2GRAY);
//cvSaveImage("C:\\Users\\CDZHYF\\Desktop\\301.png", ttt);
//int height = src->height;
//int width = src->width;
// int lenth = height * width;
//int i, j;
//int his[256] = {0.0};
//int s[256] = {0.0};
//double p[256] = {0.0};
//double st[256] = {0.0};
Fft2(src,Fourier); //傅里叶变换
CalcFftFrequencySpectrum(Fourier, Image, Fourier); //计算FFT频率
cvDFT(Fourier,dst,CV_DXT_INV_SCALE);//实现傅里叶逆变换,并对结果进行缩放
cvSplit(dst,ImageRe,ImageIm,0,0);
//对数组每个元素平方并存储在第二个参数中
cvPow(ImageRe,ImageRe,2);
cvPow(ImageIm,ImageIm,2);
cvAdd(ImageRe,ImageIm,ImageRe,NULL);
cvPow(ImageRe,ImageRe,0.5);
cvMinMaxLoc(ImageRe,&m,&M,NULL,NULL);
scale = 255/(M - m);
shift = -m * scale;
//将shift加在ImageRe各元素按比例缩放的结果上,存储为ImageDst
cvConvertScale(ImageRe,ImageDst,scale,shift);
//for (i = 0; i < ImageDst->height * ImageDst->width; i++)
//{
// int t = ImageDst->imageData[i];
// his[t]++;
//}
//for (i = 0; i < 256; i++)
//{
// p[i] = (double)his[i] / (double)(height*width);
//}
//st[0] = p[0];
//for (i = 1; i < 256; i++)
//{
// st[i] = st[i - 1] + p[i];
//}
//for (i = 0; i < 256; i++)
//{
// s[i] = (int)((256 - 1) * st[i] +0.5);
// //printf("%d %d \n", i, s[i]);
//}
memset(data_out, 0, WIDTH_IN*HEIGHT_IN*sizeof(int));
//for (i = 0; i < height*width; i++)
//{
// int t = ImageDst->imageData[i];
// //printf("t: %d \n", t);
// ImageDst->imageData[i] = s[t];
//}
cvNamedWindow("源图像",0);
cvNamedWindow("傅里叶谱",0);
cvNamedWindow("傅里叶逆变换",0);
cvShowImage("源图像",src);
cvShowImage("傅里叶谱",Image);
cvShowImage("傅里叶逆变换",ImageDst);
cvWaitKey(0);
cvReleaseImage(&src);
cvReleaseImage(&Image);
cvReleaseImage(&ImageIm);
cvReleaseImage(&ImageRe);
cvReleaseImage(&Fourier);
cvReleaseImage(&dst);
cvReleaseImage(&ImageDst);
}
void LogTest()
{
cv::Mat img = cv::imread("C:\\Users\\CDZHYF\\Desktop\\101.png", CV_LOAD_IMAGE_GRAYSCALE);
if (img.empty())
{
printf("img load error \n");
system("pause");
exit(-1);
}
cv::Mat doubleImg;
img.convertTo(doubleImg, CV_64FC1);
doubleImg += cv::Scalar(1,1,1);
double minVal = 0;
double maxVal = 0;
cv::minMaxLoc(doubleImg, &minVal, &maxVal);
printf("doubleImg: minVal: %f, maxVal: %f \n", minVal, maxVal);
cv::Mat logImg;
cv::log(doubleImg, logImg);
minVal = 0;
maxVal = 0;
cv::minMaxLoc(logImg, &minVal, &maxVal);
printf("logImg: minVal: %f, maxVal: %f \n", minVal, maxVal);
double avgVal = cv::mean(logImg).val[0];
printf("avgVal: %f \n", avgVal);
double tVal = 0;
for (int y=0; y<logImg.rows; y++)
{
for (int x=0; x<logImg.cols; x++)
{
tVal = logImg.at<double>(y,x);
if (tVal < avgVal)
{
tVal += 1;
}
else
{
tVal -= 1;
}
logImg.at<double>(y,x) = tVal;
}
}
cv::Mat expImg;
cv::exp(logImg, expImg);
minVal = 0;
maxVal = 0;
cv::minMaxLoc(expImg, &minVal, &maxVal);
printf("expImg: minVal: %f, maxVal: %f \n", minVal, maxVal);
double scale = 255/(maxVal - minVal);
double shift = -minVal * scale;
//将shift加在ImageRe各元素按比例缩放的结果上,存储为ImageDst
//cvConvertScale(expImg,resImg,scale,shift);
cv::Mat resImg;
cv::convertScaleAbs(expImg, resImg, scale, shift);
//expImg.convertTo(resImg, CV_8UC1);
cv::imshow("img", img);
cv::imshow("exp", resImg);
cv::waitKey(0);
}
int main(int argc, char ** argv)
{
//FftTest();
Fft2Test();
//LogTest();
system("pause");
return 0;
}
fft代码
最新推荐文章于 2023-10-15 22:56:05 发布