/*CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC*/
/*                                                                      */
/* REMOVES NOISE OF AN IMAGE PRESERVING TEXTURE AND ENHANCING SHARPNESS */
/* USAGE: smoother  image M smoothing sharpening                        */
/*         smoother is the name of the executable                       */
/*         image is the file name of the raw ppm image (with no .ppm)   */
/*         2<=M<=128 is the local window size (8 is recommended)        */
/*         smoothing>=0 is the desired smoothing (150 is recommended)   */
/*         1<=sharpness is the sharpness enhancement (1.1 is rec.)      */
/*                                                                      */
/*         Copyright (c) 2000, Santiago I. Betelu                       */
/*         Registration Number:  TXu-966-423                            */ 
/*CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC*/ 

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

/************************************************************************/
/* THIS ROUTINE SMOOTHES A TWO DIMENSIONAL IMAGE OF DIMENSIONS          */ 
/* width x height STORED SEQUENTIALLY IN THE ARRAY inp[width * height]; */
/* THE COORDINATES (i,j) ARE IN inp[i*width+j].                         */
/* colorlevel IS THE MAXIMUM PERMISIBLE INTENSITY VALUE OF THE IMAGE    */
/* M IS THE SIZE OF THE LOCAL SMOOTHING WINDOW, WHICH HAS THE SIZE OF   */
/* THE TEXTURE LENGTHSCALE (M=8 IS RECOMMENDED IF THE TEXTURE           */
/* CHARACTERISTICS ARE UNKNOWN).                                        */
/* smoothing IS THE DEGREE OF SMOOTHING, AND IS ABOUT 10 TIMES THE      */
/* STANDARD DEVIATION OF THE NOISE.                                     */
/* IF THE NOISE IS UNKNOWN, TRY TENTATIVELY WITH noise=1, 10, 100 ...   */
/* FOR A FIZED VALUE sharpening=1.0.                                    */
/* sharpening IS THE SHARPENING ENHANCEMENT. sharpening>1               */
/* SHARPENS THE IMAGE AND sharpening<1 BLURS IT.                        */
/************************************************************************/
void smooth(unsigned char inp[], int width, int height, int colorlevel,
            int &M, float smoothing, float sharpening, int compo);

/***********************************************************************/
/* ANY STANDARD                                                        */
/* TWO DIMENSIONAL MxM NON-NORMALIZED COMPLEX FOURIER TRANSFORM OF     */
/* THE TWO-DIMENSIONAL ARRAY window WITH REAL PART IN window[i][j] AND */
/* IMAGINARY PART IN window [i][j+M].                                  */
/* SIZE OF window IS 128x256. RETURNS THE RESULT IN THE SAME ARRAY.    */
/* M= 2**logM AND direction=1 IS FORWARD TRANSFORM, direction=-1 IS    */
/* THE INVERSE TRANSFORM                                               */
/***********************************************************************/
void ft(float window[128][256],
             int M, int logM, int direction);

/**********************************************************************/
/* ANY STANDARD                                                       */
/* FOURIER TRANSFORM FOR ARRAYS OF SIZE M=2**logM                     */
/* REAL PART IS STORED IN a[0], a[1], ... a[M-1]                      */
/* IMAGINARY PART IS IN   a[M], a[M+1], ... a[2*M-1]                  */
/* RETURNS THE NON-NORMALIZED COMPLEX FOURIER TRANSFORM IN a[]        */
/**********************************************************************/
void fastfourier(float a[], int M, int logM, int direction);

