dwt.c

基于嵌入式零树小波编码的C语言程序实现。

	t->SubbandDir = (int *)calloc(t->nSubbands, sizeof(int));

	t->SubbandPtr = (REAL **)calloc(t->nSubbands, sizeof(REAL *));

	LowXSize = (int *)calloc(t->nLevels, sizeof(int));
	LowYSize = (int *)calloc(t->nLevels, sizeof(int));
	HighXSize = (int *)calloc(t->nLevels, sizeof(int));
	HighYSize = (int *)calloc(t->nLevels, sizeof(int));

	LowXSize[t->nLevels-1] = (t->XSize+1)/2;
	LowYSize[t->nLevels-1] = (t->YSize+1)/2;
	HighXSize[t->nLevels-1] = t->XSize/2; 
	HighYSize[t->nLevels-1] = t->YSize/2;
	
	/* set the sizes for each blocks */
	for (i = t->nLevels-2; i >= 0; i--) {
		LowXSize[i] = (LowXSize[i+1]+1)/2;
		LowYSize[i] = (LowYSize[i+1]+1)/2;
		HighXSize[i] = LowXSize[i+1]/2;
		HighYSize[i] = LowYSize[i+1]/2; 
	}
	
	/* subband[0] */
	t->SubbandPtr[0] = t->Coeff;
	t->SubbandXSize[0] = LowXSize[0];
	t->SubbandYSize[0] = LowYSize[0];
	t->SubbandSize[0] = t->SubbandXSize[0]*t->SubbandYSize[0];
	t->SubbandDir[0] = 0;
	t->SubbandScale[0] = 0;

	/* all other high pass subbands, do each scale at once */
	for (i = 0; i < t->nLevels; i++) {
		t->SubbandDir[3*i+1] = 0;
		t->SubbandScale[3*i+1] = i+1;
		t->SubbandXSize[3*i+1] = HighXSize[i];
		t->SubbandYSize[3*i+1] = LowYSize[i];

		t->SubbandDir[3*i+2] = 1;
		t->SubbandScale[3*i+2] = i+1;
		t->SubbandXSize[3*i+2] = LowXSize[i];
		t->SubbandYSize[3*i+2] = HighYSize[i];

		t->SubbandDir[3*i+3] = 2;
		t->SubbandScale[3*i+3] = i+1;
		t->SubbandXSize[3*i+3] = HighXSize[i];
		t->SubbandYSize[3*i+3] = HighYSize[i];
	}
	
	/* set the data pointers */
	for (i = 1; i < t->nSubbands; i++) {
		t->SubbandSize[i] = t->SubbandXSize[i]*t->SubbandYSize[i];
		t->SubbandPtr[i] = t->SubbandPtr[i-1] + t->SubbandSize[i-1];
	}
	
	/* release allocated memory */
	free(LowXSize);
	free(LowYSize);
	free(HighXSize);
	free(HighYSize);

	return t;
}

/*--------------------------------------------------------------------------------------------*/
void DWTransformDealloc(DWTransform *t)
{
	int i;

	if (t!=NULL){
		free(t->Coeff);
		free(t->SubbandSize);
		free(t->SubbandXSize);
		free(t->SubbandYSize);
		free(t->SubbandDir);
		free(t->SubbandScale);
		free(t->SubbandPtr);
		free(t);
	}
}

/*--------------------------------------------------------------------------------------------*/
REAL DWTransformGetValue(DWTransform *t, int scale, int dir, int x, int y)
{
	int idx;
	if (scale==0){
		idx = 0;
	}
	else{
		idx = 3*scale-2+dir;

	}
	return t->SubbandPtr[idx][y*t->SubbandXSize[idx]+x];
}

/*--------------------------------------------------------------------------------------------*/
void DWTransformSetValue(DWTransform *t, int scale, int dir, int x, int y, REAL val)
{
	int idx;
	if (scale==0){
		idx = 0;
	}
	else{
		idx = 3*scale-2+dir;

	}
	
	t->SubbandPtr[idx][y*t->SubbandXSize[idx]+x] = val;
	
}

