/*****************************************************************************
File:           filter.c   July 10
Authors:        Alin Alecu    
                Ioan Tabus   (tabus@cs.tut.fi)
                Doina Petrescu  

Purpose:        Fast stack and Boolean filtering 


Copyright 1997   Alin Alecu, Ioan Tabus, Doina Petrescu. All Rights 
                 Reserved.

These programs are supplied free of charge for research purposes only,
and may not sold or incorporated into any commercial product.  There is
ABSOLUTELY NO WARRANTY of any sort, nor any undertaking that they are
fit for ANY PURPOSE WHATSOEVER.  Use them at your own risk.  If you do
happen to find a bug, or have modifications to suggest, please report
it to  tabus@cs.tut.fi.  The copyright notice above and this statement of conditions 
must remain an integral part of each and every copy made of this file.

/****************************************************************************/
/*                                                                          */
/*                     Filtering of Corrupted Image                         */ 
/*                                                                          */
/*  Inputs:                                                                 */
/*    - desired image                                                       */
/*    - noisy image                                                         */
/*    - filter file                                                         */ 
/*  Output:                                                                 */ 
/*    - image file                                                          */
/*                                                                          */
/*  Requires insertion of: - desired window size                            */
/*                       - window grid format                               */
/*                                                                          */
/*                                                                          */
/****************************************************************************/


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

#define N  21   /* Max. size of window.
                      Eg:  5   9   10     17     19       21  */
 
#define LG 2097152  /* Max. no. of coefficients.
                      Eg: 32 512 1024 131072 524288  2097152 */
 
#define IR 512     /* Dimensions in pixels                              */
#define IC 512     /*    of picture                                     */ 

#define NUM 258064

#define StartPR(i,j) (i==2 && j==2)
#define NewLinePR(j) (j==2)

#define CompGT(a,b) (a > b)


int f_opt[LG];     /* Boolean filter                                    */
double MAE;        /* Mean average error                                */

unsigned char imxr[IR][IC], /* Pixels read from noisy image             */
              imdr[IR][IC], /* Pixels read from desired image           */
              imfr[IR][IC]; /* Filtered image pixels                    */
unsigned int num[N];   /* Pozition of selected pixel                    */
                       /* relative to window location                   */

int winsize;       /* Selected size of window                           */
int coefno;        /* No. of necessary coeff.                           */

int sel[N],        /* Selected pixels                                   */
    ih[N],iv[N],   /* Common pixel indexes belonging to old window      */
                   /* vector, corresponding to horizontal and vertical  */
                   /* movement  of window                               */
    corh[N],corv[N],  /* Common pixel indexes belonging to new window   */
                   /* vector, corresponding to horizontal and vertical  */
                   /* movement  of window                               */
    nohnew,novnew, /* Total no. of new pixels introduced in new window  */
                   /* vector, corresponding to horizontal and vertical  */
                   /* movement  of window                               */ 
    ihnew[N-1],ivnew[N-1];    /* New pixels introduced in new window    */    
                   /* vector, corresponding to horizontal and vertical  */
                   /* movement  of window                               */


void Filter_N0_N1(); /* Performs image filtering                        */
void GetWindow();  /* Function to obtain user selected window grid      */


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


void main( argc, argv)
int argc;
char *argv[];