/*********************************************************************/
/* MAIN READS AN IMAGE, FILTERS IT AND SAVES THE RESULT              */
/*********************************************************************/
int main(int numberargs, char** args) {
  unsigned char *inp, *comp;
  int width,height,colorlevel,M,i;
  FILE  *fp;
  char renglon[82],name[82];
  float smoothing,sharpening;
     
/* READS THE ARGUMENTS */
  if(numberargs!=5){
    printf("Arguments: \n");
    printf("           image name (with no .ppm termination)");
    printf("           window size (power of 2, rec 8)\n");
    printf("           smoothing strength 0..infinity, (rec 150)\n");
    printf("           sharpening, (rec 1.1)\n");
    exit(1);
  }
  M=atoi(args[2]);
  if(M>128){
    printf("M must be <=128 \n");
    exit(1); 
  }
  if(M==-1)
     printf("Copyright (c) 2000, Santiago I. Betelu, all rights reserved\n");
  smoothing= atof(args[3]);
  sharpening=  atof(args[4]);
/* READS THE IMAGE */
  strcpy(name,args[1]);
  strcat(name,".ppm");  
  fp = fopen(name, "r");
  if (fp==NULL) { 
    printf("Error opening file %s\n",name);
    return 0;
  }
  fgets(renglon, 82, fp);
  if(strncmp(renglon,"P6",2)==0){
     fgets(renglon, 82, fp);
     while (renglon[0] == '#') fgets(renglon, 82, fp);
     sscanf(renglon, "%d %d\n", &width, &height);
     printf("width=%d height=%d   ", width, height);
     fscanf(fp, "%d", &colorlevel);
     printf("max colorlevel is %d\n", colorlevel);
     fgets(renglon, 82, fp);
  }
  else {
     printf("the image is not in ppm format\n");
     fclose(fp);
  }
  comp= new unsigned char[width*height];
  inp= new unsigned char[3*width*height];
  if (inp==NULL||comp==NULL) {
    printf("No enough memory\n");
    fclose(fp);
    return 0;
  }
  fread(inp, sizeof(unsigned char), 3*width*height, fp);
  fclose(fp); 
/*END READING IMAGE */
      
  for(i=0;i<width*height;i++) comp[i]=inp[3*i];
  smooth(comp,width,height,colorlevel,M,smoothing,sharpening,0);
  for(i=0;i<width*height;i++) inp[3*i]=comp[i];
  
  for(i=0;i<width*height;i++) comp[i]=inp[3*i+1];
  smooth(comp,width,height,colorlevel,M,smoothing,sharpening,1);
  for(i=0;i<width*height;i++) inp[3*i+1]=comp[i];

  for(i=0;i<width*height;i++) comp[i]=inp[3*i+2];
  smooth(comp,width,height,colorlevel,M,smoothing,sharpening,2);
  for(i=0;i<width*height;i++) inp[3*i+2]=comp[i];
      
/* OUTPUT */
  strcpy(name,args[1]);
  strcat(name,"_s.ppm");
  fp = fopen(name, "w");
  if(fp == NULL) {
    printf("Error opening output file\n");
    return 0;
  }  
  fprintf(fp, "P6\n");
  fprintf(fp, "%d %d\n", width, height);
  fprintf(fp, "%d\n", colorlevel);
  fwrite(inp, sizeof(unsigned char), 3*width*height, fp);
  fclose(fp);
/*END OUTPUT */
  printf("                     *** DONE ***\n");
  return 0;
}
/********************************************************************/
void smooth(unsigned char inp[], int width, int height, int colorlevel, 
            int &M, float smoothing, float sharpening,int compo)
{
 int i,j,i1,j1,i2,j2,kx,ky,num,in,sib,logM,ky0;
 float norm,q,s,w;
 float *result;
 float window[128][256];

/* INITIALIZATION FOR FFT */
  logM=  int((log( (double) M )/log( 2.0 ))+0.5);
  M= int(pow(2,logM));   

  result= new float[width*height];
  if (result==NULL) {
    printf("No enough memory\n");
    exit(1);
  }
  for(in=0;in<width*height;in++)result[in]=0;

/* SCAN IMAGE */
  for(i=0; i<=height; i=i+M/2) {
    printf("\rCompleted %4.1f percent",compo*33+34.*i/(M/2*(height/(M/2))));
    for(j=0; j<=width; j=j+M/2) {
/*         FIRST COMPUTES THE LOCAL FOURIER COEFFICIENTS AND STORES */ 
/*         THEM IN window                                           */
           for(i1=0;i1<M;i1++){
/*           REFLECTION AT BORDERS */
             i2=abs(i+i1-M/2);
             if(i2>=height)i2= height-2-(i2-height);
             for(j1=0;j1<M;j1++){
                j2=abs(j+j1-M/2);
                if(j2>=width)j2= width-2-(j2-width);
                window[i1][j1]=inp[i2*width+j2];
                window[i1][j1+M]=0;
             }
           }
           ft(window,M,logM,1);

/*         COMPUTES SUPER-SIMPLYFIED PSEUDO-LOCAL WIENER FILTER */
           q=smoothing*M*M+1.e-6;
           s=sharpening;
/*         ky0 IS TO AVOID PROCESSING kx=ky=0  */
           ky0=1;
           for(kx=0; kx<M; kx++){
              for(ky=ky0; ky<M; ky++){
                 norm= window[kx][ky]*window[kx][ky]+
                       window[kx][ky+M]*window[kx][ky+M];
                 window[kx][ky]=window[kx][ky]*norm/(q+norm)*s;
                 window[kx][ky+M]=window[kx][ky+M]*norm/(q+norm)*s;
              }
              ky0=0;
           }
/*         RETURNS TO THE SPATIAL DOMAIN */
           ft(window,M,logM,-1);
           sib= M/2;
           for(i1=0;i1<M;i1++){
             for(j1=0;j1<M;j1++){
                i2= i+i1-M/2;
                j2= j+j1-M/2;
                sib=sib/2;
/*              ITERATES IN THE WEIGHTED RESULT */
                if(i2>=0&&i2<height&&j2>=0&&j2<width){
                  result[i2*width+j2]= window[i1][j1]/M/M*(M/2-abs(i1-M/2))*
                                     (M/2-abs(j1-M/2))+result[i2*width+j2];
                }
             }
           }
/*  END LOOPS i,j */
    }         
  }
/* NORMALIZES THE RESULT */
  for(in=0;in<width*height;in++){
    w=4*result[in]/M/M+0.5;
    if(w>colorlevel)w=colorlevel;
    if(w<0)w=0;
    inp[in]=int(w);
  }
}

