}
	}
		
	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; i<WaveletTransform->nSubbands; i++){
      /* go through each subband */
      sMax=-MaxReal;
      sMin= MaxReal;

      /* find dynamic range */
      for (j=0; j<WaveletTransform->subbandHSize[i]*WaveletTransform->subbandVSize[i]; j++){
         tempValue=WaveletTransform->subbandPtr[i][j];
         if (tempValue>sMax){
            sMax=tempValue;
         }
         else if (tempValue<sMin){
            sMin=tempValue;
         }
      }

      sRange=sMax-sMin;

      /* scale */
      for (j=0; j<WaveletTransform->subbandHSize[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; k<WaveletTransform->nsteps;k++){
			/* draw the lines */
			for (j=0; j<vLength; j++){
				image->pixel[j][lowHSize]=0;
			}

			for (j=0; j<hLength; j++){
				image->pixel[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; i<WaveletTransform->nSubbands; i++){
      /* go through each subband */
		fprintf(in, "Subband #%d\n", i);
		fprintf(in, "---------------------------------------------------------------------------\n");

      for (j=0; j<WaveletTransform->subbandHSize[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 );
}