wavelet.c

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

#include "wavelet.h"
#include "transform.h"
#include "quantize.h"
#include "coder.h"
#include "mathutil.h"

void Wavelet_SetupQuantizers(Wavelet *w,float baseQ)
{
int i,j,p;
SubBand * sb;
	
	w->quantizer = baseQ;
	w->dpcmLoss = (uint)((baseQ + 1) * 0.5f);

	// !! NOTEZ BIEN ! this normalization of the quantizers is critical!
	//	we've failed to do this in all previous wavers, which hurts the
	//		non-normalized stuff big time!
	//	of our transforms, only daub97 is normalized

	for(p=0;p<3;p++)
	{
	SubBandTree * sbt;
	int levels;

		sbt = w->trees + p;		

		levels = w->levels;
		if ( p < 2 )
			levels--;

		sb = &(sbt->LL);
		sb->quantizer = 1.0f / SubBand_TotalUntransformScale(sb,levels);

		for(i=0;i<3;i++)
		{
			for(j=0;j<levels;j++)
			{
				sb = &(sbt->bands[i][j]);

				sb->quantizer = baseQ / SubBand_TotalUntransformScale(sb,levels);
			}
		}
	}
}

Wavelet * Wavelet_Create(void * planes[3],int width,int height)
{
Wavelet * w;
int i,j,b;
SubBand * sb;

	w = new(Wavelet);
	assert(w);

	w->ari = arithInit();

	w->dpcmLoss = 0;
	w->quantizer = 0.0f;

	w->width = width;
	w->height = height;
	w->levels = chooseLevels(width,height);

	for(i=0;i<4;i++)
	{
		w->untransformedMSE[i] = untransformedMSE(i);
	}

	for(i=0;i<3;i++)
	{
		w->CodeBuffer[i] = RowBuffer_Create(width/2,CODEBUFFER_ROWS);
		assert(w->CodeBuffer[i]);
	}
	w->TransBufL = RowBuffer_Create(width,vTransform->bufRowsL);
	w->TransBufH = RowBuffer_Create(width,vTransform->bufRowsH);

	for(b=0;b<3;b++)
	{
	SubBandTree * sbt;

		sbt = w->trees + b;

		sbt->iplane = planes[b];
		sbt->fplane = NULL;

		for(i=0;i<3;i++)
		{
			sbt->bands[i] = calloc(1,sizeof(SubBand)*w->levels);
			assert(sbt->bands[i]);
		}

		// init the LL band :
		
		sb = &(sbt->LL);
		
		sb->type = 0;
		sb->wavelet = w;
		sb->stride = width;
		sb->width  = width  >> (w->levels);
		sb->height = height >> (w->levels);
		sb->band = sbt->iplane;
		sb->level = 0;

		if ( b < 2 )
		{
			// it's half the size!
			sb->stride >>= 1;
		}

		// init the other subbands :

		for(i=0;i<3;i++)
		{
			for(j=0;j<(w->levels);j++)
			{
				sb = &(sbt->bands[i][j]);
				
				if ( b < 2 && j == (w->levels - 1) )
					continue;

				sb->type = i + 1;
				sb->wavelet = w;
				sb->stride = width;

				if ( b < 2 )
				{
					// it's half the size!
					sb->stride >>= 1;
				}

				sb->width  = width  >> (w->levels - j);
				sb->height = height >> (w->levels - j);
				
				sb->level = j;

				switch(sb->type)
				{
					case SBT_LH:
						sb->band = sbt->iplane + sb->width;
						break;
					case SBT_HL:
						sb->band = sbt->iplane + sb->height * sb->stride;
						break;
					case SBT_HH:
						sb->band = sbt->iplane + sb->width + sb->height * sb->stride;
						break;
					default:
						assert(0);
						break;
				}
			}
		}
	}

return w;
}

void SubBand_DeInit(SubBand * sb)
{
	if ( sb->lastBand1 )
		free(sb->lastBand1);
	if ( sb->lastBand2 )
		free(sb->lastBand2);
}

void Wavelet_Destroy(Wavelet *w)
{
int i,b,l;

	for(b=0;b<3;b++)
	{
	SubBandTree * sbt;
		sbt = w->trees + b;

		SubBand_DeInit(&(sbt->LL));

		for(i=0;i<3;i++)
		{
			for(l=0;l<w->levels;l++)
			{
				SubBand_DeInit(&(sbt->bands[i][l]));
			}
			destroy(sbt->bands[i]);
		}

		if ( sbt->fplane )
			destroy(sbt->fplane);
	}

	for(i=0;i<3;i++)
	{
		RowBuffer_Destroy(w->CodeBuffer[i]);
	}
	RowBuffer_Destroy(w->TransBufL);
	RowBuffer_Destroy(w->TransBufH);

	if ( w->compAlloc )
		destroy(w->compAlloc);

	arithFree(w->ari);

	destroy(w);
}

