/*****************************************************************************
File:           coeff1.c   November 9
Authors:        Alin Alecu    
                Ioan Tabus   (tabus@cs.tut.fi)
                Doina Petrescu  
		Bogdan Dumitrescu

Purpose:        Fast computation of cost coefficients


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


Reference:
@ARTICLE{TFI,
   author = { I. T\u{a}bu\c{s} and D. Petrescu and M. Gabbouj  },
   title = {A Training Framework for  Stack and {Boolean} Filtering-{F}ast optimal design procedures and robustness case study},
   journal = TRIP # " Special Issue on Nonlinear Image Processing",
   volume = "IP-5",
   year = "1996",
   month = "June",
   pages = "809--826",
}



These programs are supplied free of charge for research purposes only,
and may not be 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.

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

/****************************************************************************/
/*                                                                          */
/*                  Computation of cost coefficients                        */ 
/*                                                                          */
/*  Inputs:                                                                 */
/*    - desired image                                                       */
/*    - noisy image or images (if several images, their names should end    */
/*      with the suffix "_00", "_01" etc.; the number of such images will   */
/*      be automatically detected)                                          */
/*  Output:                                                                 */ 
/*    - costs file                                                          */
/*                                                                          */
/*  Allows insertion of: - desired window size                              */
/*                       - window grid format                               */
/*                                                                          */
/*                                                                          */
/****************************************************************************/

#define		MERGE_SORT		/* undefined: insert_sort */

#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 StartPR(i,j) (i==2 && j==2)
#define NewLinePR(j) (j==2)

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

long C[LG];        /* Coefficients                                      */
long C0;           /* Coefficient C0                                    */

unsigned char imxr[IR][IC], /* Pixels read from noisy image             */
         imdr[IR][IC];      /* Pixels read from desired image           */
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                               */

clock_t clt;


void Estimate_N0_N1(); /* Actual algorithm to obtain coefficients       */
void GetWindow();  /* Function to obtain user selected window grid      */

#ifdef MERGE_SORT
void	Merge( int *vec1, int *ind1, int len1,		/* first vector  */
	       int *vec2, int *ind2, int len2,		/* second vector */
	       int *vec,  int *ind );			/* merged vector */
#endif

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


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

{

   FILE *fpx,*fpd,*fp;
   char name[40], namex[40], named[40], namexn[40];
   int i,j,one_image, ni;


/* Read images from disk                                                */   

 printf("\n ===========================================");
 printf("\n Compute cost coefficients for filter design");
 printf("\n ===========================================\n");  

 if(argc==4) strcpy(named, argv[1]);
 else {   
         printf("\n Name of desired image: ");
         scanf("%s",named);   
      }  
 if(argc==4) strcpy(namex, argv[2]);
 else {  
         printf("\n Name of (basic) noisy image: ");
         scanf("%s",namex);  
      }
   
if(argc==4) strcpy(name, argv[3]);
else { 
         printf("\n Name of coefficients file: ");
         scanf("%s",name);    
      }

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

  strcpy( namexn, namex );
  fpx=fopen(namex,"rb");
  if ( fpx==NULL )		/* noisy file not found, maybe it's numbered */
  {
    strcat( namexn, "_0" );
    fpx = fopen(namexn,"rb");
    if ( fpx == NULL )		/* this time it's an error */
    {
      printf("\n Error reading %s\n", namex);
      fclose(fpd);
      return 0;
    }
    one_image = 0;
  }
  else				/* there will be only one noisy file */
    one_image = 1;


/* Selection of window grid                                             */

  GetWindow();
  


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



/* Initialization of coefficients                                       */ 

  for(i=0;i<coefno;i++) {C[i]=0; }
  C0=0;


/* Computation of cost coefficients, need for MAE evaluation            */

  printf("\nStarting to compute...........\n"); 
  clt=clock();

  ni = 1;
  printf("File: %s\n", namexn);

  while ( 1 )				/* loop over all images */
  {
    for ( i=0 ; i<IR ; i++ )		/* read noisy image */
      fread((char *)imxr[i],IC,1,fpx);

    Estimate_N0_N1();

    fclose(fpx); 
    if ( one_image )
      break;

    /* generate next file name */
    strcpy( namexn, namex );
    strcat(namexn, "_");
    sprintf( namexn+strlen(namexn), "%d", ni );

    printf("File: %s\n", namexn);

    fpx = fopen(namexn,"rb");
    if ( fpx == NULL )			/* no more files */
      break;

    ++ni;
  }

  fclose(fpx);

  clt = clock() - clt;
  printf("Files processed: %d\n", ni); 
  printf("time = %f\n", (double) clt/CLOCKS_PER_SEC);

  
/* Writing coefficients to disk file                                    */

  fp=fopen(name,"w");
  if(fp==NULL) {printf("\n Error writing coefficients\n"); return 0;}
   
  for (i=0 ; i<coefno ; i++)
    fprintf(fp, "%d\n", C[i]); 

  fprintf(fp,"%d\n", C0);

  fclose(fp);

}



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

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    print_long_vector( long *v, int len, char *mes )
{
  int i;
 
  printf("%s", mes);
  for ( i=0 ; i < len ; i++ )
    printf("%6ld ", v[i]);
  printf("\n");
}
 

