www.pudn.com > EZW_.rar > wtransform.c

/*---------------------------------------------------------------------------*/ 
/*---------------------------------------------------------------------------*/ 
/* Port to C from C++ 
 * 
 * Mow-Song, Ng 2/9/2002 
 * msng@mmu.edu.my 
 * http://www.pesona.mmu.edu.my/~msng 
 * 
 * I do not claim copyright to the code, but if you use them or modify them, 
 * please drop me a mail. 
 * 
 */ 
/*---------------------------------------------------------------------------*/ 
/*---------------------------------------------------------------------------*/ 
 
/*---------------------------------------------------------------------------*/ 
/*---------------------------------------------------------------------------*/ 
/* Origianl copyright info */ 
/*---------------------------------------------------------------------------*/ 
// Baseline Wavelet Transform Coder Construction Kit 
// 
// Geoff Davis 
// gdavis@cs.dartmouth.edu 
// http://www.cs.dartmouth.edu/~gdavis 
// 
// Copyright 1996 Geoff Davis 9/11/96 
// 
// Permission is granted to use this software for research purposes as 
// long as this notice stays attached to this software. 
// 
/*---------------------------------------------------------------------------*/ 
/*---------------------------------------------------------------------------*/ 
 
#include "wtransform.h" 
 
 
/*---------------------------------------------------------------------------*/ 
/*---------------------------------------------------------------------------*/ 
WTRANSFORM *WaveletTransformAlloc(WAVELET *wavelet, FIMAGE *image, int nsteps,  
											 int symmetric) 
{ 
	WTRANSFORM *WTransform; 
 
	if ((WTransform=(WTRANSFORM *)malloc(sizeof(WTRANSFORM)))==NULL){ 
		return NULL; 
	} 
	 
	WTransform->wavelet=wavelet; 
	WTransform->nsteps=nsteps; 
	WTransform->symmetric=symmetric; 
 
	WTransform->value=NULL; 
 
	if (image!=NULL){ 
		WTransform->hsize=image->xsize; 
		WTransform->vsize=image->ysize; 
		WaveletTransformForward(WTransform, image, wavelet, nsteps, symmetric); 
	} 
	else{ 
		WTransform->hsize=WTransform->vsize=0; 
	} 
	return WTransform; 
} 
 
/*---------------------------------------------------------------------------*/ 
/*---------------------------------------------------------------------------*/ 
WTRANSFORM *WaveletTransformAllocBlank(WAVELET *wavelet,  
								int hsize, int vsize, int nsteps, int symmetric) 
{	 
	//int i; 
	WTRANSFORM *WTransform; 
 
	if ((WTransform=(WTRANSFORM *)malloc(sizeof(WTRANSFORM)))==NULL){ 
		return NULL; 
	} 
	 
	WTransform->hsize=hsize; 
	WTransform->vsize=vsize; 
	WTransform->wavelet=wavelet; 
	WTransform->nsteps=nsteps; 
	WTransform->symmetric=symmetric; 
 
	//- Strange, now we set them back to 0 & -1 
	//WTransform->nsteps=0; 
	//WTransform->symmetric=-1; 
	WaveletTransformInit(WTransform); 
	 
	//for (i=0; ivalue[i]=0; 
	//} 
 
	return WTransform; 
} 
 
 
/*---------------------------------------------------------------------------*/ 
/*---------------------------------------------------------------------------*/ 
WTRANSFORM *WaveletTransformAllocCopy(WTRANSFORM *WTransformSrc) 
{ 
	int i; 
	WTRANSFORM *WTransform; 
 
	if ((WTransform=(WTRANSFORM *)malloc(sizeof(WTRANSFORM)))==NULL){ 
		return NULL; 
	} 
	 
	WTransform->wavelet = WTransformSrc->wavelet; 
	WTransform->hsize = WTransformSrc->hsize; 
	WTransform->vsize = WTransformSrc->vsize; 
	WTransform->nsteps = WTransformSrc->nsteps; 
	WTransform->symmetric = WTransformSrc->symmetric; 
	 
	if (WTransformSrc->value == NULL) { 
		WTransform->value = NULL; 
	}  
	else { 
		WaveletTransformInit(WTransform); 
		for (i = 0; i < WTransform->hsize*WTransform->vsize; i++) 
			WTransform->value[i] = WTransformSrc->value[i]; 
	} 
	 
	return WTransform; 
 
} 
 
