idwt_aux.cpp

《Visual C++小波变换技术与工程实践》靳济芳编著的光盘程序。

    }
    /* perform inverse vertical SA-DWT */
    ret=iSADWT1dInt(InBuf, InMaskBuf, OutBuf, OutMaskBuf, height, Filter, 
		    DWT_VERTICAL,  ((i>=(width>>1) && ZeroHigh == DWT_ZERO_HIGH)?DWT_ALL_ZERO:ZeroHigh));
    if(ret!=DWT_OK) {
      free(InBuf);free(OutBuf);free(InMaskBuf);free(OutMaskBuf);
      return(ret);
    }
    /* put the transformed column back */
    for(a=OutBuf, e=OutCoeff+i, c= OutMaskBuf, d= OutMask+i;
	a<OutBuf+height; a++, c++, e+=Width, d+=Width) {
      /* Scale and Leave 3 extra bits for precision reason */
      *a = ROUNDDIV(*a <<3, Filter->Scale); 
      if(*a > MaxCoeff || *a < MinCoeff) {
	free(InBuf);free(OutBuf);free(InMaskBuf);free(OutMaskBuf);
	return(DWT_COEFF_OVERFLOW);
      }
      *e= (Int) *a;
      *d = *c;
    }
  }    
  /* then horizontal synthesis */
  for(i=0,k=0;i<height;i++,k+=Width) {
    /* get a line of coefficients */
    for(a=InBuf, e=OutCoeff+k;a<InBuf+width;a++,e++) {
      *a = (Int) *e;
    }
    /* get a line of mask */
    memcpy(InMaskBuf, OutMask+k, sizeof(UChar)*width);
    /* Perform horizontal inverse SA-DWT */
    ret=iSADWT1dInt(InBuf, InMaskBuf, OutBuf, OutMaskBuf, width, Filter, 
		   DWT_HORIZONTAL, ZeroHigh);
    if(ret!=DWT_OK) {
      free(InBuf);free(OutBuf);free(InMaskBuf);free(OutMaskBuf);
      return(ret);
    }
    /* put the transformed line back */
    for(a=OutBuf,e=OutCoeff+k;a<OutBuf+width;a++,e++) {
      /*Scale and leave the sqrt(2)*sqrt(2) factor in the scale */
      *a = ROUNDDIV(*a , Filter->Scale<<2);
      if(*a > MaxCoeff || *a < MinCoeff) {
	free(InBuf);free(OutBuf);free(InMaskBuf);free(OutMaskBuf);
	return(DWT_COEFF_OVERFLOW);
      }
      *e = (Int) *a;
    }
    memcpy(OutMask+k, OutMaskBuf, sizeof(UChar)*width);
  }


  free(InBuf);free(OutBuf);free(InMaskBuf);free(OutMaskBuf);
  return(DWT_OK);
}
#else
/* Function:   SynthesizeOneLevelDbl()
   Description: Synthesize one level of wavelet coefficients
   Input:
   Width  -- Width of Image; 
   Height -- Height of Image;
   level -- current  level
   Filter     -- Filter Used.
   MaxCoeff, MinCoeff -- bounds of the output data;
   ZeroHigh -- Highpass bands are all zeros 
   Input/Output:

   OutCoeff   -- Ouput wavelet coefficients
   OutMask    -- Output mask corresponding to wavelet coefficients

   Return: DWT_OK if successful.
*/
 Int VTCIDWT:: SynthesizeOneLevelDbl(double *OutCoeff, UChar *OutMask, Int Width,
				Int Height, Int level, FILTER *Filter, Int ZeroHigh)
{
  double *InBuf, *OutBuf ;
  UChar *InMaskBuf, *OutMaskBuf;
  Int width = Width>>(level-1);
  Int height = Height>>(level-1);
  Int MaxLength = (width > height)?width:height;
  Int i,k,ret;
  double *a;
  UChar *c,*d;
  double *e;

  /* double check filter type */
  if(Filter->DWT_Type != DWT_DBL_TYPE) return(DWT_INTERNAL_ERROR);
  /* allocate line buffers */
  InBuf = (double *) malloc(sizeof(double)*MaxLength);
  InMaskBuf = (UChar *) malloc(sizeof(UChar)*MaxLength);
  OutBuf = (double *) malloc(sizeof(double)*MaxLength);
  OutMaskBuf = (UChar *) malloc(sizeof(UChar)*MaxLength);
  if(InBuf==NULL || InMaskBuf ==NULL || OutBuf == NULL || OutMaskBuf==NULL) 
    return(DWT_MEMORY_FAILED);
    

  /* vertical synthesis first */
  for(i=0;i<width;i++) {
    /* get a column of coefficients and mask */
    for(a=InBuf, e=OutCoeff+i, c= InMaskBuf, d= OutMask+i;
	a<InBuf+height; a++, c++, e+=Width, d+=Width) {
      *a = *e;
      *c = *d;
    }
    /* perform vertical inverse SA-DWT */
    ret=iSADWT1dDbl(InBuf, InMaskBuf, OutBuf, OutMaskBuf, height, Filter, 
		    DWT_VERTICAL, ((i>=(width>>1) && ZeroHigh==DWT_ZERO_HIGH)?DWT_ALL_ZERO:ZeroHigh));
    if(ret!=DWT_OK) {
      free(InBuf);free(OutBuf);free(InMaskBuf);free(OutMaskBuf);
      return(ret);
    }
    /* put the transformed column back */
    for(a=OutBuf, e=OutCoeff+i, c= OutMaskBuf, d= OutMask+i;
	a<OutBuf+height; a++, c++, e+=Width, d+=Width) {
      *e = *a;
      *d = *c;
    }
  }

  /* horizontal decomposition first*/
  for(i=0,k=0;i<height;i++,k+=Width) {
    /* get a line of coefficients */
    for(a=InBuf, e=OutCoeff+k;a<InBuf+width;a++,e++) {
      *a =  *e;
    }
    /* get a line of mask */
    memcpy(InMaskBuf, OutMask+k, sizeof(UChar)*width);
    /* Perform horizontal inverse SA-DWT */
    ret=iSADWT1dDbl(InBuf, InMaskBuf, OutBuf, OutMaskBuf, width, Filter, 
		    DWT_HORIZONTAL, ZeroHigh);
    if(ret!=DWT_OK) {
      free(InBuf);free(OutBuf);free(InMaskBuf);free(OutMaskBuf);
      return(ret);
    }
    /* put the transformed line back */
    for(a=OutBuf,e=OutCoeff+k;a<OutBuf+width;a++,e++) {
      *e =  *a;
    }
    memcpy(OutMask+k, OutMaskBuf, sizeof(UChar)*width);
  }

  free(InBuf);free(OutBuf);free(InMaskBuf);free(OutMaskBuf);
  return(DWT_OK);
}
#endif

