owavelet2.c

LastWave

/*..........................................................................*/
/*                                                                          */
/*      L a s t W a v e    P a c k a g e 'owtrans2d' 2.1                    */
/*                                                                          */
/*      Copyright (C) 1998-2002 Geoff Davis, Emmanuel Bacry, Jerome Fraleu. */
/*                                                                          */
/*      The original program was written in C++ by Geoff Davis.             */
/*      Then it has been translated in C and adapted to LastWave by         */
/*      J. Fraleu and E. Bacry.                                             */
/*                                                                          */
/*      If you are interested in the C++ code please go to                  */
/*          http://www.cs.dartmouth.edu/~gdavis                             */
/*                                                                          */
/*      emails : [e-mail]                                       */
/*               [e-mail]                               */
/*               [e-mail]                             */
/*                                                                          */
/*..........................................................................*/
/*                                                                          */
/*      This program is a free software, you can redistribute it and/or     */
/*      modify it under the terms of the GNU General Public License as      */
/*      published by the Free Software Foundation; either version 2 of the  */
/*      License, or (at your option) any later version                      */
/*                                                                          */
/*      This program is distributed in the hope that it will be useful,     */
/*      but WITHOUT ANY WARRANTY; without even the implied warranty of      */
/*      MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the       */
/*      GNU General Public License for more details.                        */
/*                                                                          */
/*      You should have received a copy of the GNU General Public License   */
/*      along with this program (in a file named COPYRIGHT);                */
/*      if not, write to the Free Software Foundation, Inc.,                */
/*      59 Temple Place, Suite 330, Boston, MA  02111-1307  USA             */
/*                                                                          */
/*..........................................................................*/


#include "lastwave.h"
#include "owtrans2d.h"



/* Misc functions for copying arrays */
static void Copy3(LWFLOAT *p1,LWFLOAT *p2,int length)
{
  memcpy(p2,p1,length*sizeof(LWFLOAT));
}
static void Copy4(LWFLOAT *p1,int stride1, LWFLOAT *p2,int length)
{
 int temp = length; while(temp--) {*p2++ = *p1; p1 += stride1;}
}
static void Copy41 (LWFLOAT *p1, LWFLOAT *p2,  int stride2,int length)
{
  int temp = length; while(temp--) {*p2 = *p1++; p2 += stride2;}
}

/*---------------------------------------------------------------------------
* 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
*/

