idwt_aux.cpp

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

				     InBuf+(Length>>1)+(start>>1), odd, 
				     SegLength, Filter,ZeroHigh);
#endif
      if(ret!=DWT_OK) return(ret);
    }
  }
  return(DWT_OK);
}

/* 
   Function:  SynthesizeSegmentOddSymInt() or  SynthesizeSegmentOddSymDbl()
   Description: SA-Synthesize a 1D segment based on its InitPos and Length using
                Odd-Symmetric Filters
   Input:

   InL -- Input lowpass data;
   InH -- Input highpass data;
   PosFlag -- Start position of this segment (ODD or EVEN);
   Length -- Length of this Segment;
   Filter -- Filter used;
   ZeroHigh -- Highpass band is zero or not.

   Output:

   Out -- Output synthesized data segment;

   Return: Return DWT_OK if Successful
*/
#ifdef _DWT_INT_
 Int VTCIDWT::  SynthesizeSegmentOddSymInt(DWTDATA *Out, DWTDATA *InL, DWTDATA *InH, 
				     Int PosFlag, Int Length, FILTER *Filter, Int ZeroHigh)
#else
 Int VTCIDWT::  SynthesizeSegmentOddSymDbl(DWTDATA *Out, DWTDATA *InL, DWTDATA *InH, 
				     Int PosFlag, Int Length, FILTER *Filter, Int ZeroHigh)
#endif

{


#ifdef _DWT_INT_
  Short *LPCoeff = (Short *)Filter->LPCoeff, *HPCoeff = (Short *)Filter->HPCoeff; 
  Short *fi;
#else
  double *LPCoeff =(double *) Filter->LPCoeff, *HPCoeff =(double *) Filter->HPCoeff; 
  double *fi; 
#endif
  Int ltaps = Filter->LPLength, htaps = Filter->HPLength; /* filter lengths*/
  Int loffset = ltaps/2, hoffset = htaps/2;  /* Filter offset */
  Int border = (ltaps>htaps)?ltaps:htaps; /*the larger of ltaps and htaps */
  Int m, n;
  DWTDATA *f,*pixel, *pixel1, val, *buf, *a, *tmp_out;
  Int r = Length-1;
  if(Length == 1) {
    PosFlag = 0;
	ZeroHigh = DWT_ZERO_HIGH; // hjlee 0928
  }
  /* allocate proper buffer */
  buf = (DWTDATA *) malloc((Length+2*border)*sizeof(DWTDATA));
  if(buf==NULL)  return(DWT_MEMORY_FAILED);
  
  for(m=0;m<Length;m++) Out[m]= (DWTDATA) 0;
  for(m=0;m<Length+2*border;m++) buf[m]= (DWTDATA) 0;
  /*  fill in the low pass coefficients by upsampling by 2*/
  a = buf+border;
  for(m=PosFlag; m< Length; m+=2) {
    a[m] = InL[m>>1];
  }
  /* symmetric extension */
  for (m=1 ; m<=border; m++) {
    a[-m] =   a[m];  /* to allow Shorter seg */
    a[r+m] =  a[r-m]; 
  }

  a = buf + border;
  f = buf + border + Length;
  tmp_out = Out;
  for (; a<f; a ++) {
    /* filter the pixel with lowpass filter */
    for( fi=LPCoeff, pixel=a-loffset, pixel1=pixel+ltaps-1, val=0, n=0; 
	 n<(ltaps >>1); 
	 n++, fi++, pixel++, pixel1--) {
      val += ( *fi * (*pixel + *pixel1));
    }
    val += ( *fi * *pixel);
    *tmp_out++ = val;
  }

  /* if Highpass are zeros, skip */
  if(ZeroHigh== DWT_NONZERO_HIGH) {
    for(m=0;m<Length+2*border;m++) buf[m]= (DWTDATA) 0;
    a = buf+border;
    for(m=1-PosFlag;m< Length; m+=2) {
      a[m] = InH[m>>1];
    }
    /* symmetric extension */
    for (m=1 ; m<=border; m++) {
      a[-m] =   a[m];  /* to allow Shorter seg */
      a[r+m] =  a[r-m]; 
    } 
    a = buf + border;
    tmp_out = Out;
    for (; a<f; a ++){
      /* filter the pixel with highpass filter */
      for(fi=HPCoeff, pixel=a-hoffset, pixel1=pixel+htaps-1, val=0, n=0; 
	  n<(htaps>>1); n++, fi++, pixel++, pixel1--) {
	val += ( *fi *  (*pixel + *pixel1));
      }
      val += ( *fi *  *pixel);
      *tmp_out++ += val;
    }
  }
  free(buf);
  return(DWT_OK);
}


#ifdef _DWT_INT_
 Int VTCIDWT:: iSADWT1dEvenSymInt(DWTDATA *InBuf, UChar *InMaskBuf, DWTDATA *OutBuf, 
			    UChar *OutMaskBuf, Int Length, FILTER *Filter, 
			    Int Direction, Int ZeroHigh)