/* Function: iSADWT1dInt() or iSADWT1dDbl()
   Description:  inverse 1-d SA-DWT
   Input:

   InBuf -- Input 1d data buffer
   InMaskBuf -- Input 1d mask buffer
   Length -- length of the input data
   Filter -- filter used
   Direction -- vertical or horizontal synthesis (used for inversible 
   mask synthsis)
   ZeroHigh -- DWT_ZERO_HIGH indicates the Highpass bands are all zeros
               DWT_ALL_ZERO indicates all the coefficients are zeros
	       DWT_NONZERO_HIGH indicates not all Highpass bands are zeros
   Output:
   
   OutBuf -- Synthesized 1d Data
   OutMask -- Mask for the synthesized 1d Data

   Return: return DWT_OK if successful
  
*/
#ifdef _DWT_INT_
 Int VTCIDWT:: iSADWT1dInt(Int *InBuf, UChar *InMaskBuf, Int *OutBuf, 
	       UChar *OutMaskBuf, Int Length, FILTER *Filter, 
		Int Direction, Int ZeroHigh)
{

  switch(Filter->DWT_Class){
  case DWT_ODD_SYMMETRIC:
    return(iSADWT1dOddSymInt(InBuf, InMaskBuf, OutBuf, OutMaskBuf, 
			    Length, Filter, Direction, ZeroHigh));
  case DWT_EVEN_SYMMETRIC:
    return(iSADWT1dEvenSymInt(InBuf, InMaskBuf, OutBuf, OutMaskBuf, 
			     Length, Filter, Direction, ZeroHigh));
    /*  case ORTHOGONAL:
	return(iSADWT1dOrthInt(InBuf, InMaskBuf, OutBuf, OutMaskBuf, 
			  Length, Filter, Direction, ZeroHigh)); */
  default:
    return(DWT_FILTER_UNSUPPORTED);
  }
}
#else
 Int VTCIDWT:: iSADWT1dDbl(double *InBuf, UChar *InMaskBuf, double *OutBuf, 
		UChar *OutMaskBuf, Int Length, FILTER *Filter, 
		Int Direction, Int ZeroHigh)
{

  switch(Filter->DWT_Class){
  case DWT_ODD_SYMMETRIC:
    return(iSADWT1dOddSymDbl(InBuf, InMaskBuf, OutBuf, OutMaskBuf, 
			    Length, Filter, Direction, ZeroHigh));
  case DWT_EVEN_SYMMETRIC:
    return(iSADWT1dEvenSymDbl(InBuf, InMaskBuf, OutBuf, OutMaskBuf, 
			     Length, Filter, Direction, ZeroHigh));
    /*  case ORTHOGONAL:
	return(iSADWT1dOrthDbl(InBuf, InMaskBuf, OutBuf, OutMaskBuf, 
			  Length, Filter, Direction, ZeroHigh)); */
  default:
    return(DWT_FILTER_UNSUPPORTED);
  }
}
#endif