static void SymmetricExtension (OWAVELET2 w,LWFLOAT *output, int size, int left_ext, 
				    int right_ext, int symmetry)
{
  int i;
  int npad = w->npad;
  int first = npad, last = npad + size-1;
int nextend;
  int originalFirst;
  int originalLast ;
  int originalSize ;
int period;

  if (symmetry == -1) {
    if (left_ext == 1)
      output[--first] = 0;
    if (right_ext == 1)
      output[++last] = 0;
  }
  originalFirst = first;
  originalLast = last;
  originalSize = originalLast-originalFirst+1;

  period = 2 * (last - first + 1) - (left_ext == 1) - (right_ext == 1);

  if (left_ext == 2)
    output[--first] = symmetry*output[originalFirst];
  if (right_ext == 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*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*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)
*/

static void PeriodicExtension (OWAVELET2 w, LWFLOAT *output, int size)
{
  int npad = w->npad;
  int first = npad, last = npad + size-1;

  /* extend left periodically */
  while (first > 0) {
    first--;
    output[first] = output[first+size];
  }

  /* extend right periodically */
  while (last < 2*npad+size-1) {
    last++;
    output[last] = output[last-size];
  }
}

/*---------------------------------------------------------------------------*/

static void MallatToLinear (OWTRANS2 wtrans, LWFLOAT *mallat)
{
  int i, j, k;
  
  int *lowHsize = IntAlloc(wtrans->noct);
  int *lowVsize = IntAlloc(wtrans->noct);
  
  lowHsize[wtrans->noct-1] = (wtrans->hsize+1)/2;
  lowVsize[wtrans->noct-1] = (wtrans->vsize+1)/2;
  

  for (i = wtrans->noct-2; i >= 0; i--) {
    lowHsize[i] = (lowHsize[i+1]+1)/2;
    lowVsize[i] = (lowVsize[i+1]+1)/2;
  }
  
  /* move transformed image (in Mallat order) into linear array structure 
   special case for LL subband */
  for (j = 0; j < wtrans->subimage[0]->nrow; j++)
    for (i = 0; i < wtrans->subimage[0]->ncol; i++)
      wtrans->subimage[0]->pixels[j*(wtrans->subimage[0]->ncol)+i] = mallat[j*wtrans->hsize+i];
 
  for (k = 0; k < wtrans->noct; k++) {
    for (j = 0; j < wtrans->subimage[k*3+1]->nrow; j++)
      for (i = 0; i < wtrans->subimage[k*3+1]->ncol; i++)
	wtrans->subimage[k*3+1]->pixels[j*wtrans->subimage[k*3+1]->ncol+i] = 
	  mallat[j*wtrans->hsize+(lowHsize[k]+i)];

    for (j = 0; j < wtrans->subimage[k*3+2]->nrow; j++)
      for (i = 0; i < wtrans->subimage[k*3+2]->ncol; i++) 
	wtrans->subimage[k*3+2]->pixels[j*wtrans->subimage[k*3+2]->ncol+i] = 
	  mallat[(lowVsize[k]+j)*wtrans->hsize+i];

    for (j = 0; j < wtrans->subimage[k*3+3]->nrow; j++)
      for (i = 0; i < wtrans->subimage[k*3+3]->ncol; i++)
	wtrans->subimage[k*3+3]->pixels[j*wtrans->subimage[k*3+3]->ncol+i] = mallat[(lowVsize[k]+j)*wtrans->hsize+(lowHsize[k]+i)];
  }

  Free(lowHsize);
  Free(lowVsize);
}

/*****************************************************************************/

static void LinearToMallat (OWTRANS2 wtrans, LWFLOAT *mallat)
{
  int i, j, k;

  int *lowHsize = IntAlloc(wtrans->noct);
  int *lowVsize = IntAlloc(wtrans->noct);
  
  lowHsize[wtrans->noct-1] = (wtrans->hsize+1)/2;
  lowVsize[wtrans->noct-1] = (wtrans->vsize+1)/2;
  
  for (i = wtrans->noct-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 < wtrans->subimage[0]->nrow; j++)
    for (i = 0; i < wtrans->subimage[0]->ncol; i++)
      mallat[j*wtrans->hsize+i] = wtrans->subimage[0]->pixels[j*wtrans->subimage[0]->ncol+i];

  for (k = 0; k < wtrans->noct; k++) {
    for (j = 0; j < wtrans->subimage[k*3+1]->nrow; j++)
      for (i = 0; i < wtrans->subimage[k*3+1]->ncol; i++)
	mallat[j*wtrans->hsize+(lowHsize[k]+i)] = 
	  wtrans->subimage[k*3+1]->pixels[j*wtrans->subimage[k*3+1]->ncol+i];

    for (j = 0; j < wtrans->subimage[k*3+2]->nrow; j++)
      for (i = 0; i < wtrans->subimage[k*3+2]->ncol; i++)
	mallat[(lowVsize[k]+j)*wtrans->hsize+i] = wtrans->subimage[k*3+2]->pixels[j*wtrans->subimage[k*3+2]->ncol+i];

    for (j = 0; j < wtrans->subimage[k*3+3]->nrow; j++)
      for (i = 0; i < wtrans->subimage[k*3+3]->ncol; i++)
	mallat[(lowVsize[k]+j)*wtrans->hsize+(lowHsize[k]+i)] = 
	  wtrans->subimage[k*3+3]->pixels[j*wtrans->subimage[k*3+3]->ncol+i];
  }
  
  Free(lowHsize);
  Free(lowVsize);
}



/***********************************************
 * 
 * Decomposition
 *
 ***********************************************/ 

static void TransformStep (OWAVELET2 w, LWFLOAT *input, LWFLOAT *output, int size, int sym_ext)
{
  int i, j;
  int npad = w->npad;
  OFILTER2 analysisLow = w->analysisLow;
  OFILTER2 analysisHigh = w->analysisHigh;

  int lowSize = (size+1)/2;
  int left_ext, right_ext;

  if (analysisLow->size %2) {
    /* odd filter length */
    left_ext = right_ext = 1;
  } else {
    left_ext = right_ext = 2;
  }

  if (sym_ext)
    SymmetricExtension (w,input, size, left_ext, right_ext, 1);
  else
    PeriodicExtension (w,input, size);
    
  /*                      coarse  detail
   xxxxxxxxxxxxxxxx --> HHHHHHHHGGGGGGGG */
  for (i = 0; i < lowSize; i++)  {
    output[npad+i] = 0.0;
    for (j = 0; j < analysisLow->size; j++)  {
      output [npad+i] += 
	input[npad + 2*i + analysisLow->firstIndex + j] *
	analysisLow->coeff[j];
    }
  }
  
  for (i = lowSize; i < size; i++)  {
    output[npad+i] = 0.0;
    for (j = 0; j < analysisHigh->size; j++)  {
      output [npad+i] += 
	input[npad + 2*(i-lowSize) + analysisHigh->firstIndex + j] * 
	analysisHigh->coeff[j];
    }
  }
}

/* 1d decomp */
static void OWtrans1dDecomp(OWAVELET2 w, LWFLOAT *input, LWFLOAT *output, int size, int noct, int sym_ext)
{
   int i;
   int currentIndex = 0;
   LWFLOAT *data[2];
   int lowSize = size, highSize;
   int npad = w->npad;
   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;

   /* data[0] and data[1] are padded with npad entries on each end */
   data [0] = FloatAlloc(2*npad+size);
   data [1] = FloatAlloc(2*npad+size);

   for (i = 0; i < size; i++)
     data[currentIndex][npad+i] = input[i];

   while (noct--)  {
     if (lowSize <= 2 && symmetric == 1) 
       Errorf("OWtrans1dDecomp() : Warning ! Reduce # of transform steps or increase signal size or switch to periodic extension Low pass subband is too small");

     /* Transform */
     Printf ("transforming, size = %d\n", lowSize);
     TransformStep (w, data[currentIndex], data[1-currentIndex],lowSize, sym_ext);
     
     highSize = lowSize/2;
     lowSize = (lowSize+1)/2;
     
     /* Copy high-pass data to output signal */ 
     Copy3 (data[1-currentIndex] + npad + lowSize, output + lowSize, highSize);

     for (i = 0; i < lowSize+highSize; i++) Printf ("%5.2f ", data[1-currentIndex][npad+i]);
     Printf ("\n\n");
     
     /* Now pass low-pass data (first 1/2 of signal) back to transform routine  */
     currentIndex = 1 - currentIndex;
   }
   
   /* Copy low-pass data to output signal */ 
   Copy3 (data[currentIndex] + npad, output, lowSize);
   
   Free(data[1]);
   Free(data[0]); /* ???*/
}

/* 2d decomp */
static void OWtrans2dDecomp (OWAVELET2 w, LWFLOAT *input, LWFLOAT *output, int hsize, int vsize,int noct, int sym_ext)
{
   int j;
   int hLowSize = hsize, hHighSize;
   int vLowSize = vsize, vHighSize;
   int npad = w->npad;
   LWFLOAT *temp_in;
   LWFLOAT *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;