fft代码

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