fft2.c

数字图像处理二维傅里叶变换代码

/*****************************************************************************/
/* 文件名：      FFT.c                                                       */
/* 程序语言：    C                                                           */
/* 功能：        Bmp图像（实验三）程序——傅立叶变换                         */
/* 说明：        无                                                          */
/* 完成日期：    2007-7-2                                                    */
/* 修改日期：    2007-7-2                                                    */
/* 作者：        李磊                                                        */
/* 版本：        1.0                                                         */
/*                                                                           */
/*  Copyright (c) 2007, Neural Network and Image Recognition Research Group, */
/*                  Beijing University of Technology,                        */
/*                        All rights reserved.                               */
/*****************************************************************************/


#include	"BmpImageIo.h"
#include	"BmpImageProcess.h"
#include    <stdio.h>
#include	<math.h>

#define		PI				3.1415926
#define		MAX				512

#define		FORWARD			1
#define		INVERSE			0

#define		FFTSHIFT		1
#define		NON_FFTSHIFT	0

typedef struct DCOMPLEX
{
	double r;
	double i;
}dcomplex;



//建立二维dcomplex类型数组，大小：dim1 * dim2
void newDcomplexMatrix(size_t dim1, size_t dim2, dcomplex ***pMatrix)
{
	size_t i;
	*pMatrix = (dcomplex **)malloc(dim1 * sizeof(dcomplex*));
	for(i = 0; i < dim1; ++i)
		(*pMatrix)[i] = (dcomplex*)malloc(dim2 * sizeof(dcomplex));
}


//删除二维dcomplex类型数组
void deleteDcomplexMatrix(size_t dim1, dcomplex*** pMatrix)
{
	size_t i;
	for(i = 0; i < dim1; ++i)
		free((*pMatrix)[i]);
	free(*pMatrix);
	*pMatrix = NULL;
}



//计算一维DFT，输入序列a[n]为dcomplex数组，长度为N，输出覆盖到原序列上
void DFT(dcomplex a[], int N, int flag)
{
    int j, k, sign;
    dcomplex temp[MAX];
	double alpha, COS, SIN;

    if(flag == 1) sign = -1;		// 正变换
        else sign = 1;				// 反变换

    for(j=0; j<N; j++)
        temp[j] = a[j];

    for(k=0; k<N; k++)
    {
        a[k].r = a[k].i =0;
        for(j=0; j<N; j++)
        {
			alpha = 2*PI*k*j/N;
			COS = cos(alpha);
			SIN = sin(alpha);
            a[k].r += (temp[j].r * COS - sign * temp[j].i * SIN);
            a[k].i += (sign * temp[j].r * SIN + temp[j].i * COS);
        }
    }

    if(sign==1)
    {
        for(k=0; k<N; k++)
        {
            a[k].r /= N;
            a[k].i /= N;
        }
    }

    return;
}


//统计能量分布
void EnergyStat(int Height, int Width, dcomplex** pImage)
{
    
	int j,k;
    double r,temp;
	double energy[7]={0};

	for(j=0; j<Height; j++)
	{
		for(k=0; k<Width; k++)
		{
			 temp = ((pImage[j][k].r)/64) * ((pImage[j][k].r)/64);
			 energy[0] += temp;
			 r = sqrt((j - (int)Height/2) * (j - (int)Height/2) + (k - (int)Width/2) * (k - (int)Width/2));
			
			if(r <= 80)
			{
			 	energy[6] += temp;
			   	if(r <= 50)
			   	{
				 	energy[5] += temp;
			      	if(r <= 35)
				  	{
					 	energy[4] += temp;
					 	if(r <= 25)
					 	{
							energy[3] += temp;
					    	if(r <= 15)
							{
								energy[2] += temp;
						   		if(r <= 5)
								{
									energy[1] += temp;
								} 
					 		}
				  		}
			   		}
				}
			}
	  	}
	}
	
	printf("\n能量分布统计结果:\n");
	printf("半径小于 5范围内的能量:  %.2f\n", 100 * energy[1] / energy[0]);
	printf("半径小于15范围内的能量:  %.2f\n", 100 * energy[2] / energy[0]);
	printf("半径小于25范围内的能量:  %.2f\n", 100 * energy[3] / energy[0]);
	printf("半径小于35范围内的能量:  %.2f\n", 100 * energy[4] / energy[0]);
	printf("半径小于50范围内的能量:  %.2f\n", 100 * energy[5] / energy[0]);
	printf("半径小于80范围内的能量:  %.2f\n", 100 * energy[6] / energy[0]);
	
	/*
	printf("\n能量分布统计结果:\n");
	printf("半径小于 5范围内的能量:  0.902\n");
	printf("半径小于15范围内的能量:  0.911\n");
	printf("半径小于25范围内的能量:  0.926\n");
	printf("半径小于35范围内的能量:  0.938\n");
	printf("半径小于50范围内的能量:  0.951\n");
	printf("半径小于80范围内的能量:  0.966\n");
	*/
}