/*---------------------------------------------------------------------------*/ 
/*---------------------------------------------------------------------------*/ 
void WaveletTransformDealloc(WTRANSFORM *WTransform) 
{ 
	WaveletTransformFreeAll(WTransform); 
	free(WTransform); 
} 
 
/*---------------------------------------------------------------------------*/ 
/*---------------------------------------------------------------------------*/ 
Real WaveletTransformGetValue(WTRANSFORM *WaveletTransform, int scale,  
										int orientation, int x, int y) 
{ 
	int idx; 
	if (scale==0){ 
		/* special case */ 
		idx = 0; 
	} 
	else{ 
		idx = 3*scale - 2 + orientation; 
	} 
	 
	return (WaveletTransform->subbandPtr[idx] 
			[y*WaveletTransform->subbandHSize[idx] + x]); 
} 
 
/*---------------------------------------------------------------------------*/ 
/*---------------------------------------------------------------------------*/ 
void WaveletTransformSetValue(WTRANSFORM *WaveletTransform, int scale,  
										int orientation, int x, int y, Real val) 
{ 
	int idx; 
	if (scale==0){ 
		/* special case */ 
		idx = 0; 
	} 
	else{ 
		idx = 3*scale - 2 + orientation; 
	} 
	 
	WaveletTransform->subbandPtr[idx] 
			[y*WaveletTransform->subbandHSize[idx] + x] = val; 
	 
	return; 
} 
 
/*---------------------------------------------------------------------------*/ 
/*---------------------------------------------------------------------------*/ 
void WaveletTransformForward(WTRANSFORM *WTransform, FIMAGE *image,  
									  WAVELET *wavelet, int nsteps, int symmetric) 
{ 
	Real *temp; 
 
	// clear out old info and set up subband pointers 
	WaveletTransformFreeAll(WTransform); 
	WTransform->hsize = image->xsize; 
	WTransform->vsize = image->ysize; 
	WTransform->wavelet = wavelet; 
	WTransform->nsteps = nsteps; 
	WTransform->symmetric = symmetric; 
	WaveletTransformInit (WTransform); 
	 
	temp = (Real *)calloc(WTransform->hsize*WTransform->vsize, sizeof(Real)); 
 
	WaveletTransform2D(WTransform->wavelet, image->pixelLinear,  
		temp, WTransform->hsize, WTransform->vsize, 
		WTransform->nsteps, WTransform->symmetric); 
 
	// linearize data 
	MallatToLinear(WTransform, temp); 
	free(temp); 
 
} 
 
 
/*---------------------------------------------------------------------------*/ 
/*---------------------------------------------------------------------------*/ 
void WaveletTransformInvert(WTRANSFORM *WTransform, FIMAGE *invertedImage) 
{ 
	Real *temp; 
	 
	temp = (Real *)calloc(WTransform->hsize*WTransform->vsize, sizeof(Real)); 
	 
	/* put data in Mallat format*/ 
	/* note that the transform->value is always in linear format */ 
	LinearToMallat (WTransform, temp); 
	 
	WaveletInvert2D(WTransform->wavelet, temp, invertedImage->pixelLinear, 
		WTransform->hsize, WTransform->vsize, WTransform->nsteps, WTransform->symmetric); 
	free(temp); 
} 
 
