transform.c

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

#include "transform.h"

/*}{********* the transforms **************/

extern Transform transformDaub97;
extern Transform transformCDF22;
extern Transform transformHaar;

/*********

Lena512 , with D97 H:
Vert	PSNR
Haar	35.65	(35.45 with cdf22 H)
CDF22	36.95	(36.59 with cdf22 H)
D97		37.17

*********/

#ifndef FAST_TRANSFORM //{ // high quality
const Transform * hTransform = &transformDaub97;
const Transform * vTransform = &transformCDF22;
#else //}{// high speed
const Transform * hTransform = &transformCDF22;
const Transform * vTransform = &transformHaar;
#endif //}

//const Transform * hTransform = &transformDaub97;
//const Transform * hTransform = &transformCDF22;
// const Transform * hTransform = &transformHaar;

//const Transform * vTransform = &transformHaar;
//const Transform * vTransform = &transformCDF22;
//const Transform * vTransform = &transformDaub97;

/*}{********* functions **************/

int chooseLevels(int w,int h)
{
int l;

	l = 0;

	for(;;)
	{
//		if ( (w&1) || (h&1) )
		if ( (w&3) || (h&3) ) // must be a multiple of 4 cuz we will subsamp the U&V !
			break;
		if ( w < 16 || h < 16 )
			break;

		l++;
		w >>= 1;
		h >>= 1;
	}

return l;
}

void transformPlane(float * plane,int width,int height,int levels,bool inverse)
{
int w,h,l;
transformFunc hTranser,vTranser;
float *data,*low,*high;
float *work,*column;

	l = max(width,height) + 9;
	work = malloc(2*l*sizeof(float));
	column = work + l;

	if ( ! inverse ) // forward
	{

		hTranser = hTransform->forward;
		vTranser = vTransform->forward;

		assert( hTranser && vTranser );

		for(l=0;l<levels;l++)
		{
		int x,y,half;

			w = (width >>l);
			h = (height>>l);

			//columns
			
			half = h>>1;
			low  = work;
			high = low + half;

			data = plane;
			for(x=0;x<w;x++)
			{
				for(y=0;y<h;y++)
					column[y] = data[y*width];
				vTranser(column,low,high,half);
				for(y=0;y<h;y++)
					data[y*width] = work[y];
				data ++;
			}

			//rows

			half = w>>1;
			low  = work;
			high = low + half;

			data = plane;
			for(y=0;y<h;y++)
			{
				hTranser(data,low,high,half);
				memcpy(data,work,w*4);
				data += width;
			}
		}
	}
	else // inverse
	{
		hTranser = hTransform->inverse;
		vTranser = vTransform->inverse;

		assert( hTranser && vTranser );

		for(l=levels-1;l>=0;l--)
		{
		int x,y,half;
		float *wlo,*whi;

			w = (width >>l);
			h = (height>>l);

			//rows

			half = w>>1;
			low  = work;
			high = low + half;
			wlo = work + 4;
			whi = column + 4;

			data = plane;
			for(y=0;y<h;y++)
			{
				memcpy(wlo,data     ,half*sizeof(float));
				memcpy(whi,data+half,half*sizeof(float));
				hTranser(data,wlo,whi,half);
				data += width;
			}

			//columns
			
			half = h>>1;
			low  = work;
			high = low + half;
			wlo = work + 4;
			whi = column + 4;

			data = plane;
			for(x=0;x<w;x++)
			{
				for(y=0;y<h;y++)
					work[y] = data[y*width];
				vTranser(column,low,high,half);
				for(y=0;y<h;y++)
					data[y*width] = column[y];
				data ++;
			}
		}
	}

	free(work);
}

//-----------------------------------------------------
// untransformedMSE
//		measure how much an error in band N affects the untransformed band

#define WORK_SIZE	(32)	// <> can be 16?
#define WORK_14		(WORK_SIZE/4)
#define WORK_34		(WORK_SIZE*3/4)

float	untransformedMSE(subBandType sbt)
{
float work[WORK_SIZE*WORK_SIZE];
float mse;
int i;

	memset(work,0,WORK_SIZE*WORK_SIZE*sizeof(float));

	switch(sbt)
	{
		case SBT_LL:
			work[ WORK_14 * WORK_SIZE + WORK_14 ] = 1.0f;
			break;
		case SBT_LH:
			work[ WORK_14 * WORK_SIZE + WORK_34 ] = 1.0f;
			break;
		case SBT_HL:
			work[ WORK_34 * WORK_SIZE + WORK_14 ] = 1.0f;
			break;
		case SBT_HH:
			work[ WORK_34 * WORK_SIZE + WORK_34 ] = 1.0f;
			break;
	}

	transformPlane(work,WORK_SIZE,WORK_SIZE,1,true);

	mse = 0.0f;
	for(i=0;i<(WORK_SIZE*WORK_SIZE);i++)
	{
		mse += work[i] * work[i];
	}

return mse;
}

/*}{**** EOF *******/