{

   FILE *fpx,*fpd,*fp;
   char name[40],namex[40],named[40],namef[40];
   int i,j,nr;
   clock_t cl1=clock();


/* Read images from disk                                                */   

 printf("\n ===================== ");
 printf("\n Program for filtering ");
 printf("\n ===================== \n");  


/* Reading images */   

 printf("\n For computing MAE of result you'll need the desired image \n");

 if(argc==5) strcpy(named,argv[1]);
 else {  
         printf("\n Name of desired image: ");
         scanf("%s",named);
      }   
 if(argc==5) strcpy(namex, argv[2]);
 else {  
         printf("\n Name of noisy image: ");
         scanf("%s",namex);
      }
   
 if(argc==5) strcpy(name,argv[3]);
 else { 
         printf("\n Name of vector form of Boolean filter file: ");
         scanf("%s",name);
      }
  
 if(argc==5) strcpy(namef,argv[4]);
 else { 
         printf("\n The RESULT will be here: ");
         scanf("%s",namef);
      }

 fpd=fopen(named,"rb");
 if(fpd==NULL) {printf("\n Error reading %s\n",named); return;}

 fpx=fopen(namex,"rb");
 if(fpx==NULL) {printf("\n Error reading %s\n",namex); fclose(fpd);return;}


 for(i=0;i<IR;i++)
   { 
     fread((char *)imxr[i],IC,1,fpx);
     fread((char *)imdr[i],IC,1,fpd);
   }
    
 fclose(fpd);
 fclose(fpx); 


/* Selection of window grid                                             */

  GetWindow();
  
/* Read filter from disk                                                */

 fpd=fopen(name,"r");
 if(fpd==NULL) {printf("\n Can't open %s\n",name); return;}


 for(i=0;i<coefno;i++)
   { if(fscanf(fpd,"%d",&f_opt[i])!=1)
          { fclose(fpd); printf("\n Error reading %s\n",name); exit(-1);}
   }

 fclose(fpd);




/* Filter the image                                                     */

  printf("\nStarting to compute..........."); 
  MAE=0;
  Filter_N0_N1();
  printf("All done\n"); 

  
/* Print time ellapsed                                                  */ 
  printf("\nCPU time was %.3lf\n",((double)clock()-cl1)/CLOCKS_PER_SEC);  

  printf("\nMAE was %.5lf\n\n",MAE/NUM);  



/* Writing image to disk file                                          */

  fpd=fopen(namef,"wb");
    if(fpd==NULL) {
      printf("\n Error writing  %s\n",namef); 
      fclose(fpx); 
      return;
    }
  for(i=0;i<IR;i++)
          if(fwrite(imfr[i],IC,1,fpd)!=1) { 
            printf("\n Error writing  %s\n",namef);
            break;
          }

   fclose(fpd);
  

}



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

void GetWindow()
{ 
  int i,j,x,chosen[25];    
  char ch;
  div_t qx;

  int matchh,matchv;
  
  system("clear");

  printf("Now select the type of window you wish to use:\n");


/* Display window grid from which to choose                            */
 
  printf("\n     |----|----|----|----|----|"); 
  for(i=0;i<25;i++) 
    if ((i%5) == 0) 
      printf("\n     | %2d ",i);
    else 
      if ((i+1)%5 == 0) {
        printf("| %2d |",i);
        printf("\n     |----|----|----|----|----|"); 
      }
    else printf("| %2d ",i);

  
  do {
  printf("\n\n How many pixels do you wish to use ? (min 5, max %d) : ", N);
  scanf("%d",&winsize);
  } while( (winsize<5) || (winsize>N) );
 
  coefno=pow(2,winsize); /* Compute no. of necessary coefficents        */

  printf("\nPixel 0 : 12 (by default)\n");
  qx=div(12,5);
  num[0]=(qx.quot*IC) + qx.rem;
  sel[0]=12;

  for(i=0;i<=24;i++) chosen[i]=0;
  chosen[12]=1;

/* Input desired pixels                                                 */
    
  for(i=1;i<winsize;i++) {

    do {
      printf("\nPixel %d : ",i);
      scanf("%d",&x);
    } while( (x<0) || (x>24) || (chosen[x]));

    chosen[x]=1;
    qx=div(x,5); 
    num[i]= (qx.quot*IC) + qx.rem;
    sel[i]=x;
  }

  scanf("%c",&ch);
  printf("\nIs it OK ? (y/n) : ");
  scanf("%c",&ch);

  if(ch !='y') GetWindow();

/* If all selections are OK, determine which pixels will remain and      */  
/* which will leave the window while it moves                            */

  else {

    nohnew=0;novnew=0;

    for(i=0;i<winsize;i++) {

      matchh=0;matchv=0;;
      for(j=0;j<winsize;j++) {

        /*if (sel[i]==(sel[j]+1)) { ih[i]=1;corh[i]=j; }      
        if (sel[i]==(sel[j]+5)) { iv[i]=1;corv[i]=j; }
        if (!matchh && (sel[i]==(sel[j]-1))) matchh=1;
        if (!matchv && (sel[i]==(sel[j]-5))) matchv=1; */
 
        if ((sel[i]==(sel[j]+1))&&(sel[i]%5!=0)) { ih[i]=1;corh[i]=j; }
        if (sel[i]==(sel[j]+5)) { iv[i]=1;corv[i]=j; }
        if ((!matchh && (sel[i]==(sel[j]-1)))&&(sel[j]%5!=0)) matchh=1;
        if (!matchv && (sel[i]==(sel[j]-5))) matchv=1;

      }

      if (!matchh) { ihnew[nohnew]=i;nohnew++; }
      if (!matchv) { ivnew[novnew]=i;novnew++; }

    }

  }

  
}



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