/*---------------------------------------------------------------------------*/ 
/*---------------------------------------------------------------------------*/ 
FIMAGE *WaveletTransformInvertThis(WTRANSFORM *WTransform) 
{ 
	FIMAGE *image; 
	 
	image = FImageAlloc(WTransform->hsize, WTransform->vsize); 
 
   WaveletTransformInvert(WTransform, image); 
 
   return image; 
} 
 
 
/*---------------------------------------------------------------------------*/ 
/*---------------------------------------------------------------------------*/ 
void MallatToLinear(WTRANSFORM *WTransform, Real *mallat) 
{ 
	int i, j, k; 
	int *lowHSize, *lowVSize; 
	 
   /* arrays of top left corner coordinates for subbands */ 
	lowHSize = (int *) calloc(WTransform->nsteps, sizeof(int)); 
	lowVSize = (int *) calloc(WTransform->nsteps, sizeof(int)); 
	 
	/* highsed scale size */ 
	lowHSize[WTransform->nsteps-1] = (WTransform->hsize+1)/2; 
	lowVSize[WTransform->nsteps-1] = (WTransform->vsize+1)/2; 
	 
	for (i = WTransform->nsteps-2; i >= 0; i--) { 
		lowHSize[i] = (lowHSize[i+1]+1)/2; 
		lowVSize[i] = (lowVSize[i+1]+1)/2; 
	} 
	 
	// move transformed image (in Mallat order) into linear array structure 
	// special case for LL subband 
	for (j = 0; j < WTransform->subbandVSize[0]; j++){ 
		for (i = 0; i < WTransform->subbandHSize[0]; i++){ 
			WTransform->subbandPtr[0][j*WTransform->subbandHSize[0]+i]  
            = mallat[j*WTransform->hsize+i]; 
		} 
	} 
		 
	for (k = 0; k < WTransform->nsteps; k++) { 
		for (j = 0; j < WTransform->subbandVSize[k*3+1]; j++){ 
			for (i = 0; i < WTransform->subbandHSize[k*3+1]; i++){ 
					WTransform->subbandPtr[k*3+1][j*WTransform->subbandHSize[k*3+1]+i] =  
					mallat[j*WTransform->hsize+(lowHSize[k]+i)]; 
			} 
		} 
		 
		for (j = 0; j < WTransform->subbandVSize[k*3+2]; j++){ 
			for (i = 0; i < WTransform->subbandHSize[k*3+2]; i++){ 
				WTransform->subbandPtr[k*3+2][j*WTransform->subbandHSize[k*3+2]+i] =  
						mallat[(lowVSize[k]+j)*WTransform->hsize+i]; 
			} 
		} 
		 
		for (j = 0; j < WTransform->subbandVSize[k*3+3]; j++){ 
			for (i = 0; i < WTransform->subbandHSize[k*3+3]; i++){ 
				WTransform->subbandPtr[k*3+3][j*WTransform->subbandHSize[k*3+3]+i] =  
					mallat[(lowVSize[k]+j)*WTransform->hsize+(lowHSize[k]+i)]; 
			} 
		} 
	} 
		 
	free(lowHSize); 
	free(lowVSize); 
} 
 
 
 
