owavelet2.c

LastWave

   temp_in = FloatAlloc(2*npad+MAX(hsize,vsize));
   temp_out = FloatAlloc(2*npad+MAX(hsize,vsize));

   Copy3 (input, output, hsize*vsize);

   while (noct--)  {
      if ((hLowSize <= 2 || vLowSize <= 2) && sym_ext == 1) 
	    Errorf ("Reduce # of transform steps or increase signal size or switch to periodic extension Low pass subband is too small");

      /* Do a convolution on the low pass portion of each row */
      for (j = 0; j < vLowSize; j++)  {
	   /* Copy row j to data array */ 
	   Copy3 (output+(j*hsize), temp_in+npad, hLowSize);
	 
	   /* Convolve with low and high pass filters */
	   TransformStep (w,temp_in, temp_out, hLowSize, sym_ext);

	   /* Copy back to image */ 
	   Copy3 (temp_out+npad, output+(j*hsize), hLowSize);
      }

      /* Now do a convolution on the low pass portion of  each column */
      for (j = 0; j < hLowSize; j++)  {
	   /* Copy column j to data array */ 
	   Copy4 (output+j, hsize, temp_in+npad, vLowSize);
	 
	  /* Convolve with low and high pass filters */
	  TransformStep (w,temp_in, temp_out, vLowSize, sym_ext);

	 /* Copy back to image */ 
	 Copy41 (temp_out+npad, output+j, hsize, vLowSize);
      }

      /* Now convolve low-pass portion again */
      hHighSize = hLowSize/2;
      hLowSize = (hLowSize+1)/2;
      vHighSize = vLowSize/2;
      vLowSize = (vLowSize+1)/2;
   }

  Free(temp_out);
  Free(temp_in);
}


/* C-function to perform decomposition */
void OWt2d (OWTRANS2 wtrans, int noct)
{  
  LWFLOAT *temp;
  
  SetWaveletOWtrans2(wtrans,NULL);

  if (wtrans->wavelet == NULL) Errorf("OWt2d() : Weird error : Wavelet is not initialized");
  if (wtrans->original->ncol == 0 || wtrans->original->nrow == 0) Errorf("OWt2d() : No image to analyze");
   
  SetNOctOWtrans2(wtrans,noct,wtrans->original->ncol,wtrans->original->nrow);
         
  temp = FloatAlloc(wtrans->hsize*wtrans->vsize);
  TempPtr(temp);
  
   
  OWtrans2dDecomp (wtrans->wavelet,wtrans->original->pixels, temp, wtrans->hsize, wtrans->vsize, wtrans->noct,-1);

  /* linearize data */
  MallatToLinear (wtrans, temp);
}


/* The corresponding command */
 void C_OWt2d(char **argv)
{
  OWTRANS2 wtrans;
  int noct;
  
  argv = ParseArgv(argv,tOWTRANS2_,NULL,&wtrans,tINT,&noct,0);
  
  if (wtrans ==NULL) wtrans= GetOWtrans2Cur(); 
  
  OWt2d(wtrans,noct);
}



/***********************************************
 * 
 * Reconstruction
 *
 ***********************************************/ 