/* ******************************************************************** */
 
void Estimate_N0_N1()
 
{  int nb,numar,cccb;
 
   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         */
       *vect1,       /* = vect+1                                        */
       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                              */
 
  #ifdef MERGE_SORT
   int vect_old[N+2];
   int vect_new[N+2];  /* used in conjunction with vect */
   int ind_old[N];
   int ind_new[N];
   int no_old;         /* number of old points after window sliding     */
  #endif
 
   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;
   vect1 = vect + 1;
 
#ifdef MERGE_SORT
   no_old = winsize - nohnew;
#endif

/* Start scanning whole image                                           */
 
   for(i=2;i<IR-2;i++)
 
       for(j=2;j<IC-2;j++)
         {
 
           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=0;
             for(k=0;k<winsize;k++) {
               vect1[k]=(int)(*(pim + num[k]));
               indici[k]=k;
             }
 
            #ifdef MERGE_SORT
             InsertSort(vect1, winsize, indici, 0);
            #endif
           }
 
           else
 
/* If newline, means we have vertical shift of window                   */
 
             if (NewLinePR(j)) {
 
               t=0;
               for(k=0;k<winsize;k++) {
                 tmp=lineindici[k];
                 if(iv[tmp]) {
                   vect1[t] = linevec[k];
                   indici[t] = corv[tmp];
                   t++;
                 }
               }
               start=t;
               for(k=0;k<novnew;k++) {
                 tmp=ivnew[k];
                 vect1[t]=(int)(*(pim + num[tmp]));
                 indici[t]=tmp;
                 t++;
               }
 
              #ifdef MERGE_SORT
               InsertSort(vect1, winsize, indici, start);
              #endif
             }
 
 
           else
/* Any other case, means we have horizontal shift of window             */
 
             {
 
             t=0;
            #ifndef MERGE_SORT
             for(k=0;k<winsize;k++) {
               tmp=indici[k];
               if(ih[tmp]) {
                 vect1[t]=vect1[k];
                 indici[t]=corh[tmp];
                 t++;
               }
             }
             start=t;
             for(k=0;k<nohnew;k++) {
               tmp=ihnew[k];
               vect1[t]=(int)(*(pim + num[tmp]));
               indici[t]=tmp;
               t++;
             }
            #else
             for ( k=0 ; k<winsize ; k++ )
             {
               tmp = indici[k];
               if ( ih[tmp] )
               {
                 vect_old[t] = vect1[k];
                 ind_old[t]  = corh[tmp];
                 t++;
               }
             }
             for ( k=0 ; k < nohnew ; k++ )
             {
               tmp = ihnew[k];
               vect_new[k] = (int)(*(pim + num[tmp]));
               ind_new[k]  = tmp;
             }
             InsertSort(vect_new, nohnew, ind_new, 0);
 
             Merge(vect_old, ind_old, no_old,
                   vect_new, ind_new, nohnew,
                   vect1, indici);
            #endif
 
           }
 
/* After loading pixel values into vect, now sort them                  */
 
          #ifndef MERGE_SORT
           InsertSort(vect1, winsize, indici, start);
          #endif
 
/* If newline or start case, keep in mind vect and indici vectors for   */
/* futher use                                                           */
 
           if (StartPR(i,j) || NewLinePR(j))
             for(k=0 ; k<winsize ; k++) {
               linevec[k] = vect1[k];
               lineindici[k] = indici[k];
             }
 
/* Obtain desired value for pixel                                       */
 
           dfromimd(i,j,&d);
 
/* Now compute C0                                                       */
 
           C0+=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];
 
 
               if(d >= vect[ix]) C[nb]-=numar;
               else
                 {
                   if (d < vect[ix -1]) C[nb]+=numar;
                   else
                     {
                       C[nb]-=(2*d - vect[ix-1] - vect[ix]);
                     }
                 }
 
               cccb = coefno / 2;
               /* cccb >>= indici[ix-1]; */
               cccb >>= (winsize - indici[ix-1] - 1);
               nb ^=cccb;
             }
 
 
     }
 
 }
 
 

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

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

{  

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

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

/* simplified version: same indices in vector and index vector */

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];

/* 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+1] = Index[j];
        }

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

    }
}

/* ******************************************************************** */
  
#ifdef MERGE_SORT

/* merge two increasing sequences */
void	Merge( int *vec1, int *ind1, int len1,		/* first vector */
	       int *vec2, int *ind2, int len2,		/* second vector */
	       int *vec,  int *ind )			/* merged vector */
{
  int k1, k2, k;
  int *v, *i;

  k1 = k2 = 0;
  while ( k1 < len1 && k2 < len2 )
  {
    if ( CompGT( *vec1, *vec2 ) )
    {
      *vec++ = *vec2++;
      *ind++ = *ind2++;
      k2++;
    }
    else
    {
      *vec++ = *vec1++;
      *ind++ = *ind1++;
      k1++;
    }
  }

  if ( k1 == len1 )
  {
    v = vec2;
    i = ind2;
    k = len2 - k2;
  }
  else
  {
    v = vec1;
    i = ind1;
    k = len1 - k1;
  }

  while ( k-- )
  {
    *vec++ = *v++;
    *ind++ = *i++;
  }
}

#endif