#else
 Int VTCIDWT:: iSADWT1dEvenSymDbl(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;
  /* double check filter class and type */
  if(Filter->DWT_Class != DWT_EVEN_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=SynthesizeSegmentEvenSymInt(OutBuf+start, InBuf+(start>>1), 
				     InBuf+(Length>>1)+(start>>1), odd, 
				     SegLength, Filter, ZeroHigh);
#else
      ret=SynthesizeSegmentEvenSymDbl(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=SynthesizeSegmentEvenSymInt(OutBuf+start, InBuf+(start>>1), 
				     InBuf+(Length>>1)+((start+1)>>1), odd, 
				     SegLength, Filter,ZeroHigh);
      
#else
      ret=SynthesizeSegmentEvenSymDbl(OutBuf+start, InBuf+(start>>1), 
				     InBuf+(Length>>1)+((start+1)>>1), odd, 
				     SegLength, Filter,ZeroHigh);
#endif
      if(ret!=DWT_OK) return(ret);
    }
  }
  return(DWT_OK);
}

/* 
   Function:  SynthesizeSegmentEvenSymInt() or  SynthesizeSegmentEvenSymDbl()
   Description: SA-Synthesize a 1D segment based on its InitPos and Length using
                Even-Symmetric Filters
   Input:

   InL -- Input lowpass data;
   InH -- Input highpass data;
   PosFlag -- Start position of this segment (ODD or EVEN);
   Length -- Length of this Segment;
   Filter -- Filter used;
   ZeroHigh -- Highpass band is zero or not.

   Output:

   Out -- Output synthesized data segment;

   Return: Return DWT_OK if Successful
*/
#ifdef _DWT_INT_
 Int VTCIDWT:: SynthesizeSegmentEvenSymInt(DWTDATA *Out, DWTDATA *InL, DWTDATA *InH, 
				     Int PosFlag, Int Length, FILTER *Filter, Int ZeroHigh)
#else
 Int VTCIDWT:: SynthesizeSegmentEvenSymDbl(DWTDATA *Out, DWTDATA *InL, DWTDATA *InH, 
				     Int PosFlag, Int Length, FILTER *Filter, Int ZeroHigh)
#endif

{


#ifdef _DWT_INT_
  Short *LPCoeff = (Short *)Filter->LPCoeff, *HPCoeff = (Short *)Filter->HPCoeff; 
  Short *fi;
#else
  double *LPCoeff = (double *) Filter->LPCoeff, *HPCoeff = (double *)Filter->HPCoeff; 
  double *fi;
#endif
  Int ltaps = Filter->LPLength, htaps = Filter->HPLength; /* filter lengths*/
  Int loffset = ltaps/2, hoffset = htaps/2;  /* Filter offset */
  Int border = (ltaps>htaps)?ltaps:htaps; /*the larger of ltaps and htaps */
  Int m, n;
  DWTDATA *f,*pixel, *pixel1, val, *buf, *a, *tmp_out;
  Int r = Length-1;
  
  if(Length == 1) {
   /*  *Out = *InL  * Filter->Scale; */
    PosFlag = 0;
	ZeroHigh = DWT_ZERO_HIGH; // hjlee 0928
    /* return(DWT_OK); */
  }
  /* allocate proper buffer */
  buf = (DWTDATA *) malloc((Length+2*border+1)*sizeof(DWTDATA));
  if(buf==NULL)  return(DWT_MEMORY_FAILED);
  
  for(m=0;m<Length;m++) Out[m]= (DWTDATA) 0;
  for(m=0;m<Length+2*border+1;m++) buf[m]= (DWTDATA) 0;
  /*  fill in the low pass coefficients by upsampling by 2*/
  a = buf+border+1;
  
  for(m=-PosFlag; m< Length; m+=2) {
    a[m] = InL[(m+1)>>1];
  }
  
  /* symmetric extension */
  for (m=1 ; m<=border; m++) {
    a[-m-1] =   a[m-1];  /* to allow Shorter seg */ /*symmetric poInt a[-1] */
    a[r+m] =  a[r-m]; /*symmetric poInt a[r] */
  }

  f = buf + border + 1 + Length;
  tmp_out = Out;
  for (; a<f; a ++) {
    /* filter the pixel with lowpass filter */
    for( fi=LPCoeff, pixel=a-loffset, pixel1=pixel+ltaps-1, val=0, n=0; 
	 n<(ltaps>>1); n++, fi++, pixel++, pixel1--)
      val += ( *fi * (*pixel + *pixel1));
    *tmp_out++ = val;
  }
  /* if Highpass are zeros, skip */
  if(ZeroHigh== DWT_NONZERO_HIGH) {
    for(m=0;m<Length+2*border+1;m++) buf[m]= (DWTDATA) 0;
    a = buf+border+1;
    for(m=PosFlag;m< Length; m+=2) {
      a[m] = InH[m>>1];
    }
    /* anti-symmetric extension */
    for (m=1 ; m<=border; m++) {
      a[-m-1] =   -a[m-1];  /* to allow Shorter seg */ /*symmetric poInt a[-1] */
      a[r+m] =  -a[r-m];  /*symmetric poInt a[r] */
    } 
    tmp_out = Out;
    for (; a<f; a ++){
      /* filter the pixel with highpass filter */
      for(fi=HPCoeff, pixel=a-hoffset, pixel1 = pixel+ htaps-1, val=0, n=0; 
	  n<(htaps>>1); n++, fi++, pixel++, pixel1--) {
	val += ( *fi *  (*pixel - *pixel1));
      }
      *tmp_out++ += val;
    }
  }
  free(buf);
  return(DWT_OK);
}

#ifdef _DWT_INT_
#undef _DWT_INT_
#define _DWT_DBL_
#include "idwt_aux.cpp"
#endif