static void InvertStep (OWAVELET2 w, LWFLOAT *input, LWFLOAT *output, int size, int sym_ext)
{
   int i, j;
   int left_ext, right_ext, symmetry;
   /* amount of low and high pass -- if odd # of values, extra will be
      low pass */
   int lowSize = (size+1)/2, highSize = size/2;
  int npad = w->npad;
  OFILTER2 analysisLow = w->analysisLow;
  OFILTER2 analysisHigh = w->analysisHigh;
  OFILTER2 synthesisLow = w->synthesisLow;
  OFILTER2 synthesisHigh = w->synthesisHigh;  
   int firstIndex;
   int lastIndex;
   LWFLOAT *temp;

   symmetry = 1;
   if (analysisLow->size % 2 == 0) {
    /* even length filter -- do (2, X) extension */
     left_ext = 2;
   } else {
     /* odd length filter -- do (1, X) extension */
     left_ext = 1;
   }

   if (size % 2 == 0) {
     /* even length signal -- do (X, 2) extension */
     right_ext = 2;
   } else {
     /* odd length signal -- do (X, 1) extension */
     right_ext = 1;
   }

   temp = FloatAlloc(2*npad+lowSize);
   for (i = 0; i < lowSize; i++) {
     temp[npad+i] = input[npad+i];
   }

   if (sym_ext)
     SymmetricExtension (w,temp, lowSize, left_ext, right_ext, symmetry);
   else
     PeriodicExtension (w,temp, lowSize);

   /* coarse  detail
    HHHHHHHHGGGGGGGG --> xxxxxxxxxxxxxxxx */
   for (i = 0; i < 2*npad+size; i++)
     output[i] = 0.0;

   firstIndex = synthesisLow->firstIndex;
   lastIndex = synthesisLow->size - 1 + firstIndex;

   for (i = -lastIndex/2; i <= (size-1-firstIndex)/2; i++)  {
      for (j = 0; j < synthesisLow->size; j++)  {
	output[npad + 2*i + firstIndex + j] +=
	  temp[npad+i] * synthesisLow->coeff[j];
      }
   }

   left_ext = 2;

   if (analysisLow->size % 2 == 0) {
    /* even length filters */
     right_ext = (size % 2 == 0) ? 2 : 1;
     symmetry = -1;
   } else {
     /* odd length filters */
     right_ext = (size % 2 == 0) ? 1 : 2;
     symmetry = 1;
   }

   for (i = 0; i < highSize; i++) {
     temp[npad+i] = input[npad+lowSize+i];
   }
   if (sym_ext)
     SymmetricExtension (w,temp, highSize, left_ext, right_ext,
			  symmetry);
   else
     PeriodicExtension (w,temp, highSize);


   firstIndex = synthesisHigh->firstIndex;
   lastIndex = synthesisHigh->size - 1 + firstIndex;

   for (i = -lastIndex/2; i <= (size-1-firstIndex)/2; i++)  {
      for (j = 0; j < synthesisHigh->size; j++)  {
	output[npad + 2*i + firstIndex + j] +=
	  temp[npad+i] * synthesisHigh->coeff[j];
      }
   }

   Free(temp);
}

static void OWtrans1dRecons (OWAVELET2 w, LWFLOAT *input, LWFLOAT *output, int size, int noct, int sym_ext)
{
   int i;
   int currentIndex = 0;
   LWFLOAT *data[2];
   int npad = w->npad;
   int *lowSize;
   int *highSize;
   int symmetric = w->symmetric;
   
   /* If form of extension unspecified, default to symmetric
    extensions for symmetrical filters and periodic extensions for
    asymmetrical filters */
   if (sym_ext == -1) sym_ext = symmetric;

   lowSize = (int *)  Malloc(sizeof(int)*noct);
   highSize = (int *)  Malloc(sizeof(int)*noct);

   lowSize[0] = (size+1)/2;
   highSize[0] = size/2;

   for (i = 1; i < noct; i++) {
     lowSize[i] = (lowSize[i-1]+1)/2;
     highSize[i] = lowSize[i-1]/2;
   }

   data [0] = FloatAlloc(2*npad+size);
   data [1] = FloatAlloc(2*npad+size);

   Copy3 (input, data[currentIndex]+npad, lowSize[noct-1]); 

   while (noct--)  {
     
     /* grab the next high-pass component */ 
     Copy3 (input + lowSize[noct], 
	 data[currentIndex]+npad+lowSize[noct], highSize[noct]);
     
     /* Combine low-pass data (first 1/2^n of signal) with high-pass
      data (next 1/2^n of signal) to get higher resolution low-pass data */
     InvertStep (w,data[currentIndex], data[1-currentIndex], 
		  lowSize[noct]+highSize[noct], sym_ext);
     
     /* Now pass low-pass data (first 1/2 of signal) back to transform routine */
     currentIndex = 1 - currentIndex;
   }

   /* Copy inverted signal to output signal */ 
   Copy3 (data[currentIndex]+npad, output, size);

   Free(highSize);
   Free(lowSize); /* ???*/

   Free(data[1]);
   Free(data[0]); /* ???*/
}

