// Please see the README file or the LEGAL NOTES-section in the manual before modifying, compiling or using 'xmrm'
//   Idea: Manfred Kopp
//   Programming: Gerhard Waldhör, Andreas Artmann

#include <stdlib.h>
#include <stdio.h>
#include <math.h>

#define coef_sm 24
#define coef_dt 24
#define nr_wave 11
#define XX 0,0,0,0,0
#define MIN_DECOMP 1
#define SQR2 1.41421356237309504880168872420969808

/* Image Data: */
/*
const wavelet_nr=10;
const nr_pixels=24;
float data_image[nr_pixels]={1,10,3,5, 99,13,8,4,2,1, 3,4,1,99, 0,8,165,2,3, 6,7,8,9,100 };

float image[nr_pixels],buffer[nr_pixels];

*/
/* WAVELETS: */

//---------------------------------------------------------  
// TRANSFORMATION FILTERS:
//---------------------------------------------------------  
float H[nr_wave][coef_sm]={
/* H: Spline 1-1 (Haar) */
{ 1.0/2.0, 1.0/2.0,  0,0,0,0,XX,XX, 0,0,0,XX },

/* H: Spline 1-3 */
{ 1.0/2.0, 1.0/2.0,  0,0,0,0,XX,XX, 0,0,0,XX  },

/* H: Spline 1-5 */
{ 1.0/2.0, 1.0/2.0,  0,0,0,0,XX,XX, 0,0,0,XX  }, 

/* H: Spline 3-3 */
{ 3.0/64.0, -9.0/64.0, -7.0/64.0, 45.0/64.0, 45.0/64.0, -7.0/64.0, -9.0/64.0, 3.0/64.0,  0,0,0,XX,  0,0,0,XX },

/* H: Spline 3-7 */
{ 
  35.0/16384.0,
 -105.0/16384.0,
 -195.0/16384.0,
  865.0/16384.0,
  363.0/16384.0,
 -3489.0/16384.0,
 -307.0/16384.0,
  11025.0/16384.0,

  11025.0/16384.0,
 -307.0/16384.0,
 -3489.0/16384.0,
  363.0/16384.0,
  865.0/16384.0,
 -195.0/16384.0,
 -105.0/16384.0,
  35.0/16384.0,
   0,0,0,XX 
},

/* H: Spline 2-4 */
{ 
  0.0,
  3.0/128.0,
  -3.0/64.0,
  -8.0/64.0,
  19.0/64.0,
  45.0/64.0,
  19.0/64.0,
  -8.0/64.0,
  -3.0/64.0,
  3.0/128.0,
  0,0,0,0,XX,XX 
},
  
/* H: Spline 2-6 */
{
  0.0,
  -5.0/1024.0,
  5.0/512.0,
  17.0/512.0,
  -39.0/512.0,
  -123.0/1024.0,
  81.0/256.0,
  175.0/256.0,
  81.0/256.0,  
  -123.0/1024.0,
  -39.0/512.0,  
  17.0/512.0,  
  5.0/512.0,
  -5.0/1024.0,
  XX,XX
},    
  
/* H: Pseudociflet_4 */
{
  0.0,
 -1.0 / 512.0,
  0.0,
  18.0 / 512.0,
 -16.0 / 512.0,
 -63.0 / 512.0,
  144.0 / 512.0,
  348.0 / 512.0,
  144.0 / 512.0,
 -63.0 / 512.0,
 -16.0 / 512.0,
  18.0 / 512.0,
  0.0,
 -1.0 / 512.0,

  0,0,0,XX 
},

/* H: Battle Lemarie */
{
  0.0,
 -0.002,
 -0.003,
  0.006,
  0.006,
 -0.013,
 -0.012,
  0.030,  /*  5 and 6 sign change from Mallat's paper */
  0.023,
 -0.078,
 -0.035,
  0.307,
  0.542,
  0.307,
 -0.035,
 -0.078,
  0.023,
  0.030,
 -0.012,
 -0.013,
  0.006,
  0.006,
 -0.003,
 -0.002
},

/* H: Spline 3-3inv */
{  1.0/8.0, 3.0/8.0,  3.0/8.0, 1.0/8.0, 0,XX,  0,0,0,0,XX,XX  },

/* H: Spline 3-1inv */
{  1.0/8.0, 3.0/8.0,  3.0/8.0, 1.0/8.0, 0,XX,  0,0,0,0,XX,XX  }
};