/*--------------------------------------------------------------------------------------------*/
void DWTransformInputFromMallat(DWTransform *t, REAL *in)
{
	int i, j, k;
	int *LowXSize, *LowYSize;
	
   /* arrays of top left corner coordinates for subbands */
	LowXSize = (int *) calloc(t->nLevels, sizeof(int));
	LowYSize = (int *) calloc(t->nLevels, sizeof(int));
	
	/* highsed scale size */
	LowXSize[t->nLevels-1] = (t->XSize+1)/2;
	LowYSize[t->nLevels-1] = (t->YSize+1)/2;
	
	for (i = t->nLevels-2; i >= 0; i--) {
		LowXSize[i] = (LowXSize[i+1]+1)/2;
		LowYSize[i] = (LowYSize[i+1]+1)/2;
	}
	
	/* move transformed image (in Mallat order) into linear array structure */
	/* special case for LL subband */
	for (j = 0; j < t->SubbandYSize[0]; j++){
		for (i = 0; i < t->SubbandXSize[0]; i++){
			t->SubbandPtr[0][j*t->SubbandXSize[0]+i] = in[j*t->XSize+i];
		}
	}
		
	for (k = 0; k < t->nLevels; k++) {
		for (j = 0; j < t->SubbandYSize[3*k+1]; j++){
			for (i = 0; i < t->SubbandXSize[3*k+1]; i++){
				t->SubbandPtr[3*k+1][j*t->SubbandXSize[3*k+1]+i] =
					in[j*t->XSize+(LowXSize[k]+i)];
			}
		}
		
		for (j = 0; j < t->SubbandYSize[3*k+2]; j++){
			for (i = 0; i < t->SubbandXSize[3*k+2]; i++){
				t->SubbandPtr[3*k+2][j*t->SubbandXSize[3*k+2]+i] = 
					in[(LowYSize[k]+j)*t->XSize+i];
			}
		}
		
		for (j = 0; j < t->SubbandYSize[k*3+3]; j++){
			for (i = 0; i < t->SubbandXSize[k*3+3]; i++){
				t->SubbandPtr[3*k+3][j*t->SubbandXSize[k*3+3]+i] = 
					in[(LowYSize[k]+j)*t->XSize+(LowXSize[k]+i)];
			}
		}
	}
		
	free(LowXSize);
	free(LowYSize);

}

/*--------------------------------------------------------------------------------------------*/
void DWTransformOutputToMallat(DWTransform *t, REAL *out)
{
	int i, j, k;
	int *LowXSize, *LowYSize;
	
   /* arrays of top left corner coordinates for subbands */
	LowXSize = (int *) calloc(t->nLevels, sizeof(int));
	LowYSize = (int *) calloc(t->nLevels, sizeof(int));
	
	/* highsed scale size */
	LowXSize[t->nLevels-1] = (t->XSize+1)/2;
	LowYSize[t->nLevels-1] = (t->YSize+1)/2;
	
	for (i = t->nLevels-2; i >= 0; i--) {
		LowXSize[i] = (LowXSize[i+1]+1)/2;
		LowYSize[i] = (LowYSize[i+1]+1)/2;
	}
	
	/* put linearized image in Mallat format */
	/* special case for LL subband */
	for (j = 0; j < t->SubbandYSize[0]; j++){
		for (i = 0; i < t->SubbandXSize[0]; i++){
			out[j*t->XSize+i] = t->SubbandPtr[0][j*t->SubbandXSize[0]+i];
		}
	}

	for (k = 0; k < t->nLevels; k++) {
		
		for (j = 0; j < t->SubbandYSize[3*k+1]; j++){
			for (i = 0; i < t->SubbandXSize[3*k+1]; i++){
				out[j*t->XSize+(LowXSize[k]+i)] = 
					t->SubbandPtr[3*k+1][j*t->SubbandXSize[3*k+1]+i];
			}
		}
				
		for (j = 0; j < t->SubbandYSize[3*k+2]; j++){
			for (i = 0; i < t->SubbandXSize[3*k+2]; i++){
				out[(LowYSize[k]+j)*t->XSize+i] = 
					t->SubbandPtr[3*k+2][j*t->SubbandXSize[3*k+2]+i];
			}
		}
		
		for (j = 0; j < t->SubbandYSize[k*3+3]; j++){
			for (i = 0; i < t->SubbandXSize[k*3+3]; i++){
				out[(LowYSize[k]+j)*t->XSize+(LowXSize[k]+i)] = 
					t->SubbandPtr[3*k+3][j*t->SubbandXSize[3*k+3]+i];
			}
		}
	}

	free(LowXSize);
	free(LowYSize);
}

