trans_d97.c

傅立叶变换和小波变换是图像压缩的重要工具。该代大戏是利用小波变换进行图像压缩。

/*********************

	I. 9 multiplies for every two pels pel, for each level
		this could be reduced by going explicility two-dimensional
		(since you can merge multiplies) but that's a real hassle
		9 * (1/2) * (10 cycles) * log2(4 levels) =
			100 clocks per pixel , is a lot of time!

		 the inverse is four multiplies per pel
		this hurts huge on my K6, not too bad on a P3 or K7
		<> could try 3dnow

	II. we're not quite normalized; eg. a 1 in transform space
			-> 0.96 , 1.02 or 1.08 depending on the bande
		(we compensate with the quantizers)

	III. should be able to do this fast in MMX; the wave coefficients are 16-bit
			R = (r[0] - r[4])
			C = coefficients on 0 - 4
		super-fast and inaccurate : result = R * C ; multiply, high words, & accumulate
			(one or a few instructions)
		more accurate is to multiply to 32 bits, accumulate, and then take the high word
		note that we need signs, which mmx is not good for...

***********************/

#include "transform.h"

#define F_EXTPAD 4
#define D_EXTPAD 2

void forwardDaub97(float *data,float *low,float *high,int half)
{
int i, N;
float *x;
static float wrapBuf[MAX_LINE_SIZE + F_EXTPAD*2 + 2];

	N = half + half;
	assert( N < MAX_LINE_SIZE );

	x = wrapBuf + F_EXTPAD;
	memcpy(x,data,sizeof(float)*N);
	for(i=1;i<=F_EXTPAD;i++)
	{
		x[-i] = x[i];
		x[(N-1) + i] = x[(N-1) - i];
	}

	while(half--)
	{
		*low++ =	0.0378285f * (x[ 4] + x[-4]) - 0.0238495f * (x[ 3] + x[-3]) - 0.110624f * (x[ 2] + x[-2]) + 
					0.3774030f * (x[ 1] + x[-1]) + 0.8526990f * x[0];

		x++;

		*high++ =  -0.064539f * (x[ 3] + x[-3]) + 0.040689f * (x[ 2] + x[-2]) + 
					0.418092f * (x[ 1] + x[-1]) - 0.788486f * x[ 0] ;

		x++;
	}
}

/***
void inverseDaub97(float *data,float *low,float *high,int half)
{
int i;
float *r, *d, *x;

	r = wrapBuf  + D_EXTPAD;  
	d = r + half + D_EXTPAD*2;
	memcpy(r,low ,half*sizeof(float));
	memcpy(d,high,half*sizeof(float));

	for(i=1;i<=D_EXTPAD;i++)
	{
		r[-i] = r[i];
		r[(half-1)+i] = r[half - i];
		d[-i] = d[i-1];
		d[(half-1)+i] = d[(half-1) - i];
	}
 
	x = data;
	while(half--)
	{
		x[0] = 0.788486f * r[0]				- 0.0406894f * (r[1] + r[-1])
		     - 0.023849f * (d[1] + d[-2])	+ 0.3774030f * (d[0] + d[-1]);

		x[1] = 0.418092f * (r[1] + r[0])  - 0.0645389f * (r[2] + r[-1])
		     - 0.037829f * (d[2] + d[-2]) + 0.1106240f * (d[1] + d[-1]) - 0.852699f * d[0];

		d++; r++;
		x += 2;
	}
}
*****/
void inverseDaub97(float *data,float *low,float *high,int half)
{

	assert( half >= 4 );
	half -= 4;

	// unroll to explicitly handle the four boundary pixels:

		data[0] = 0.788486f * low[0]				- 0.0406894f * (low[1] + low[1])
			     - 0.023849f * (high[1] + high[1])	+ 0.3774030f * (high[0] + high[0]);

		data[1] = 0.418092f * (low[1] + low[0])  - 0.0645389f * (low[2] + low[1])
			     - 0.037829f * (high[2] + high[1]) + 0.1106240f * (high[1] + high[0]) - 0.852699f * high[0];

		low++; high++; data += 2;

		data[0] = 0.788486f * low[0]				- 0.0406894f * (low[1] + low[-1])
			     - 0.023849f * (high[1] + high[-1])	+ 0.3774030f * (high[0] + high[-1]);

		data[1] = 0.418092f * (low[1] + low[0])  - 0.0645389f * (low[2] + low[-1])
			     - 0.037829f * (high[2] + high[-1]) + 0.1106240f * (high[1] + high[-1]) - 0.852699f * high[0];

		low++; high++; data += 2;

	while(half--)
	{
		data[0] = 0.788486f * low[0]				- 0.0406894f * (low[1] + low[-1])
			     - 0.023849f * (high[1] + high[-2])	+ 0.3774030f * (high[0] + high[-1]);

		data[1] = 0.418092f * (low[1] + low[0])  - 0.0645389f * (low[2] + low[-1])
			     - 0.037829f * (high[2] + high[-2]) + 0.1106240f * (high[1] + high[-1]) - 0.852699f * high[0];

		low++; high++; data += 2;
	}
	
		data[0] = 0.788486f * low[0]				- 0.0406894f * (low[1] + low[-1])
			     - 0.023849f * (high[1] + high[-2])	+ 0.3774030f * (high[0] + high[-1]);

		data[1] = 0.418092f * (low[1] + low[0])  - 0.0645389f * (low[1] + low[-1])
			     - 0.037829f * (high[0] + high[-2]) + 0.1106240f * (high[1] + high[-1]) - 0.852699f * high[0];

		low++; high++; data += 2;
		
		data[0] = 0.788486f * low[0]				- 0.0406894f * (low[0] + low[-1])
			     - 0.023849f * (high[-1] + high[-2])	+ 0.3774030f * (high[0] + high[-1]);

		data[1] = 0.418092f * (low[0] + low[0])  - 0.0645389f * (low[-1] + low[-1])
			     - 0.037829f * (high[-2] + high[-2]) + 0.1106240f * (high[-1] + high[-1]) - 0.852699f * high[0];

		low++; high++; data += 2;
}



Transform transformDaub97 =
{
	forwardDaub97,
	inverseDaub97,
	"Daub 97",
	4,5,
	NULL
};