//---------------------------------------------------------  

float G[nr_wave][coef_dt]={
/* G: Spline 1-1 (Haar) */
{ 1.0/2.0, -1.0/2.0,  0,0,0,XX, 0,0,0,0,XX,XX },

/* G: Spline 1-3 */
{ -1.0/16.0, -1.0/16.0, 1.0/2.0, -1.0/2.0, 1.0/16.0, 1.0/16.0,  0,0,0,0,  0,0,0,0,XX,XX  },

/* G: Spline 1-5 */
{ 3.0/256.0, 3.0/256.0, -11.0/128.0, -11.0/128.0, 0.5, -0.5, 11.0/128.0, 11.0/128.0,-3.0/256.0,-3.0/256.0,
  0,0,0,0,XX,XX  },

/* G: Spline 3-3 */
{ -1.0/8.0, 3.0/8.0, -3.0/8.0, 1.0/8.0, 0,XX,  0,0,0,0,XX,XX  },

/* G: Spline 3-7 */
{ -1.0/8.0, 3.0/8.0, -3.0/8.0, 1.0/8.0, 0,XX,  0,0,0,0,XX,XX  },
  
/* G: Spline 2-4 */
{ -1.0/4.0, 2.0/4.0, -1.0/4.0, 0.0,  XX,XX,XX,XX  },

/* G: Spline 2-6 */
{ -1.0/4.0, 2.0/4.0, -1.0/4.0, 0.0,  XX,XX,XX,XX  },

/* G: Pseudocoiflet_4 */
{
  0.0,
  1.0 / 32.0,
  0.0,
 -9.0 / 32.0,
 16.0 / 32.0,
 -9.0 / 32.0,
  0.0,
  1.0 / 32.0,  

  0,0,
  0,0,0,0,XX,XX   
},

/* G: Battle Lemarie */
{
  0.0,
  0.002,
 -0.003,
 -0.006,
  0.006,
  0.013,
 -0.012,
 -0.030,  /*  5 and 6 sign change from Mallat's paper */
  0.023,
  0.078,
 -0.035,
 -0.307,
  0.542,
 -0.307,
 -0.035,
  0.078,
  0.023,
 -0.030,
 -0.012,
  0.013,
  0.006,
 -0.006,
 -0.003,
  0.002
},

/* G: Spline 3-3inv */
{-3.0/64.0, -9.0/64.0,  7.0/64.0, 45.0/64.0,-45.0/64.0, -7.0/64.0,  9.0/64.0, 3.0/64.0,  0,0,0,XX,  0,0,0,XX },

/* G: Spline 3-1inv */
{ 1.0/4.0, 3.0/4.0, -3.0/4.0, -1.0/4.0, XX,XX,XX,XX }
};    

//---------------------------------------------------------  
// INVERSE TRANSFORMATION FILTERS:
//---------------------------------------------------------  
float Hinv[nr_wave][2][coef_dt/2];
float Ginv[nr_wave][2][coef_sm/2];

//---------------------------------------------------------  
// FILTER SIZE:
//---------------------------------------------------------  
const int h_size[nr_wave]={
/* h_size: Spline 1-1 (Haar) */
2,
/* h_size: Spline 1-3 */
2,
/* h_size: Spline 1-5 */
2,
/* h_size: Spline 3-1 */
8,
/* h_size: Spline 3-7 */
16,
/* h_size: Spline 2-4 */
10,
/* h_size: Spline 2-6 */
14,
/* h_size: Pseudocoiflet_4 */
14,
/* h:size; Battle Lemarie */
24,
/* h_size: Spline 3-3inv */
4,
/* h_size: Spline 3-1inv */
4

};