/*---------------------------------------------------------------------------*/ 
/*---------------------------------------------------------------------------*/ 
void LinearToMallat(WTRANSFORM *WTransform, Real *mallat) 
{ 
	int i, j, k; 
	int *lowHSize, *lowVSize; 
	 
	lowHSize = (int *) calloc(WTransform->nsteps, sizeof(int)); 
	lowVSize = (int *) calloc(WTransform->nsteps, sizeof(int)); 
	 
	lowHSize[WTransform->nsteps-1] = (WTransform->hsize+1)/2; 
	lowVSize[WTransform->nsteps-1] = (WTransform->vsize+1)/2; 
	 
	for (i = WTransform->nsteps-2; i >= 0; i--) { 
		lowHSize[i] = (lowHSize[i+1]+1)/2; 
		lowVSize[i] = (lowVSize[i+1]+1)/2; 
	} 
	 
	// put linearized image in Mallat format 
	// special case for LL subband 
	for (j = 0; j < WTransform->subbandVSize[0]; j++){ 
		for (i = 0; i < WTransform->subbandHSize[0]; i++){ 
			mallat[j*WTransform->hsize+i] = WTransform->subbandPtr[0][j*WTransform->subbandHSize[0]+i]; 
		} 
	} 
 
	for (k = 0; k < WTransform->nsteps; k++) { 
		 
		for (j = 0; j < WTransform->subbandVSize[k*3+1]; j++){ 
			for (i = 0; i < WTransform->subbandHSize[k*3+1]; i++){ 
				mallat[j*WTransform->hsize+(lowHSize[k]+i)] =  
					WTransform->subbandPtr[k*3+1][j*WTransform->subbandHSize[k*3+1]+i]; 
			} 
		} 
				 
		for (j = 0; j < WTransform->subbandVSize[k*3+2]; j++){ 
			for (i = 0; i < WTransform->subbandHSize[k*3+2]; i++){ 
				mallat[(lowVSize[k]+j)*WTransform->hsize+i] =  
					WTransform->subbandPtr[k*3+2][j*WTransform->subbandHSize[k*3+2]+i]; 
			} 
		} 
		 
		for (j = 0; j < WTransform->subbandVSize[k*3+3]; j++){ 
			for (i = 0; i < WTransform->subbandHSize[k*3+3]; i++){ 
				mallat[(lowVSize[k]+j)*WTransform->hsize+(lowHSize[k]+i)] =  
					WTransform->subbandPtr[k*3+3][j*WTransform->subbandHSize[k*3+3]+i]; 
			} 
		} 
	} 
 
	free(lowHSize); 
	free(lowVSize); 
 
} 
 
 
/*---------------------------------------------------------------------------*/ 
/*---------------------------------------------------------------------------*/ 
/* NOTE : Some bug due to free() exists. If nsteps=0, then we got problem. */ 
void WaveletTransformInit(WTRANSFORM *WTransform) 
{ 
	int i; 
	int *lowHSize, *lowVSize, *highHSize, *highVSize; 
 
	assert(WTransform->nsteps > 0); 
 
	WTransform->value=(Real *)calloc(WTransform->hsize*WTransform->vsize, sizeof(Real)); 
	// TO DO : error checking 
 
	WTransform->nSubbands = 3*WTransform->nsteps+1; 
	WTransform->subbandSize = (int *)calloc(WTransform->nSubbands, sizeof(int)); 
	WTransform->subbandHSize = (int *)calloc(WTransform->nSubbands, sizeof(int)); 
	WTransform->subbandVSize = (int *)calloc(WTransform->nSubbands, sizeof(int)); 
	WTransform->subbandPtr = (Real **)calloc(WTransform->nSubbands, sizeof(Real *)); 
 
	lowHSize = (int *)calloc(WTransform->nsteps, sizeof(int)); 
	lowVSize = (int *)calloc(WTransform->nsteps, sizeof(int)); 
	highHSize = (int *)calloc(WTransform->nsteps, sizeof(int)); 
	highVSize = (int *)calloc(WTransform->nsteps, sizeof(int)); 
 
	lowHSize[WTransform->nsteps-1] = (WTransform->hsize+1)/2; 
	lowVSize[WTransform->nsteps-1] = (WTransform->vsize+1)/2; 
	highHSize[WTransform->nsteps-1] = WTransform->hsize/2;  
	highVSize[WTransform->nsteps-1] = WTransform->vsize/2; 
	 
	/* set the sizes for each blocks */ 
	for (i = WTransform->nsteps-2; i >= 0; i--) { 
		lowHSize[i] = (lowHSize[i+1]+1)/2; 
		lowVSize[i] = (lowVSize[i+1]+1)/2; 
		highHSize[i] = lowHSize[i+1]/2; 
		highVSize[i] = lowVSize[i+1]/2;  
	} 
	 
	/* subband[0] */ 
	WTransform->subbandPtr[0] = WTransform->value; 
	WTransform->subbandHSize[0] = lowHSize[0]; 
	WTransform->subbandVSize[0] = lowVSize[0]; 
	WTransform->subbandSize[0] = WTransform->subbandHSize[0]*WTransform->subbandVSize[0]; 
	 
	/* all other high pass subbands, do each scale at once */ 
	for (i = 0; i < WTransform->nsteps; i++) { 
		WTransform->subbandHSize[3*i+1] = highHSize[i]; 
		WTransform->subbandVSize[3*i+1] = lowVSize[i]; 
		WTransform->subbandHSize[3*i+2] = lowHSize[i]; 
		WTransform->subbandVSize[3*i+2] = highVSize[i]; 
		WTransform->subbandHSize[3*i+3] = highHSize[i]; 
		WTransform->subbandVSize[3*i+3] = highVSize[i]; 
	} 
	 
	/* set the data pointers */ 
	for (i = 1; i < WTransform->nSubbands; i++) { 
		WTransform->subbandSize[i] = WTransform->subbandHSize[i]*WTransform->subbandVSize[i]; 
		WTransform->subbandPtr[i] = WTransform->subbandPtr[i-1] + WTransform->subbandSize[i-1]; 
	} 
	 
	/* release allocated memory */ 
	free(lowHSize); 
	free(lowVSize); 
	free(highHSize); 
	free(highVSize); 
 
} 
	 
 
/*---------------------------------------------------------------------------*/ 
/*---------------------------------------------------------------------------*/ 
void WaveletTransformFreeAll(WTRANSFORM *WTransform) 
{ 
	if (WTransform->value!=NULL){ 
		free(WTransform->value); 
		free(WTransform->subbandSize); 
		free(WTransform->subbandHSize); 
		free(WTransform->subbandVSize); 
		free(WTransform->subbandPtr); 
	} 
} 
 
 
/*---------------------------------------------------------------------------*/ 
/*---------------------------------------------------------------------------*/ 
/* Scale the coefficients in each subband so that we have values from 0..255 */ 
FIMAGE *WaveletTransformCodeCoeff(WTRANSFORM *WaveletTransform, int DrawLine) 
{ 
   FIMAGE *image; 
   int i, j, k, lowHSize, lowVSize, hLength, vLength; 
   Real sMax, sMin, sRange; 
   Real tempValue; 
   WTRANSFORM *wtransform; 
       
   /* allocate image */ 
   image=FImageAlloc(WaveletTransform->hsize, WaveletTransform->vsize); 
 
   if (image==NULL){    
      WaveletTransformWarning("Cannot allocate image.\n"); 
      return NULL; 
   } 
 
   wtransform=WaveletTransformAllocBlank(WaveletTransform->wavelet, WaveletTransform->hsize, 
      WaveletTransform->vsize, WaveletTransform->nsteps, WaveletTransform->symmetric); 
 
	k=0; 
   for (i=0; inSubbands; i++){ 
      /* go through each subband */ 
      sMax=-MaxReal; 
      sMin= MaxReal; 
 
      /* find dynamic range */ 
      for (j=0; jsubbandHSize[i]*WaveletTransform->subbandVSize[i]; j++){ 
         tempValue=WaveletTransform->subbandPtr[i][j]; 
         if (tempValue>sMax){ 
            sMax=tempValue; 
         } 
         else if (tempValuesubbandHSize[i]*WaveletTransform->subbandVSize[i]; j++){ 
         wtransform->value[k++]=RINT(255*(WaveletTransform->subbandPtr[i][j]-sMin)/sRange);       
      } 
   } 
 
   /* convert to mallat format */    
   LinearToMallat(wtransform, image->pixelLinear); 
 
 
	if (DrawLine){ 
		hLength=WaveletTransform->hsize; 
		vLength=WaveletTransform->vsize; 
		lowHSize= (WaveletTransform->hsize+1)/2; 
		lowVSize= (WaveletTransform->vsize+1)/2; 
	 
		for (k=0; knsteps;k++){ 
			/* draw the lines */ 
			for (j=0; jpixel[j][lowHSize]=0; 
			} 
 
			for (j=0; jpixel[lowVSize][j]=0; 
			} 
 
			/* next levels */ 
			hLength=(hLength+1)/2; 
			vLength=(vLength+1)/2; 
			lowHSize=(lowHSize+1)/2; 
			lowVSize=(lowVSize+1)/2; 
		} 
	} 
 
	WaveletTransformDealloc(wtransform); 
   return image; 
} 
 
/*---------------------------------------------------------------------------*/ 
/*---------------------------------------------------------------------------*/ 
void WaveletTransformWriteCoeffToText(WTRANSFORM *WaveletTransform, char *filename) 
{ 
	FILE *in; 
   int i, j; 
       
	if ( (in=fopen(filename, "w")) == NULL ){ 
		WaveletTransformWarning("Unable to open text file.\n"); 
		return ; 
	} 
 
	for (i=0; inSubbands; i++){ 
      /* go through each subband */ 
		fprintf(in, "Subband #%d\n", i); 
		fprintf(in, "---------------------------------------------------------------------------\n"); 
 
      for (j=0; jsubbandHSize[i]*WaveletTransform->subbandVSize[i]; j++){ 
			fprintf(in, "%-8d %lf\n", j, WaveletTransform->subbandPtr[i][j]);			 
		} 
		fprintf(in, "\n"); 
   } 
 
	fclose(in); 
	return; 
} 
 
/*---------------------------------------------------------------------------*/ 
/*---------------------------------------------------------------------------*/ 
void WaveletTransformError(char *fmt, ...) 
{ 
	va_list argptr; 
	 
	va_start( argptr, fmt ); 
	fprintf(stderr, "WaveletTransformError: " ); 
	vprintf( fmt, argptr ); 
	va_end( argptr ); 
	exit( -1 ); 
 
} 
 
 
/*----------------------------------------------------------------------------*/	 
/*----------------------------------------------------------------------------*/	 
void WaveletTransformWarning(char *fmt, ...) 
{ 
	va_list argptr; 
	 
	va_start( argptr, fmt ); 
	fprintf( stderr, "WaveletTransformWarning: " ); 
	vprintf( fmt, argptr ); 
	va_end( argptr ); 
}