#ifdef _DEBUG
#include "image.h"
#endif

void Wavelet_Transform(Wavelet *w)
{
int b;

	for(b=0;b<3;b++)
	{
	SubBandTree * sbt;
	int width,height,levels;

		sbt = w->trees + b;

		width  = w->width;
		height = w->height;
		levels = w->levels;

		if ( b < 2 )
		{
			width  >>= 1;
			height >>= 1;
			levels --;
		}

		if ( ! sbt->fplane )
			sbt->fplane = malloc(width * height * sizeof(float));
		memcpy(sbt->fplane,sbt->iplane,width * height * sizeof(float));

		transformPlane(sbt->fplane,width,height,levels,false);
	}
}

void Wavelet_Quantize(Wavelet *w,float q)
{
int i,l;
int b;
int width,height;

	Wavelet_SetupQuantizers(w,q);

	for(b=0;b<3;b++)
	{
	SubBandTree * sbt;
		sbt = w->trees + b;

		width = w->width;
		height = w->height;

		if ( b < 2 )
		{
			width  >>= 1;
			height >>= 1;
		}

		memcpy(sbt->iplane,sbt->fplane,width * height * sizeof(float));

		SubBand_Quantize(&(sbt->LL),true);

		for(l=0;l<(w->levels);l++)
		{
			if ( b < 2 && l == (w->levels - 1) )
			{
				continue;
			}
			for(i=0;i<3;i++)
			{
				SubBand_Quantize(&(sbt->bands[i][l]),true);
			}
		}
	}
}

void Wavelet_UnTransformAndQuantize(Wavelet *w)
{
int i,l,b;
int width,height,levels;

	for(b=0;b<3;b++)
	{
	SubBandTree * sbt;
		sbt = w->trees + b;

		SubBand_Quantize(&(sbt->LL),false);

		width  = w->width;
		height = w->height;
		levels = w->levels;

		if ( b < 2 )
		{
			width  >>= 1;
			height >>= 1;
			levels --;
		}

		for(l=0;l<levels;l++)
		{
			for(i=0;i<3;i++)
			{
				SubBand_Quantize(&(sbt->bands[i][l]),false);
			}
		}

		transformPlane((float *)sbt->iplane,width,height,levels,true);
	}
}

uint Wavelet_CompLen(Wavelet *w)
{
return w->compLen;
}

void Wavelet_WritePacketsToBuf(Wavelet *w,ubyte *buf)
{
	memcpy(buf, w->compPtr, w->compLen );
}

void Wavelet_CodeToPackets(Wavelet *w)
{
	if ( ! w->compAlloc )
	{
		w->compAlloc = malloc(w->width * w->height);
		w->compPtr = w->compAlloc + 4;
	}

	arithEncodeInit(w->ari,w->compPtr);

	codeWavelet(w,true);
	
	w->compLen = arithEncodeDoneMinimal(w->ari);
}

void Wavelet_DecodeFromBuf(Wavelet *w,ubyte *buf,float q)
{
	Wavelet_SetupQuantizers(w,q);

	w->compPtr = buf;

	arithDecodeInit(w->ari,w->compPtr);

	codeWavelet(w,false);
	
	arithDecodeDone(w->ari);
}

void SubBand_MakeLastBands(SubBand *sb)
{
	assert( sb );

	if ( ! sb->lastBand1 )
	{
		sb->lastBand1 = calloc(1,sb->width * sb->height * 4);
		sb->lastBand2 = calloc(1,sb->width * sb->height * 4);
	}
}

float SubBand_TotalUntransformScale(SubBand *sb,int levels)
{
float scale,LLscale;
int level;

	scale = 1.0f;
	level = levels - sb->level - 1;
	assert( level >= 0 );

	LLscale = sb->wavelet->untransformedMSE[SBT_LL];
	while( level--)
	{
		scale *= LLscale;	
	}
	scale *= sb->wavelet->untransformedMSE[sb->type];

return fsqrt(scale);
}

void SubBand_Quantize(SubBand *sb,bool forward)
{
	if ( sb->type == SBT_LL )
	{
		if ( forward )
		{
			doquantizenopack(sb->band,sb->width,sb->height,sb->stride,sb->quantizer);
		}
		else
		{
			dequantizenopack(sb->band,sb->width,sb->height,sb->stride,sb->quantizer);
		}
	}
	else
	{
		if ( forward )
		{
			doquantize(sb->band,sb->width,sb->height,sb->stride,sb->quantizer);
		}
		else
		{
			dequantize(sb->band,sb->width,sb->height,sb->stride,sb->quantizer);
		}
	}
}