const int g_size[nr_wave]={
/* g_size: Spline 1-1 (Haar) */
2,
/* g_size: Spline 1-3 */
6,
/* g_size: Spline 1-5 */
10,
/* g_size: Spline 3-3 */
4,
/* g_size: Spline 3-7 */
4,
/* g_size: Spline 2-4 */
4,
/* g_size: Spline 2-6 */
4,
/* g_size: Pseudocoiflet_4 */
8,
/* g:size; Battle Lemarie */
24,
/* g_size: Spline 3-3inv */
8,
/* g_size: Spline 3-1inv */
4
};

//---------------------------------------------------------  
// FILTER SYMETRY
//---------------------------------------------------------  
const int wave_sym[nr_wave]={
/* Spline 1-1 (Haar) */
0,
/* Spline 1-3 */
0,
/* Spline 1-5 */
0,
/* Spline 3-1 */
0,
/* Spline 3-7 */
0,
/* Spline 2-4 */
1,
/* Spline 2-6 */
1,
/* Pseudocoiflet_4 */
1,
/* Battle Lemarie */
1,
/* Spline 3-3inv */
0,
/* Spline 3-1inv */
0
};

//---------------------------------------------------------  
// COEF-POSITIONS: h_pos,g_pos | hinv_g_pos,hinv_u_pos | ginv_g_pos,ginv_u_pos
//---------------------------------------------------------  
/*const*/ int start_pos[nr_wave][6]={
/* 0: Coef-Positions Spline 1-1 (Haar) */
{ 0,0, 0,0, 0,0 },

/* 1: Coef-Positions Spline 1-3 */
{ 0,-2, -1,-1, 0,0 },

/* 2: Coef-Positions Spline 1-5 */
{ 0,-4, -2,-2, 0,0 },

/* 3: Coef-Positions Spline 3-3 */  
{ -3,-1, -1,-0, -2,-1 },

/* 4: Coef-Positions Spline 3-7 */  
{ -7,-1, -1,0, -4,-3 },

/* 5: Coef-Positions Spline 2-4 */  
{ -4,-1, -1,-0, -2,-1 },

/* 6: Coef-Positions Spline 2-6 */  
{ -6,-1, -1,-0, -3,-2 },

/* 7: Coef-Positions Pseudocoiflet_4 */
{ -6,-4,-2,-2,-3,-2 },

/* 8: Coef-Positione Battle Lemarie */
{ -11,-12,-6,-6,-5,-5 },

/* 9: Coef-Positions Spline 3-3inv */  
{ -1,-3, -2,-1, -1, 0 },

/* 10: Coef-Positions Spline 2-4inv */  
{ -1,-1, -1,0, -1,0 }

};
  
//---------------------------------------------------------  
// INDEX-MIRROR: 
//
// Examples: min=3, max=6: ...012|345|678... -> ...45543|345|54334...
//           min=2, max=4:    ...01|23|45... -> ...2332|23|3223...
//---------------------------------------------------------  
float mirror(int k,int min,int max, int testend, int mi, int alter_sign, int long adr,int i_mult,float *data)
{
  int m,d,p,l,newmax;
  float sign;

  sign=1.0; 
  if ((k<min) || (k>=max))
  {
    newmax=max-min;
    
    // Calculate mirror position:
    if (k<min) 
      p=-(k-min)-(1-mi); 
    else 
    { 
      p=k-min; 
      sign=-sign;
      if (testend) mi=(mi+1)%2;
    }

    if (newmax>mi)
    {
    
      l=p%(newmax-mi); // Scale position area to data field
      d=p/(newmax-mi); // If (d%2 !=0, then mirroring:
    
      if (d%2 != 0) 
        m=(newmax-1)-l+min;
      else 
      {  m=l+min; sign=-sign; }
    }
    else m=min;
  }
  else m=k;

  if (alter_sign && !testend) 
    return data[adr+i_mult*m]*sign;
  else
    return data[adr+i_mult*m];
}

