// MyFourial.cpp: implementation of the CMyFourial class.
//
//////////////////////////////////////////////////////////////////////

#include "stdafx.h"
#include "PhotoGrammetry.h"
#include "MyFourial.h"

#ifdef _DEBUG
#undef THIS_FILE
static char THIS_FILE[]=__FILE__;
#define new DEBUG_NEW
#endif

//////////////////////////////////////////////////////////////////////
// Construction/Destruction
//////////////////////////////////////////////////////////////////////

CMyFourial::CMyFourial()
{

}

CMyFourial::~CMyFourial()
{

}

BOOL CMyFourial::Fourier(LPSTR lpDIBBits, LONG lWidth, LONG lHeight)
{
	
	unsigned char*	lpSrc;	// 指向源圖像的指針
	double	dTemp;	
	LONG	i;				// 循環變量
	LONG	j;	
	LONG	w;				// 進行付立葉變換的寬度（2的整數次方）
	LONG	h;				// 進行付立葉變換的高度（2的整數次方）	
	int		wp;
	int		hp;
		
	LONG	lLineBytes;		// 圖像每行的字節數
	lLineBytes = WIDTHBYTES(lWidth * 8);	// 計算圖像每行的字節數
	
	// 賦初值
	w = 1;
	h = 1;
	wp = 0;
	hp = 0;
	
	// 計算進行付立葉變換的寬度和高度（2的整數次方）
	while(w * 2 <= lWidth)
	{
		w *= 2;
		wp++;
	}	
	while(h * 2 <= lHeight)
	{
		h *= 2;
		hp++;
	}
	
	complex<double> *TD = new complex<double>[w * h];
	complex<double> *FD = new complex<double>[w * h];
		
	for(i = 0; i < h; i++)			// 行
	{		
		for(j = 0; j < w; j++)		// 列
		{
			// 指向DIB第i行，第j個象素的指針
			lpSrc = (unsigned char*)lpDIBBits + lLineBytes * (lHeight - 1 - i) + j;
			
			// 給時域賦值
			TD[j + w * i] = complex<double>(*(lpSrc), 0);
		}
	}
	
	for(i = 0; i < h; i++)
	{
		// 對y方向進行快速付立葉變換
		FFT(&TD[w * i], &FD[w * i], wp);
	}
	
	// 保存變換結果
	for(i = 0; i < h; i++)
	{
		for(j = 0; j < w; j++)
		{
			TD[i + h * j] = FD[j + w * i];
		}
	}
	
	for(i = 0; i < w; i++)
	{
		// 對x方向進行快速付立葉變換
		FFT(&TD[i * h], &FD[i * h], hp);
	}
		
	for(i = 0; i < h; i++)			// 行
	{		
		for(j = 0; j < w; j++)		// 列
		{
			// 計算頻譜
			dTemp = sqrt(FD[j * h + i].real() * FD[j * h + i].real() + 
				         FD[j * h + i].imag() * FD[j * h + i].imag()) / 100;
			if (dTemp > 255)
			{
				dTemp = 255;
			}
			
			// 指向DIB第(i<h/2 ? i+h/2 : i-h/2)行，第(j<w/2 ? j+w/2 : j-w/2)個
			// 象素的指針，此處不直接取i和j，是為了將變換后的原點移到中心
			// lpSrc = (unsigned char*)lpDIBBits + lLineBytes * (lHeight-1-i) + j;
			lpSrc = (unsigned char*)lpDIBBits + lLineBytes * 
				(lHeight - 1 - (i<h/2 ? i+h/2 : i-h/2)) + (j<w/2 ? j+w/2 : j-w/2);
			
			// 更新源圖像
			* (lpSrc) = (BYTE)(dTemp);
		}
	}
	
	delete TD;
	delete FD;

	return TRUE;
}

/*************************************************************************
 * 函數名稱：
 *   FFT()
 * 參數:
 *   complex<double> * TD	- 指向時域數組的指針
 *   complex<double> * FD	- 指向頻域數組的指針
 *   r						－2的冪數，即迭代次數
 * 返回值:
 *   無。
 * 說明:
 *   該函數用來實現快速付立葉變換。
 ************************************************************************/
VOID CMyFourial::FFT(complex<double> * TD, complex<double> * FD, int r)
{	
	const double PI = 3.1415926535;
	
	LONG	count;				// 付立葉變換點數	
	int		i,j,k;				// 循環變量
	int		bfsize,p;	
	double	angle;				// 角度	
	complex<double> *W,*X1,*X2,*X;
	
	count = 1 << r;				// 計算付立葉變換點數
	
	// 分配運算所需存儲器
	W  = new complex<double>[count / 2];
	X1 = new complex<double>[count];
	X2 = new complex<double>[count];
	
	// 計算加權系數
	for(i = 0; i < count / 2; i++)
	{
		angle = -i * PI * 2 / count;
		W[i] = complex<double> (cos(angle), sin(angle));
	}
	
	// 將時域點寫入X1
	memcpy(X1, TD, sizeof(complex<double>) * count);
	
	// 采用蝶形算法進行快速付立葉變換
	for(k = 0; k < r; k++)
	{
		for(j = 0; j < 1 << k; j++)
		{
			bfsize = 1 << (r-k);
			for(i = 0; i < bfsize / 2; i++)
			{
				p = j * bfsize;
				X2[i + p] = X1[i + p] + X1[i + p + bfsize / 2];
				X2[i + p + bfsize / 2] = (X1[i + p] - X1[i + p + bfsize / 2]) 
					* W[i * (1<<k)];
			}
		}
		X  = X1;
		X1 = X2;
		X2 = X;
	}
	
	// 重新排序
	for(j = 0; j < count; j++)
	{
		p = 0;
		for(i = 0; i < r; i++)
		{
			if (j&(1<<i))
			{
				p+=1<<(r-i-1);
			}
		}
		FD[j]=X1[p];
	}
	
	delete W;
	delete X1;
	delete X2;

}