/* Function: iSADWT1dOddSymInt() or iSADWT1dOddSymDbl()
   Description: 1D  iSA-DWT using Odd Symmetric Filter
   Input:

   InBuf -- Input 1d data buffer
   InMaskBuf -- Input 1d mask buffer
   Length -- length of the input data
   Filter -- filter used
   Direction -- vertical or horizontal synthesis (used for inversible 
   mask synthesis)
   ZeroHigh -- DWT_ZERO_HIGH indicates the Highpass bands are all zeros
               ALL_ZERO indicates all the coefficients are zeros
	       DWT_NONZERO_HIGH indicates not all Highpass bands are zeros

   Output:
   
   OutBuf -- Synthesized 1d Data
   OutMask -- Synthesized 1d Mask

   Return: return DWT_OK if successful

*/



#ifdef _DWT_INT_
 Int VTCIDWT:: iSADWT1dOddSymInt(DWTDATA *InBuf, UChar *InMaskBuf, DWTDATA *OutBuf, 
			    UChar *OutMaskBuf, Int Length, FILTER *Filter, 
			    Int Direction, Int ZeroHigh)
#else
 Int VTCIDWT:: iSADWT1dOddSymDbl(DWTDATA *InBuf, UChar *InMaskBuf, DWTDATA *OutBuf, 
			    UChar *OutMaskBuf, Int Length, FILTER *Filter, 
			    Int Direction, Int ZeroHigh)
#endif

{
  Int i;
  Int SegLength = 0;
  Int odd;
  Int start, end;
  UChar *a, *b, *c;
  Int ret;
  /* static Int count=0; */
  /* double check filter class and type */
  if(Filter->DWT_Class != DWT_ODD_SYMMETRIC) return(DWT_INTERNAL_ERROR);
#ifdef _DWT_INT_
  if(Filter->DWT_Type != DWT_INT_TYPE) return(DWT_INTERNAL_ERROR);
#else
  if(Filter->DWT_Type != DWT_DBL_TYPE) return(DWT_INTERNAL_ERROR);
#endif
  /* double check if Length is even */
  if(Length & 1) return(DWT_INTERNAL_ERROR);
  /* synthesize mask first */
  for(a=OutMaskBuf, b = InMaskBuf, c= InMaskBuf+(Length>>1); a<OutMaskBuf+Length;b++,c++) {
    if(Direction == DWT_VERTICAL) {
      if(*c == DWT_OUT2) { 
	*a++=DWT_OUT0;
	*a++=DWT_IN;
      }
      else if(*c == DWT_OUT3) {
	*a++=DWT_OUT1;
	*a++=DWT_IN;
      }
      else {
	*a++=*b;
	*a++=*c;
      }
    }
    else {
      if(*c == DWT_OUT1) { 
	*a++=DWT_OUT0;
	*a++=DWT_IN;
      }
      else {
	*a++ = *b;
	*a++ = *c;
      }
    }
  }

  /* initial OutBuf to zeros */
  memset(OutBuf, (UChar)0, sizeof(DWTDATA)*Length);
  if(ZeroHigh == DWT_ALL_ZERO) return(DWT_OK);
  
  i = 0;
  a = OutMaskBuf;
  while(i<Length) {
    /* search for a segment */
    while(i<Length && (a[i])!=DWT_IN) i++;
    start = i;
    if(i >= Length) break;
    while(i<Length && (a[i])==DWT_IN) i++;
    end = i;
    SegLength = end-start;
    odd = start%2;
    if(SegLength==1) { /* special case for single poInt */
#ifdef _DWT_INT_
      ret=SynthesizeSegmentOddSymInt(OutBuf+start, InBuf+(start>>1), 
				     InBuf+(Length>>1)+(start>>1), odd, 
				     SegLength, Filter, ZeroHigh);
#else
      ret=SynthesizeSegmentOddSymDbl(OutBuf+start, InBuf+(start>>1), 
				     InBuf+(Length>>1)+(start>>1), odd, 
				     SegLength, Filter,ZeroHigh);
#endif
      if(ret!=DWT_OK) return(ret);
    }
    else {
#ifdef _DWT_INT_
      ret=SynthesizeSegmentOddSymInt(OutBuf+start, InBuf+((start+1)>>1), 
				     InBuf+(Length>>1)+(start>>1), odd, 
				     SegLength, Filter,ZeroHigh);
      
#else
      ret=SynthesizeSegmentOddSymDbl(OutBuf+start, InBuf+((start+1)>>1),