//---------------------------------------------------------  
// WAVELET-INIT
//---------------------------------------------------------  
void Wave_Init(int print_inverse_filters)
{
  int i,j,k,pot,g_u;
  
//  j=wavelet_nr;
  for (j=0; j<nr_wave; j++)  
  {
    for (i=0; i<h_size[j]; i++) H[j][i]=H[j][i]*SQR2;
    for (i=0; i<g_size[j]; i++) G[j][i]=G[j][i]*SQR2;

    for (i=(coef_dt)/2-1; i>=(g_size[j]/2); i--)
    {
      Hinv[j][0][i]=0.0;
      Hinv[j][1][i]=0.0;
    }
    for (i=(coef_sm)/2-1; i>=(h_size[j]/2); i--)
    {
      Ginv[j][0][i]=0.0;
      Ginv[j][1][i]=0.0;      
    }

    if (!wave_sym[j])
    {
      pot=-start_pos[j][1]+1; //-g_pos+1
      g_u=pot%2; pot=-((g_u*2)-1);
      k=g_size[j]-1;
      for (i=0; i<(g_size[j]/2); i++)
      {
        Hinv[j][g_u][i]=G[j][k]*pot;
        pot=-pot; k--;
        Hinv[j][1-g_u][i]=G[j][k]*pot;
        pot=-pot; k--;
      }
  
      pot=-start_pos[j][0]+1; //-h_pos+1
      g_u=pot%2; pot=-((g_u*2)-1);
      k=h_size[j]-1;
      for (i=0; i<(h_size[j]/2); i++)
      {
        Ginv[j][g_u][i]=H[j][k]*pot;
        pot=-pot; k--;
        Ginv[j][1-g_u][i]=H[j][k]*pot;
        pot=-pot; k--;
      }
    }
    else
    {
      pot=-start_pos[j][1]; //-g_pos
      g_u=pot%2; pot=-((g_u*2)-1);
      k=0;
      for (i=0; i<(g_size[j]/2); i++)
      {
        Hinv[j][1-g_u][i]=G[j][k]*pot;
        pot=-pot; k++;
        Hinv[j][g_u][i]=G[j][k]*pot;
        pot=-pot; k++;
      }
  
      pot=-start_pos[j][0]; //-h_pos
      g_u=pot%2; pot=-((g_u*2)-1);
      k=h_size[j]-1;
      for (i=0; i<(h_size[j]/2); i++)
      {
        Ginv[j][g_u][i]=H[j][k]*pot;
        pot=-pot; k--;
        Ginv[j][1-g_u][i]=H[j][k]*pot;
        pot=-pot; k--;
      }
    }

    if (print_inverse_filters)
    {
      printf("Wavelet #%d inverse Filters:\n",j);
      printf("Hinv_g: ");
      for (i=0; i<g_size[j]/2; i++) printf("%f ",Hinv[j][0][i]);
      printf("\nHinv_u: ");
      for (i=0; i<g_size[j]/2; i++) printf("%f ",Hinv[j][1][i]);
      printf("\nGinv_g: ");
      for (i=0; i<h_size[j]/2; i++) printf("%f ",Ginv[j][0][i]);
      printf("\nGinv_u: ");
      for (i=0; i<h_size[j]/2; i++) printf("%f ",Ginv[j][1][i]);
      printf("\n");
    }
  }
  if (print_inverse_filters) printf("\n");
}



//---------------------------------------------------------  
// WAVELET-CONSTRUCTION-STEP
// in:  w_nr, nr_data, *data
// out: *wave
//---------------------------------------------------------  
void Construction(int w_nr, int nr_data, long adr, int i_mult, float *data, float *wave)
{
  int j,k,u,p,sym;
  long i;
  double sum;
  
  sym=wave_sym[w_nr];
  u=nr_data%2;

  // Smooth decomposition:  
  p=0;
  for (i=0; i<(nr_data>>1); i++)
  {
    sum=0;
    k=p+start_pos[w_nr][0]; // h_pos
    for (j=0; j<h_size[w_nr]; j++)
    {
      sum+=mirror(k,0,nr_data-u, 0, sym , 0, adr,i_mult,data)*H[w_nr][j]; 
//      printf("sum: %d %f\n",k,sum);
      k++;
    }
    wave[i]=sum;
    p+=2;
  }
  if (u) wave[(nr_data>>1)]=data[adr+(nr_data-1)*i_mult];


  // Detail decomposition:    
  p=0;
  for (i=(nr_data>>1)+u; i<nr_data; i++)
  {
    sum=0;
    k=p+start_pos[w_nr][1]; // g_pos
    for (j=0; j<g_size[w_nr]; j++)
    {
      sum+=mirror(k,0,nr_data-u, 0, sym, 0, adr,i_mult,data)*G[w_nr][j]; 
//      printf("sum: %d %f\n",k,sum);      
      k++;
    }
    wave[i]=sum;
    p+=2;
  }
   
  // Copy back to data:
  for (i=0; i<nr_data; i++) data[adr+i*i_mult]=wave[i];
}