/******************************************************************/
void fastfourier(float f[], int M, int logM, int direction){
/* THIS FFT WORKS FOR ARRAYS
   f[0], f[1], ..., f[M-1];     M=2**logM
*/
 int i,j,k,l,nl,halfnl,in;
 float ur,ui,wr,wi,tr,ti;
 const float pi=3.141592654;

 j=0;
 for(i=0;i<M-1;i++){   
/*  PUTS THE VECTOR f IN THE ORDER OF REVERSED BITS OF i */
    if(i<j){
      tr=f[j];
      f[j]=f[i];
      f[i]=tr;
      ti=f[j+M];
      f[j+M]=f[i+M];
      f[i+M]=ti;
    }
    k=M/2;
/*  COUNTS, AND WRITES THE COUNT IN REVERSE ORDER IN j (j=reverse i) */
    while(k<=j){
      j=j-k;
      k=k/2;
    }
    j=j+k;
 }
/* RECURSIVE FOURIER SUMS*/
 nl=1;
 for(l=0;l<logM;l++){
    nl=nl*2;
    halfnl= nl/2;
    ur=1.0;
    ui=0.0;
    wr= cos(-pi/halfnl);
    wi= sin(direction*pi/halfnl);
    for(j=0;j<halfnl;j++){
       for(i=j;i<M;i=i+nl){
          in=i+halfnl;
          tr=-f[in]*ur+f[in+M]*ui;
          ti=-f[in+M]*ur-f[in]*ui;
          f[in]=f[i]+tr;
          f[in+M]=f[i+M]+ti;  
          f[i]=f[i]-tr;
          f[i+M]=f[i+M]-ti;
       }
       tr=ur*wr-ui*wi;
       ui=ui*wr+ur*wi;
       ur=tr;
    }
 }
 return;
}
/**************************************************************/
void ft(float window[128][256],
             int M, int logM, int direction){
/* VERY CRUDE AND DIRECT COMPUTATION OF THE 2D FFT */
 int  i, j;
 float aux[256];

/* INTEGRAL IN x FOR y CONSTANT */
 for(i=0; i<M;i++){
   for(j=0;j<M;j++){
        aux[j]  = window[i][j];
        aux[j+M]= window[i][j+M];
   }
   fastfourier(aux,M,logM,direction);
   for( j = 0; j < M; j++ ){
        window[i][j] = aux[j];
        window[i][j+M]=aux[j+M];
   }
 }
/* INTEGRAL IN y OF THE INTEGRAL IN x */
 for(i=0;i<M;i++){
     for(j=0;j<M;j++){   
        aux[j]  = window[j][i];
        aux[j+M]= window[j][i+M];
     }
     fastfourier(aux,M,logM,direction);
     for( j = 0; j < M; j++ ){
        window[j][i] = aux[j];
        window[j][i+M]=aux[j+M];
     }
 }
}
/******************************************************************/