/*****************************************************************************
File:           train_nsip97.c   November 9
Authors:        Bogdan Dumitrescu 
                Ioan Tabus   (tabus@cs.tut.fi)

Purpose:        Design of a stack filter starting from costs.


Copyright 1998   Bogdan Dumitrescu, Ioan Tabus.
		All Rights Reserved.


Reference:
@inproceedings{yoo,
   author = { K.L.Fong and G.B-Adams III and E.J.Coyle and J.Yoo  },
   title = {Synthesis of a Parallel Optimal Stack Filter Training Algorithm},
   booktitle = "Nonlinear Signal and Image Processing",
   address = "",
   volume = "1",
   year = "1997",
   pages = "paper 222",
}



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.

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

/****************************************************************************/
/*                                                                          */
/*                  Fast (sub)optimal stack filter design                   */
/*                                                                          */
/*  Inputs:                                                                 */
/*    - name of coefficients file  (e.g. as written by coeff1.c)            */
/*    - name of output filter file                                          */
/*  Output:                                                                 */
/*    - filter file  (it will be written in boolean tabular form)           */
/*                                                                          */
/*  Allows insertion of: - desired window size                              */
/*                                                                          */
/*                                                                          */
/****************************************************************************/

#define NUMBER_ITER	100			/* number of iterations */
#define DEC_THR		1000000			/* decision threshold */


#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 NUM 258064

long D[LG];        /* Decision variables                                */

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

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

double J;	   /* mean average error				*/

clock_t clt;


void Yoo_iteration();
void Yoo();

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


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

{

   FILE *fp;
   char name[40], namef[40];
   int i,j,nr;


/* Read images from disk                                                */   

 printf("\n ===========================================");
 printf("\n Program for TRAIN_NSIP97 parallel algorithm");
 printf("\n ===========================================\n");  

  printf("\n Input window size: ");
  scanf("%d",&winsize);
  coefno=pow(2,winsize); /* Compute no. of necessary coefficents        */


if(argc==3) strcpy(name, argv[1]);
else { 
         printf("\n Name of coefficients file: ");
         scanf("%s",name);    
      }

if(argc==3) strcpy(namef, argv[2]);
else { 
         printf("\n Name of filter file: ");
         scanf("%s",namef);    
      }


/* Initialization of decision variables                                 */ 

  for(i=0;i<coefno;i++)
    D[i]=0;


/* Read cost coefficients                                               */

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

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

  fscanf(fp,"%d",&C0);
 
  fclose(fp);


/* Yoo's algorithm                                           */

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

  Yoo();
  printf("All done\n"); 
  clt = clock() - clt;
  printf("J = %f    time = %f\n", J, (double) clt/CLOCKS_PER_SEC);

  
/* writing filter description in file */

  fp=fopen(namef,"w");
  if(fp==NULL) {printf("\n Error writing filter\n"); return 0;}
   
  for (i=0 ; i<coefno ; i++)
  {
    j = ( D[i] <= 0 ) ? 1 : 0;
    fprintf(fp, "%1d\n", j); 
  }

  fclose(fp);

}



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

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");
}
 
/* ******************************************************************** */
/* One iteration of Yoo's algorithm */

void Yoo_iteration()

{  int nb,cccb;

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

   int d,            /* Desired pixel                                   */
       tmp;          /* temporary variable                              */ 


   /* update decision variables using costs */

   for ( k=0 ; k < coefno ; k++ )
   {
     D[k] += C[k];

     /* adapt decision variables to threshold */
     if ( D[k] >= DEC_THR )
       D[k] = DEC_THR - 1;
     else if ( D[k] < -DEC_THR )
       D[k] = - DEC_THR;
   }

   /* enforce stacking property */

   /* print_long_vector(D, coefno, "D_before: ");  */

   for ( k=0, t=1 ; k < winsize ; k++, t*=2 )	/* passes */
   {
     for ( i=0, j=2*t ; i < coefno ; i+=j )	/* one pass */
     {
       for ( s=0 ; s < t ; s++ )
       {
	 ix = i+s;
         if ( D[ix] < D[ix+t] )			/* stack violation */
         {
	   d = D[ix] + D[ix+t];
	   D[ix] = d / 2;
	   D[ix+t] = d - d/2;
	 }
       }
     }
   }
   /* print_long_vector(D, coefno, "D_after:  ");  */
 }


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

/* Evaluate J(c,cost) starting from  canonical terms                       */
double evalJc( cof)
unsigned char cof[];
 
{
 unsigned i;
 double Jloc;
 
 Jloc= (double) C0;
 for ( i=0 ; i < coefno ; i++ )
  if ( D[i] <= 0 )
    Jloc += (double) C[i];
 return Jloc / NUM;
}
 

/* ******************************************************************** */
/* Yoo's algorithm */

void Yoo()
{
  int k;
  clock_t cl1;

  for ( k=0 ; k < NUMBER_ITER ; k++ )
  {
    /* cl1 = clock(); */
    Yoo_iteration();			/* perform iteration */
    /* clt += clock() - cl1; */
    J = evalJc();
    printf("it:%3d   J = %f\n", k+1, J);
  }
}