//---------------------------------------------------------  
// WAVELET-RECONSTRUCTION-STEP
// in:  w_nr, nr_data, *wave
// out: *data
//---------------------------------------------------------  
void Reconstruction(int w_nr, int nr_data, long adr, int i_mult, float *wave, float *data)
{
  int j,k,l,h,p,pot,u,sym,sym_inv;
  long i;
  double sum;

  u=nr_data%2;
  h=(nr_data>>1)+u;
  sym=wave_sym[w_nr]; sym_inv=(sym+1)%2;

  p=0;
  for(i=0; i<nr_data-u; i+=2)
  {
    sum=0;
    l=start_pos[w_nr][2]+p; //hinv_g_pos
    for (j=0; j<g_size[w_nr]/2; j++)
    {
      k=j+l;    
      sum+=mirror(k,0,h-u, sym, 0, 0, adr,i_mult,wave)*Hinv[w_nr][0][j];
//      printf("sum_s: %d %f\n",k,sum);
    }

    l=start_pos[w_nr][4]+p+h; //g_inv_g_pos
    for (j=0; j<h_size[w_nr]/2; j++)
    {
      k=j+l;
      sum+=mirror(k,h,nr_data, sym, sym, sym_inv, adr,i_mult,wave)*Ginv[w_nr][0][j]; 
//      printf("sum_d: %d %f\n",k,sum);
    }
       
    data[i]=sum;
    p++;
  }
  if (u) data[nr_data-1]=wave[adr+(nr_data>>1)*i_mult];


  p=0;
  for(i=1; i<nr_data; i+=2)
  {
    sum=0;
    l=start_pos[w_nr][3]+p; //hinv_u_pos
    for (j=0; j<g_size[w_nr]/2; j++)
    {
      k=j+l;    
      sum+=mirror(k,0,h-u, sym, 0, 0, adr,i_mult,wave)*Hinv[w_nr][1][j]; 
    }

    l=start_pos[w_nr][5]+p+h; //g_inv_u_pos
    for (j=0; j<h_size[w_nr]/2; j++)
    {
      k=j+l;
      sum+=mirror(k,h,nr_data, sym, sym, sym_inv, adr,i_mult,wave)*Ginv[w_nr][1][j];
    }

    data[i]=sum;
    p++;
  }

  // Copy back to wave:
  for (i=0; i<nr_data; i++) wave[adr+i*i_mult]=data[i];
}

//---------------------------------------------------------  
// WAVELET-DECOMPOSITION
// in:  wavelet, size_x,size_y, *c
// out: *c
//---------------------------------------------------------  
void Decompose(int wavelet, int size_x, int size_y, int *wave_w, int *wave_h, float *c)
{
  int h,w,i,j,row,col,ug,l;
  unsigned long adr;
  float *buffer;
  
  if (size_x>size_y) i=size_x; else i=size_y;
  buffer=(float*)malloc(i*sizeof(float));
    
  w=size_x; h=size_y; l=0;
  wave_w[l]=w; wave_h[l]=h; l++;      
 
  while ( (w>MIN_DECOMP) && (h>MIN_DECOMP) )
  {
    for (row=0; row<h; row++)
    {
      adr=row*size_x;
      ug = w & 1;
      
      Construction(wavelet, w, adr, 1, c, buffer);
    }
    w = (w>>1)+ug;


    for (col=0; col<w; col++)
    {
      adr = col;
      ug = h & 1;

      Construction(wavelet, h, adr, size_x, c, buffer);
    }
    h = (h>>1)+ug;

//    printf("decompose: %d %d\n",w,h);
    wave_w[l]=w; wave_h[l]=h; l++;      
  }
  
  free(buffer);
}