/*--------------------------------------------------------------------------------------------*/
void DWTransformScale(DWTransform *t, int scale_all)
{
	int i, n, s;

	if (scale_all){
		for (n=0; n<t->nSubbands; n++){
			s = (1<<(t->nLevels-t->SubbandScale[n]));
			for (i=0; i<t->SubbandSize[n]; i++){
				t->SubbandPtr[n][i]/= (Real)s;
			}	
		}
	}
	else{	
		s = (1<<(t->nLevels));
		for (i=0; i<t->SubbandSize[0]; i++){
			t->SubbandPtr[0][i]/= (Real)s;
		}	
	}
}

/*--------------------------------------------------------------------------------------------*/
void DWTransformError(char *fmt, ...)
{
	va_list argptr;
	
	va_start( argptr, fmt );
	fprintf(stderr, "DWTransformError: " );
	vprintf( fmt, argptr );
	va_end( argptr );
	exit( -1 );
}

/*--------------------------------------------------------------------------------------------*/
void DWTransformWarning(char *fmt, ...)
{
	va_list argptr;
	
	va_start( argptr, fmt );
	fprintf(stderr, "DWTransformWarning: " );
	vprintf( fmt, argptr );
	va_end( argptr );
	exit( -1 );
}

/*--------------------------------------------------------------------------------------------*/
/*--------------------------------------------------------------------------------------------*/

/*--------------------------------------------------------------------------------------------*/
#define AnalDelay97		4
#define AnalDelay137		6
#define AnalDelay53		2
#define AnalDelay1018	8
#define AnalDelayAdel   4

/*--------------------------------------------------------------------------------------------*/
/* perform 1-D subband analysis - odd length LP FIR filter */
void Subband1DAnal(REAL *in, int cols, int nLevels, int WaveletIndex)
{
	int c;
	int i, j;
	int nc;
	REAL *ImgPtr, *LinePtr, *TempPtr;
	void (*Anal)(REAL *, int);

	nc = cols;

	switch(WaveletIndex){
	case D97:
		TempPtr = LineA+AnalDelay97;
		Anal = SubbandAnal97;
		break;
	case W53:
		TempPtr = LineA+AnalDelay53;
		Anal = SubbandAnal53;
		break;
	case W137:
		TempPtr = LineA+AnalDelay137;
		Anal = SubbandAnal137;
		break;
	case V1018:
		TempPtr = LineA+AnalDelay1018;
		Anal = SubbandAnal1018;
		break;
  	case ADEL:
		TempPtr = LineA+AnalDelayAdel;
		Anal = SubbandAnalAdel;
		break;

	}

	for (j=0; j<nLevels; j++){
		/* read one row */
		ImgPtr = in;
		LinePtr = TempPtr;
			
		for (i=0; i<nc; i++){
			*LinePtr++ = *ImgPtr++;
		}
			
		/* transform the line - result in LineB */
		Anal(LineA, nc);

		/* overwrite the input */
		ImgPtr = in;
		LinePtr = LineB;

		for (i=0; i<nc; i++){
			*ImgPtr++ = *LinePtr++;
		}	
			
		/* halved the size - extra sample goes to low pass */
		nc = (nc+1)>>1;

	} /* next level */
}

