#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 T_LowPass[5] =
  { 0.852699,  0.377403, -0.110624, -0.023849, 0.037829 };

float T_HighPass[5] =
  { 0.788485, -0.418092, -0.040690,  0.064539, 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;

  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");

  for(i=0;i<dimension_x;i++)
     for(j=0;j<dimension_y;j++) coeff[i][j]=(float)getc(in_file);
  image_transform();

  /*for(i=0;i<dimension_x;i++)
    for(j=0;j<dimension_y;j++) putc((unsigned char)(coeff[i][j]/2),out_file);*/

fwrite(coeff,sizeof(float),dimension_x*dimension_y,out_file);

}


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

printf("Starting image 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;

mx = dimension_x, my = dimension_y;

  for (lv = 0; lv < levels; lv++) {

    nx = mx;  mx >>= 1;  ny = my;  my >>= 1;

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

    for (i = 0; i < nx; i++) {
      memcpy(in_line, coeff[i], ny * sizeof(float));
      Reflection(in_line, in_line + ny - 1);
      for (j = 0, t = in_line; j < my; j++) {
        coeff[i][j] = Filter_L(T_LowPass, t++);
        coeff[i][j+my] = Filter_H(T_HighPass, t++); } } 
}

sc=(float)(clock()-mark)/CLOCKS_PER_SEC;
printf("Image transformed 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]; }
}