static void OWtrans2dRecons (OWAVELET2 w, LWFLOAT *input, LWFLOAT *output, int hsize, int vsize, int noct, int sym_ext)
{
   int i, j;
   int npad = w->npad;
   int *hLowSize,*hHighSize;
   int *vLowSize ,*vHighSize;
    LWFLOAT *temp_in ,*temp_out;
   int symmetric = w->symmetric;

   /* If form of extension unspecified, default to symmetric
    extensions for symmetrical filters and periodic extensions for
    asymmetrical filters */
   if (sym_ext == -1)
     sym_ext = symmetric;

   hLowSize = (int *) Malloc(sizeof(int)*noct);
   hHighSize = (int *) Malloc(sizeof(int)*noct);
   vLowSize =(int *) Malloc(sizeof(int)*noct),
   vHighSize = (int *) Malloc(sizeof(int)*noct);

   hLowSize[0] = (hsize+1)/2;
   hHighSize[0] = hsize/2;
   vLowSize[0] = (vsize+1)/2;
   vHighSize[0] = vsize/2;

   for (i = 1; i < noct; i++) {
     hLowSize[i] = (hLowSize[i-1]+1)/2;
     hHighSize[i] = hLowSize[i-1]/2;
     vLowSize[i] = (vLowSize[i-1]+1)/2;
     vHighSize[i] = vLowSize[i-1]/2;
   }

   temp_in = FloatAlloc(2*npad+MAX(hsize,vsize));
   temp_out = FloatAlloc(2*npad+MAX(hsize,vsize));

   Copy3 (input, output, hsize*vsize); 

   while (noct--)  {
      /* Do a reconstruction for each of the columns */
      for (j = 0; j < hLowSize[noct]+hHighSize[noct]; j++)  {
	 /* Copy column j to data array */ 
	 Copy4 (output+j, hsize, temp_in+npad, 
	       vLowSize[noct]+vHighSize[noct]);
	 
	 /* Combine low-pass data (first 1/2^n of signal) with high-pass
	  data (next 1/2^n of signal) to get higher resolution low-pass data */
	 InvertStep (w,temp_in, temp_out,vLowSize[noct]+vHighSize[noct], sym_ext);

	 /* Copy back to image */ 
	 Copy41 (temp_out+npad, output+j, hsize,
	       vLowSize[noct]+vHighSize[noct]);
      }

      /* Now do a reconstruction pass for each row */
      for (j = 0; j < vLowSize[noct]+vHighSize[noct]; j++)  {
	 /* Copy row j to data array */ 
	 Copy3 (output + (j*hsize), temp_in+npad, 
	       hLowSize[noct]+hHighSize[noct]);

	 /* Combine low-pass data (first 1/2^n of signal) with high-pass
	  data (next 1/2^n of signal) to get higher resolution low-pass data */
	 InvertStep (w,temp_in, temp_out,
		      hLowSize[noct]+hHighSize[noct], sym_ext);
	 
	 /* Copy back to image */ 
	 Copy3 (temp_out+npad, output + (j*hsize), 
	       hLowSize[noct]+hHighSize[noct]);
      }
   }

  Free(hLowSize);
  Free(hHighSize);
  Free(vLowSize);
  Free(vHighSize);
  Free(temp_in);
  Free(temp_out);
}

/* The C function for reconstruction */
void OWt2r (OWTRANS2 wtrans, IMAGE image)
{
  LWFLOAT *temp;
    
  /* We check the fact that the owtrans2 should not be empty */
  CheckOWtrans2(wtrans);

  /* Some Allocation */
  temp = FloatAlloc(wtrans->hsize*wtrans->vsize);
  TempPtr(temp);

  /* put data in Mallat format */
  LinearToMallat (wtrans,temp);

  /* Allocation of reconstructed image */
  SizeImage(image,wtrans->hsize,wtrans->vsize);
  
  /* Let's reconstruct */
  OWtrans2dRecons(wtrans->wavelet,temp, image->pixels, wtrans->hsize, wtrans->vsize,wtrans->noct,-1);
}

/* The corresponding command */
void C_OWt2r(char **argv)
{
  OWTRANS2 wtrans;
  IMAGE image;
  
  argv = ParseArgv(argv,tOWTRANS2_,NULL,&wtrans,-1);
  if (wtrans==NULL) wtrans= GetOWtrans2Cur(); 
  
  argv = ParseArgv(argv,tIMAGE_,wtrans->original,&image,0);

  OWt2r(wtrans,image);
}