#include <stdio.h>
#include <malloc.h>

#define dimension_x   512
#define dimension_y   512
#define levels        6

#ifndef CLOCKS_PER_SEC
#define CLOCKS_PER_SEC 1e6
#endif

int NumbTap=4;
float R_LowPass[5] =
  { 0.852699,  0.418092, -0.110624, -0.064539, 0.037829 };

float R_HighPass[5] =
  { 0.788485, -0.377403, -0.040690,  0.023849, 0.0 };

float coeff[dimension_x][dimension_y], Filter_L(), Filter_H();


main(argc, argv)
int     argc;
char    *argv[];
{
FILE    *in_file,*out_file;
int i,j;
int t,xrint();

  if (argc != 3)
    {
    printf("Usage: %s image_file out_file\n", argv[0]);
    exit(0); 
    }

  in_file=fopen(argv[1],"r");  out_file = fopen(argv[2],"w");

  fread(coeff,sizeof(float),dimension_x*dimension_y,in_file);

  image_recover();

for(i=0;i<dimension_x;i++)
  for(j=0;j<dimension_y;j++) 
     {
     t=xrint(coeff[i][j]);
     if(t<0) t=0;
     if(t>255) t=255;
     putc((unsigned char)t,out_file);
     }
}


image_recover()
{
int i,j,k;
int lv,hr,mn;
int nx,ny,mx,my,mark;
float *temp_line, *t, *in_line, *out_line,sc;

printf("Starting inverse transformation...\n");
mark=clock();
i = NumbTap, k = Max(dimension_x, dimension_y) + (i << 1);

if( (temp_line=(float *)calloc(k, sizeof(float))) == NULL ) {
    fprintf(stderr, "==> calloc() error");
    exit(-1);
    }

in_line = temp_line + i;
out_line = in_line;

nx = dimension_x>>levels, ny = dimension_y>>levels;

for (lv = 0; lv < levels; lv++) {
  mx = nx;  nx <<= 1;  my = ny;  ny <<= 1;
    for (i = 0; i < nx; i++) {
      for (j = 0, t = in_line; j < my; j++) {
        *(t++) = coeff[i][j];  *(t++) = coeff[i][j+my]; }
      Reflection(in_line, in_line + ny - 1);
      for (j = 0, t = in_line; j < ny;) {
        coeff[i][j++] = Filter_H(R_HighPass, t++);
        coeff[i][j++] = Filter_L(R_LowPass, t++); } }

    for (j = 0; j < ny; j++) {
      for (i = 0, t = in_line; i < mx; i++) {
        *(t++) = coeff[i][j];  *(t++) = coeff[i+mx][j]; }
      Reflection(in_line, in_line + nx - 1);
      for (i = 0, t = in_line; i < nx;) {
        coeff[i++][j] = Filter_H(R_HighPass, t++);
        coeff[i++][j] = Filter_L(R_LowPass, t++); } }
}

  sc=(float)(clock()-mark)/CLOCKS_PER_SEC;
  printf("Image recovered in");
  hr = (int) (sc / 3600.0);  sc -= 3600.0 * hr;
  mn = (int) (sc / 60.0);  sc -= 60.0 * mn;
  if (hr) {
    printf("%d hour", hr);
    if (hr > 1) printf("s, "); else printf(", "); }
  if ((hr) || (mn)) {
    printf("%d minute", mn);
    if (mn > 1) printf("s, and "); else printf(", and "); }
  printf("%5.2f seconds.\n", sc);

}



int Max(a,b)
int a,b;
{
return (a > b ? a : b);
}

float Filter_L(f,v)
float * f, *v; 
{
  return f[0] * v[0] +
    f[1] * (v[1] + v[-1]) + f[2] * (v[2] + v[-2]) +
    f[3] * (v[3] + v[-3]) + f[4] * (v[4] + v[-4]);
}


float Filter_H(f, v)
float * f, *v; 
{
  return f[0] * v[0] +
    f[1] * (v[1] + v[-1]) + f[2] * (v[2] + v[-2]) +
    f[3] * (v[3] + v[-3]);
}


Reflection(h,t)
float * h, * t;
{
int i;
  for (i = 1; i <= NumbTap; i++) {
    h[-i] = h[i];  t[i] = t[-i]; }
}

int xrint(x)
float x;
{
if(x>=0.) return ((int)(x+.5));
  else return ((int)(x-.5));
}