dct.cpp

本文件是离散余弦的具体实现

#define EXTRA_NAME "@dcttransform."
#include "loadbmp.h"

#define Point(x,y) lpPoints[(x)+(y)*nWidth]

void FFT(COMPLEX * TD, COMPLEX * FD, int power)
{
	int count;
	int i,j,k,bfsize,p;
	double angle;
	COMPLEX *W,*X1,*X2,*X;

	count=1<<power;
	W=(COMPLEX *)malloc(sizeof(COMPLEX)*count/2);
	X1=(COMPLEX *)malloc(sizeof(COMPLEX)*count);
	X2=(COMPLEX *)malloc(sizeof(COMPLEX)*count);
	for(i=0;i<count/2;i++)
	{
		angle=-i*pi*2/count;
		W[i].re=cos(angle);
		W[i].im=sin(angle);
	}
	memcpy(X1,TD,sizeof(COMPLEX)*count);
	for(k=0;k<power;k++)
	{
		for(j=0;j<1<<k;j++)
		{
			bfsize=1<<(power-k);
			for(i=0;i<bfsize/2;i++)
			{
				p=j*bfsize;
				X2[i+p]=Add(X1[i+p],X1[i+p+bfsize/2]);
				X2[i+p+bfsize/2]=Mul(Sub(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<power;i++)
		{
			if (j&(1<<i)) p+=1<<(power-i-1);
		}
		FD[j]=X1[p];
	}
	free(W);
	free(X1);
	free(X2);
}

void DCT_Pre(double *f, double *F, int power)
{
	int i,count;
	COMPLEX *X;
	double s;

	count=1<<power;
	X=(COMPLEX *)malloc(sizeof(COMPLEX)*count*2);
	memset(X,0,sizeof(COMPLEX)*count*2);
	for(i=0;i<count;i++)
	{
		X[i].re=f[i];
	}
	FFT(X,X,power+1);
	s=1/sqrt(count);
	F[0]=X[0].re*s;
	s*=sqrt(2);
	for(i=1;i<count;i++)
	{
		F[i]=(X[i].re*cos(i*pi/(count*2))+X[i].im*sin(i*pi/(count*2)))*s;
	}
	free(X);
}

void Dct()
{
	int w=1,h=1,wp=0,hp=0;
	while(w*2<=nWidth)
	{
		w*=2;
		wp++;
	}
	while(h*2<=nHeight)
	{
		h*=2;
		hp++;
	}
	int x,y;
	BYTE *lpPoints=new BYTE[nWidth*nHeight];
	GetPoints(lpPoints);

	double *f=new double[w*h];
	double *W=new double[w*h];

	for(y=0;y<h;y++)
	{
		for(x=0;x<w;x++)
		{
			f[x+y*w]=Point(x,y);
		}
	}

	for(y=0;y<h;y++)
	{
		DCT_Pre(&f[w*y],&W[w*y],wp);
	}
	for(y=0;y<h;y++)
	{
		for(x=0;x<w;x++)
		{
			f[x*h+y]=W[x+w*y];
		}
	}
	for(x=0;x<w;x++)
	{
		DCT_Pre(&f[x*h],&W[x*h],hp);
	}
	double a;
	memset(lpPoints,0,nWidth*nHeight);

	for(y=0;y<h;y++)
	{
		for(x=0;x<w;x++)
		{
			a=fabs(W[x*h+y]);
			if (a>255) a=255;
			Point(x,nHeight-y-1)=(BYTE)(a);
		}
	}
	delete f;
	delete W;
	PutPoints(lpPoints);
	delete lpPoints;
}

void main(int argc, char *argv[])
{
	if(argc==2)
		FileName=argv[1];
	else
		return;	
	OpenFile();
	Dct();
	SaveAs();
}