wavelet.c

该程序把数字图像处理与小波变换结合起来

			WaveletInvertStep (wavelet, temp_in, temp_out,
				hLowSize[nsteps]+hHighSize[nsteps], symExt);
			
			// Copy back to image
			copy (temp_out+wavelet->npad, output + (j*hsize), 
				hLowSize[nsteps]+hHighSize[nsteps]);
      }
   }
	
   free(hLowSize);
   free(hHighSize);
	free(vLowSize);
   free(vHighSize);
	
   free(temp_in);
   free(temp_out);
}


/*----------------------------------------------------------------------------*/	
/*----------------------------------------------------------------------------*/
// Do symmetric extension of data using prescribed symmetries
//   Original values are in output[npad] through output[npad+size-1]
//   New values will be placed in output[0] through output[npad] and in
//      output[npad+size] through output[2*npad+size-1] (note: end values may
//      not be filled in)
//   left_ext = 1 -> extension at left bdry is   ... 3 2 1 | 0 1 2 3 ...
//   left_ext = 2 -> extension at left bdry is ... 3 2 1 0 | 0 1 2 3 ...
//   right_ext = 1 or 2 has similar effects at the right boundary
//
//   symmetry = 1  -> extend symmetrically
//   symmetry = -1 -> extend antisymmetrically
void WaveletSymmetricExtension(WAVELET *wavelet, Real *output, int size, int leftExt, int
										 rightExt, int symmetry)
{
	int i;
	int first = wavelet->npad, last = wavelet->npad + size-1;
	int originalFirst, originalLast, originalSize, period;
	int nextend;
	
	if (symmetry == -1) {
		if (leftExt == 1)
			output[--first] = 0;
		if (rightExt == 1)
			output[++last] = 0;
	}
	
	originalFirst = first;
	originalLast = last;
	
	originalSize = originalLast-originalFirst+1;
	
	period = 2 * (last - first + 1) - (leftExt == 1) - (rightExt == 1);
	
	if (leftExt == 2)
		output[--first] = symmetry*output[originalFirst];
	if (rightExt == 2)
		output[++last] = symmetry*output[originalLast];
	
	// extend left end
	nextend = MIN (originalSize-2, first);
	for (i = 0; i < nextend; i++) {
		output[--first] = symmetry*output[originalFirst+1+i];
	}
	
	// should have full period now -- extend periodically
	while (first > 0) {
		first--;
		output[first] = output[first+period];
	}
	
	// extend right end
	nextend = MIN (originalSize-2, 2*wavelet->npad+size-1 - last);
	for (i = 0; i < nextend; i++) {
		output[++last] = symmetry*output[originalLast-1-i];
	}
	
	// should have full period now -- extend periodically
	while (last < 2*wavelet->npad+size-1) {
		last++;
		output[last] = output[last-period];
	}
}



/*----------------------------------------------------------------------------*/	
/*----------------------------------------------------------------------------*/
// Do periodic extension of data using prescribed symmetries
//   Original values are in output[npad] through output[npad+size-1]
//   New values will be placed in output[0] through output[npad] and in
//      output[npad+size] through output[2*npad+size-1] (note: end values may
//      not be filled in)
void WaveletPeriodicExtension(WAVELET *wavelet, Real *output, int size)
{
	int first = wavelet->npad, last = wavelet->npad + size-1;
	
	// extend left periodically
	while (first > 0) {
		first--;
		output[first] = output[first+size];
	}
	
	// extend right periodically
	while (last < 2*wavelet->npad+size-1) {
		last++;
		output[last] = output[last-size];
	}
}

/*----------------------------------------------------------------------------*/	
/*----------------------------------------------------------------------------*/	
void WaveletTransformStep (WAVELET *wavelet, Real *input, Real *output, 
									int size, int symExt)
{
	int i, j;
	int lowSize=(size+1)/2;
	int leftExt, rightExt;

	if (wavelet->analysisLow->size%2){
		// odd filter length
		leftExt=rightExt=1;
	}
	else{
		leftExt=rightExt=2;
	}

	if (symExt){
		WaveletSymmetricExtension(wavelet, input, size, leftExt, rightExt, 1);
	}
	else{
		WaveletPeriodicExtension(wavelet, input, size);
	}

	//                      coarse  detail
	// xxxxxxxxxxxxxxxx --> HHHHHHHHGGGGGGGG
	for (i = 0; i < lowSize; i++)  {
		output[wavelet->npad+i] = 0.0;
		for (j = 0; j < wavelet->analysisLow->size; j++)  {
			output [wavelet->npad+i] += 
				input[wavelet->npad + 2*i + wavelet->analysisLow->firstIndex + j] *
				wavelet->analysisLow->coeff[j];
		}
	}

	for (i = lowSize; i < size; i++)  {
		output[wavelet->npad+i] = 0.0;
		for (j = 0; j < wavelet->analysisHigh->size; j++)  {
			output [wavelet->npad+i] += 
				input[wavelet->npad + 2*(i-lowSize) + wavelet->analysisHigh->firstIndex + j] * 
				wavelet->analysisHigh->coeff[j];
		}
	}
}