//计算二维DFT
void DFT2(BmpImage *pOriginalImage, int Flag_of_Forwd_or_Inv, int Flag_of_FFTSHIFT)
{
	int j, k;
	double tmp;
	
	dcomplex** dImage;			// 复图象结构体，用于保存计算结果
	dcomplex*  dTemp;			// 复序列数组，用于计算一维DFT

	newDcomplexMatrix(MAX, MAX, &dImage);		//分配空间
	dTemp = (dcomplex*)malloc(MAX * sizeof(dcomplex));

	if(Flag_of_FFTSHIFT == FFTSHIFT && Flag_of_Forwd_or_Inv == FORWARD)
	{
		for(j=1; j<(int)pOriginalImage->height; j+=2)		//把image的内容乘以(-1)^(j+k)
		{
			for(k=0; k<(int)pOriginalImage->width; k+=2)
			{
				pOriginalImage->pB[j][k] = -(pOriginalImage->pB[j][k]);
			}
		}

		for(j=0; j<(int)pOriginalImage->height; j+=2)			
		{
			for(k=1; k<(int)pOriginalImage->width; k+=2)
			{
				pOriginalImage->pB[j][k] = -(pOriginalImage->pB[j][k]);
			}
		}
	}

	for(j=0; j<(int)pOriginalImage->height; j++)		//把image的内容复制到dImage实部
	{
		for(k=0; k<(int)pOriginalImage->width; k++)
		{
			dImage[j][k].r = ((double)pOriginalImage->pB[j][k]) / 256;
			dImage[j][k].i = 0.00;
		//	dImage[j][k].r = 256.00;
		//	dImage[j][k].i = ((double)pOriginalImage->pB[j][k])*256;
		}
	}

	/*********** 2-D FFT 开始 *************/
		for(j=0; j<(int)pOriginalImage->height; j++)		//对行数据做一维DFT
		{
			DFT(dImage[j], (int)pOriginalImage->width, Flag_of_Forwd_or_Inv);
		}

		for(k=0; k<(int)pOriginalImage->width; k++)			//对列数据做一维DFT
		{
			for(j=0; j<(int)pOriginalImage->height; j++)
			{
				dTemp[j] = dImage[j][k];
			}
			
			DFT(dTemp, (int)pOriginalImage->height, Flag_of_Forwd_or_Inv);

			for(j=0; j<(int)pOriginalImage->height; j++)
			{
				dImage[j][k] = dTemp[j];
			}
		}
	/*********** 2-D FFT 结束 *************/

	//统计能量分布
	//EnergyStat((int)pOriginalImage->height, (int)pOriginalImage->width, dImage);

	for(j=0; j<(int)pOriginalImage->height; j++)		//生成幅度谱
	{
		for(k=0; k<(int)pOriginalImage->width; k++)
		{
			tmp = sqrt(dImage[j][k].r * dImage[j][k].r + dImage[j][k].i * dImage[j][k].i);
			pOriginalImage->pB[j][k] = (BYTE)(tmp);
			pOriginalImage->pG[j][k] = pOriginalImage->pB[j][k];
			pOriginalImage->pR[j][k] = pOriginalImage->pB[j][k];
		}
	}


	if(Flag_of_Forwd_or_Inv == FORWARD)
	{
		if(Flag_of_FFTSHIFT == NON_FFTSHIFT)
		{
			writeBmp("FFT_Amp_NonShift.bmp", pOriginalImage);	//保存幅度谱图象
		}
		else if(Flag_of_FFTSHIFT == FFTSHIFT)
		{
			writeBmp("FFT_Amp_SHIFT.bmp", pOriginalImage);
		}
	}
	else if (Flag_of_Forwd_or_Inv == INVERSE)
	{
		writeBmp("IFFT.bmp", pOriginalImage);	
	}
	
	if(Flag_of_Forwd_or_Inv == FORWARD)
	{
		for(j=0; j<(int)pOriginalImage->height; j++)		//生成相位谱
		{
			for(k=0; k<(int)pOriginalImage->width; k++)
			{
				tmp = fabs(atan2(dImage[j][k].i, dImage[j][k].r));
				pOriginalImage->pB[j][k] = (BYTE)(tmp * 128);
				pOriginalImage->pG[j][k] = pOriginalImage->pB[j][k];
				pOriginalImage->pR[j][k] = pOriginalImage->pB[j][k];
			}
		}
	
		if(Flag_of_FFTSHIFT == NON_FFTSHIFT)
		{
			writeBmp("FFT_Phase_NonShift.bmp", pOriginalImage);			//保存相位谱图象
		}
		else if(Flag_of_FFTSHIFT == FFTSHIFT)
		{
			writeBmp("FFT_Phase_SHIFT.bmp", pOriginalImage);
		}
	}

	//删除临时变量的内存空间
	deleteDcomplexMatrix(pOriginalImage->height, &dImage);
	free(dTemp);

	return;
}


void main()
{
	BmpImage    image;				//建立BmpImage类型数据

    //初始化数据
    initialBmpImage(&image);
	
    //读图像
    readBmp("img12.bmp", &image);

	//做二维DFT
	DFT2(&image, FORWARD, NON_FFTSHIFT);

	//readBmp("MATLAB_FFT_Amp_SHIFT.bmp", &image);
	//EnergyStat(&image);

    //释放图像信息所占内存空间
    BmpImageDeleteRgbGray(&image);
  

}