/*--------------------------------------------------------------------------------------------*/
/* perform 2-D subband analysis - odd length LP FIR filter */
void Subband2DAnal(REAL *in, int rows, int cols, int nLevels, int WaveletIndex)
{
	int r, c;
	int i, j;
	int nr, nc;
	REAL *ImgPtr, *LinePtr, *TempPtr;
	void (*Anal)(REAL *, int);

	
	nr = rows;
	nc = cols;

	switch(WaveletIndex){
	case D97:
		TempPtr = LineA+AnalDelay97;
		Anal = SubbandAnal97;
		break;
	case W53:
		TempPtr = LineA+AnalDelay53;
		Anal = SubbandAnal53;
		break;
	case W137:
		TempPtr = LineA+AnalDelay137;
		Anal = SubbandAnal137;
		break;
	case V1018:
		TempPtr = LineA+AnalDelay1018;
		Anal = SubbandAnal1018;
		break;
	case ADEL:
		TempPtr = LineA+AnalDelayAdel;
		Anal = SubbandAnalAdel;
		break;
	}

	for (j=0; j<nLevels; j++){
		
		/* filter in the horizontal direction */
		for (r=0; r<nr; r++){
			
			/* read one row */
			ImgPtr = in+r*cols;
			LinePtr = TempPtr;
			
			for (i=0; i<nc; i++){
				*LinePtr++ = *ImgPtr++;
			}
			
			/* transform the line - result in LineB */
			Anal(LineA, nc);

			/* overwrite the input */
			ImgPtr = in+r*cols;
			LinePtr = LineB;

			for (i=0; i<nc; i++){
				*ImgPtr++ = *LinePtr++;
			}	
			
		}	/* next row */
		
		/* filter in the vertical direction */
		for (c=0; c<nc; c++){
			
			/* read one col */
			ImgPtr = in+c;
			LinePtr = TempPtr;
			
			for (i=0; i<nr; i++){
				*LinePtr++ = *ImgPtr;
				ImgPtr += cols;
			}
			
			/* transform the line - result in LineB */
			Anal(LineA, nr);

			/* overwrite the input */
			ImgPtr = in+c;
			LinePtr = LineB;
			
			for (i=0; i<nr; i++){
				*ImgPtr = *LinePtr++;
				ImgPtr+=cols;
			}				
			
		}	/* next col */

		/* halved the size - extra sample goes to low pass */
		nr = (nr+1)>>1;
		nc = (nc+1)>>1;

	} /* next level */
}

/*--------------------------------------------------------------------------------------------*/
/* perform 1-D subband synthesis - odd length LP FIR filter */
void Subband1DSynt(REAL *in, int cols, int nLevels, int WaveletIndex)
{
	int c;
	int i, j;
	int nc;
	REAL *ImgPtr, *LinePtr;
	void (*Synt)(REAL *, int);
	
	switch(WaveletIndex){
	case D97:
		Synt = SubbandSynt97;
		break;
	case W53:
		Synt = SubbandSynt53;
		break;
	case W137:
		Synt = SubbandSynt137;
		break;
	case V1018:
		Synt = SubbandSynt1018;
		break;
	case ADEL:
		Synt = SubbandSyntAdel;
		break;
	}

	for (j=nLevels-1; j>=0; j--){

		nc = cols;

		/* get the synthesized length */
		for (i=1; i<=j; i++){
			nc = (nc+1)>>1;
		}
		
		/* filter in the horizontal direction */
		Synt(in, nc);
		
	}	/* next level */
}




/*--------------------------------------------------------------------------------------------*/
/* perform 2-D subband synthesis - odd length LP FIR filter */
void Subband2DSynt(REAL *in, int rows, int cols, int nLevels, int WaveletIndex)
{
	int r,c;
	int i, j;
	int nr, nc;
	REAL *ImgPtr, *LinePtr;
	void (*Synt)(REAL *, int);
	
	switch(WaveletIndex){
	case D97:
		Synt = SubbandSynt97;
		break;
	case W53:
		Synt = SubbandSynt53;
		break;	
	case W137:
		Synt = SubbandSynt137;
		break;
	case V1018:
		Synt = SubbandSynt1018;
		break;
	case ADEL:
		Synt = SubbandSyntAdel;
		break;
	}
	
	for (j=nLevels-1; j>=0; j--){

		nr = rows;
		nc = cols;

		/* get the synthesized length */
		for (i=1; i<=j; i++){
			nr = (nr+1)>>1;
			nc = (nc+1)>>1;
		}
		
		/* filter in the vertical direction */
		for (c=0; c<nc; c++){

			ImgPtr = in + c;
			LinePtr = LineA;
			for (i=0; i<nr; i++){
				*LinePtr++ = *ImgPtr;
				ImgPtr += cols;
			}
			
			/* synthesis */
			Synt(LineA, nr);
			
			/* result in LineA */
			LinePtr = LineA;
			ImgPtr = in + c;
			for (i=0; i<nr; i++){
				*ImgPtr = *LinePtr++;
				ImgPtr+=cols;
			}
		}


		/* filter in the horizontal direction */
		for (r=0; r<nr; r++){
			Synt(in+r*cols, nc);
		}

	}	/* next level */
}