//---------------------------------------------------------  
// WAVELET-COMPOSITION
// in:  wavelet, lvls, size_x,size_y, *c
// out: *c
//---------------------------------------------------------  
void Compose(int wavelet, int lvls, int size_x, int size_y, int *wave_w, int *wave_h, float *c)
{
  int h,w,i,row,col,half,ug,level;
  unsigned long adr;
  float *buffer;
  
  if (size_x>size_y) i=size_x; else i=size_y;
  buffer=(float*)malloc(i*sizeof(float));
  
  level=lvls;
 
  w=wave_w[level]; level--; h=wave_h[level];  
  while (level>=0)
//  while ((w<=size_x) && (h<=size_y))
  {
    for (col=0; col<w; col++)
    {
      adr = col;
      half = h>>1;
      ug=h & 1;
      
      Reconstruction(wavelet, h, adr, size_x, c, buffer);
    }
    w=wave_w[level];
//    printf("compose: %d %d\n",w,h); //****

    for (row=0; row<h; row++)
    { 
      adr=row*size_x;
      half = w>>1;
      ug=w & 1;

      Reconstruction(wavelet, w, adr, 1, c, buffer);      
    }

    w=wave_w[level]; level--; h=wave_h[level];
  }
  free(buffer);
}

//*********************************************************
/*
#define MIN_PARAM -5
#define MAX_PARAM 0
#define AVOID_MIRROR 4

void SearchParam()
{
  int i[6],j,ok; // please set first two values of start_pos

  // permutate start_pos:

  for (i[2]=MIN_PARAM; i[2]<=MAX_PARAM; i[2]++)
  {
    start_pos[wavelet_nr][2]=i[2];
    for (i[3]=MIN_PARAM; i[3]<=MAX_PARAM; i[3]++)
    {
      start_pos[wavelet_nr][3]=i[3];
      for (i[4]=MIN_PARAM; i[4]<=MAX_PARAM; i[4]++)
      {
        start_pos[wavelet_nr][4]=i[4];
        for (i[5]=MIN_PARAM; i[5]<=MAX_PARAM; i[5]++)
        {
          start_pos[wavelet_nr][5]=i[5];
          
          // refresh image:
          for (j=0; j<nr_pixels; j++) image[j]=data_image[j];
          
          // transform & retransform:
          Construction(wavelet_nr,nr_pixels,0,1,image,buffer);               
          Reconstruction(wavelet_nr,nr_pixels,0,1,image,buffer);   

          // test retransformated image only in area not 
          // influenced by mirror: from  AVOID_MIRROR to (nr_pixels-AVOID_MIRROR)
          ok=1;               
          for (j=AVOID_MIRROR; j<(nr_pixels-AVOID_MIRROR); j++)
            if (image[j] != data_image[j]) ok=0;
            
          // 'good' image found: print start_pos
          if (ok)
          {
            printf("wavelet:%d -> start_pos: ",wavelet_nr);
            for (j=2; j<6; j++) printf("%d ",i[j]);
            printf("\n");
          }
          
          // ...and next permutation...
        }
      }
    }
  }
}


int main(void)
{
  // CODE:
  Wave_Init(0);
  
//  SearchParam();
//  for (int i=-10; i<8; i++) mirror2(1,i,4,8,0,0,NULL);

  for (int i=0; i<nr_pixels; i++) image[i]=data_image[i];

  printf("Image:\n");  
  for (int i=0; i<nr_pixels; i++) printf("%f ",image[i]);
  printf("\n\n");

  Construction(wavelet_nr,nr_pixels,0,1,image,buffer); 

  printf("Transformed Image:\n");
  for (int i=0; i<nr_pixels; i++) printf("%f ",image[i]);
  printf("\n\n");
  
  Reconstruction(wavelet_nr,nr_pixels,0,1,image,buffer);   

  printf("Re-Transformed Image:\n");
  for (int i=0; i<nr_pixels; i++) printf("%f ",image[i]);
  printf("\n\n");
  //(void)getchar();
  
  return 0;
}
*/