/*----------------------------------------------------------------------------*/	
/*----------------------------------------------------------------------------*/	
void WaveletInvertStep (WAVELET *wavelet, Real *input, Real *output, 
								int size, int symExt)
{
	int i, j;
	int leftExt, rightExt, symmetry;
	Real *temp;
	int firstIndex, lastIndex;
	
	// amount of low and high pass -- if odd # of values, extra will be
   // low pass
   int lowSize = (size+1)/2, highSize = size/2;
	
   symmetry = 1;
	
   if (wavelet->analysisLow->size % 2 == 0) {
		// even length filter -- do (2, X) extension
		leftExt = 2;
   } else {
		// odd length filter -- do (1, X) extension
		leftExt = 1;
   }
	
	if (size % 2 == 0) {
		// even length signal -- do (X, 2) extension
		rightExt = 2;
   } else {
		// odd length signal -- do (X, 1) extension
		rightExt = 1;
   }
	
	temp = (Real *)calloc(2*wavelet->npad+lowSize, sizeof(Real));
   for (i = 0; i < lowSize; i++) {
		temp[wavelet->npad+i] = input[wavelet->npad+i];
   }
	
	if (symExt){
		WaveletSymmetricExtension (wavelet, temp, lowSize, leftExt, rightExt, symmetry);
	}
   else{
		WaveletPeriodicExtension (wavelet, temp, lowSize);
	}
	
	// coarse  detail
   // HHHHHHHHGGGGGGGG --> xxxxxxxxxxxxxxxx
   for (i = 0; i < 2*wavelet->npad+size; i++){
		output[i] = 0.0;
	}
	
	firstIndex = wavelet->synthesisLow->firstIndex;
   lastIndex = wavelet->synthesisLow->size - 1 + firstIndex;
	
   for (i = -lastIndex/2; i <= (size-1-firstIndex)/2; i++)  {
      for (j = 0; j < wavelet->synthesisLow->size; j++)  {
			output[wavelet->npad + 2*i + firstIndex + j] +=
				temp[wavelet->npad+i] * wavelet->synthesisLow->coeff[j];
      }
   }
	
	leftExt = 2;
	
   if (wavelet->analysisLow->size % 2 == 0) {
		// even length filters
		rightExt = (size % 2 == 0) ? 2 : 1;
		symmetry = -1;
   } 
	else {
		// odd length filters
		rightExt = (size % 2 == 0) ? 1 : 2;
		symmetry = 1;
   }
	
   for (i = 0; i < highSize; i++) {
		temp[wavelet->npad+i] = input[wavelet->npad+lowSize+i];
   }
	
   if (symExt){
		WaveletSymmetricExtension (wavelet, temp, highSize, leftExt, rightExt,
			symmetry);
	}
   else{
		WaveletPeriodicExtension (wavelet, temp, highSize);
	}
	
	
   firstIndex = wavelet->synthesisHigh->firstIndex;
   lastIndex = wavelet->synthesisHigh->size - 1 + firstIndex;
	
   for (i = -lastIndex/2; i <= (size-1-firstIndex)/2; i++)  {
      for (j = 0; j < wavelet->synthesisHigh->size; j++)  {
			output[wavelet->npad + 2*i + firstIndex + j] +=
				temp[wavelet->npad+i] * wavelet->synthesisHigh->coeff[j];
      }
   }
	
   free(temp);
}

// copy length elements from p1 to p2
/*----------------------------------------------------------------------------*/	
/*----------------------------------------------------------------------------*/	
void copy (const Real *p1, Real *p2, const int length){
	int temp = length; 
	while(temp--) *p2++ = *p1++;
}

/*----------------------------------------------------------------------------*/	
/*----------------------------------------------------------------------------*/	
void copy_p1_skip (const Real *p1, const int stride1, Real *p2, const int length){
	int temp = length; 
	while(temp--){
		*p2++ = *p1; 
		p1 += stride1;
	}
}

/*----------------------------------------------------------------------------*/	
/*----------------------------------------------------------------------------*/	
void copy_p2_skip (const Real *p1, Real *p2, const int stride2, const int length){
	int temp = length; 
	while(temp--){
		*p2 = *p1++; 
		p2 += stride2;
	}
}

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

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