void Filter_N0_N1()

{  int nb,numar,cccb;

   int filtrata;

   int ix,i,j,t,k;   /* Loop indexes                                    */

   int d,            /* Desired pixel                                   */
       vect[N+2],    /* Vector which stores pixels in current window,   */
                     /* and on which computations shall be made         */
       indici[N],    /* Corresponding indexes of these pixels           */ 
       linevec[N+2], /* The pixel vector at each newline                */  
       lineindici[N],/* The corresponding indexes of these pixels       */
       start,        /* Point from which sorting shall begin            */ 
       tmp;          /* temporary variable                              */ 

   unsigned char *pim; /* Value of pixel corresponding to pixel 0 on    */
                       /* user grid                                     */

   void dfromimd() ;   /* Function to obtain desired pixel              */
   void InsertSort();  /* Sorting function                              */
         

   vect[0]=0; vect[winsize+1]=255;

/* Start scanning whole image                                           */
   
   for(i=2;i<IR-2;i++)

       for(j=2;j<IC-2;j++)
         {
           
           filtrata=0;

           pim=imxr[i-2]+j-2;
           
/* If start of window movement, load in vect and indici values          */ 
/* respectively indexes of pixels in window                             */

           if(StartPR(i,j)) {
           
             start=1;
             for(k=0;k<winsize;k++) {         
               vect[k+1]=(int)(*(pim + num[k]));
               indici[k]=k;          
             }
        
           }

           else 

/* If newline, means we have vertical shift of window                   */

             if (NewLinePR(j)) {
               
               t=1;
               for(k=0;k<winsize;k++) {
                 tmp=lineindici[k];
                 if(iv[tmp]) {
                   vect[t]=linevec[k+1];
                   indici[t-1]=corv[tmp];
                   t++;
                 }
               }
               start=t;
               for(k=0;k<novnew;k++) {
                 tmp=ivnew[k];
                 vect[t]=(int)(*(pim + num[tmp]));
                 indici[t-1]=tmp;
                 t++;
               }
               
             }


           else                
/* Any other case, means we have horizontal shift of window             */

             {
                     
             t=1;
             for(k=0;k<winsize;k++) {
               tmp=indici[k];
               if(ih[tmp]) {
                 vect[t]=vect[k+1];
                 indici[t-1]=corh[tmp];
                 t++;
               }
             }
             start=t;
             for(k=0;k<nohnew;k++) {
               tmp=ihnew[k];
               vect[t]=(int)(*(pim + num[tmp]));
               indici[t-1]=tmp;
               t++;
             }

           }

/* After loading pixel values into vect, now sort them                  */
           
           InsertSort(vect, winsize, indici, start); 

/* If newline or start case, keep in mind vect and indici vectors for   */
/* futher use                                                           */

           if (StartPR(i,j) || NewLinePR(j))
             for(k=1;k<winsize+1;k++) {
               linevec[k]=vect[k];
               lineindici[k-1]=indici[k-1];
             }

/* Obtain desired value for pixel                                       */
           
           dfromimd(i,j,&d);           


/* Now compute coefficients based on already computed vectors vect and  */
/* indici                                                               */   
  
           nb = coefno - 1;

           for(ix=1;ix<winsize+2;ix++)
             {
               numar=vect[ix] - vect[ix -1];

         
               filtrata+=numar*f_opt[nb];
    
         
               cccb = coefno / 2;
               /*cccb >>= indici[ix-1]; */
               cccb >>= (winsize - indici[ix-1] - 1);
               nb ^=cccb;              
             }

         imfr[i][j]=filtrata;
         if((d-filtrata)>0) MAE+=(double)(d-filtrata); 
         else MAE-=(double)(d-filtrata);    

                
     }
               
 }


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

void dfromimd(ilin,icol,d)
  
 int *d ,ilin,icol;

{  

  *d=(int)(*(imdr[ilin]+icol));
}

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

void InsertSort(int Vector[N], int  n, int Index[N], int start) {
  int i,j,tmp,itmp;
  
    for (i = start; i <= n; i++) {
        tmp = Vector[i];
        itmp = Index[i-1];

/* Shift elements down until insertion point found.                     */

        for (j = i-1; j >=0 && CompGT(Vector[j],tmp); j--) { 
            Vector[j+1] = Vector[j];
            Index[j] = Index[j-1];
        }

        Vector[j+1] = tmp;
        Index[j] = itmp;

    }
}