/*************************************************************************/
/*                                                                       */
/*                         Wavelet Transform                             */
/*                                                                       */
/*                                    Date  : 1991. 11. 1.               */
/*                                    Programmer : 9125M15  Yoon, Doo-man*/
/*                                                                       */
/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
/*                                                                       */
/*   This program decompose the input image with wavelet transform, and  */
/*   then recover it to the original. And in addition, this detects      */
/*   the edge with Haar's wavelet. The intermediate results are saved    */
/*   on files with extensions ".xxx" and ".edge" ,respectively.          */
/*                                                                       */
/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
/*  # Compile :                                                          */
/*                                                                       */
/*    % cc < source filename> -lpixrect -lm -o <execution file>          */
/*                                                                       */
/*  # Run  :                                                             */
/*                                                                       */
/*    % wavelet                                                          */
/*                                                                       */
/*      < screen cleared >                                               */
/*                                                                       */
/*  # output :                                                           */
/*                                                                       */
/*   Wavelet Transform :                                                 */
/*          On screen :                                                  */
/*                                                                       */
/*            < Original image   Decomposed image   Composed image  >    */
/*                                                                       */
/*          On File :  < file_name".wxxx" >                              */
/*                           here, xxx means the size of decomposed unit */
/*                                                                       */
/*   Edge detection :                                                    */
/*          On screen :                                                  */
/*                                                                       */
/*            < Original image   Edge image  >                           */
/*                                                                       */
/*          On File :  < file_name".edge" >                              */
/*                                                                       */
/*************************************************************************/


#include <stdio.h>
#include <malloc.h>
#include <math.h>
#include <pixrect/pixrect_hs.h>
#define MAXSIZE 256 
#define HSIZE 500                            /* size of histogram  */
#define PMODE 0644
#define WHITE 255
#define BLACK 0 

char save_img ;

/***********************************************************************/
/*           Write the byte image to file                              */
/***********************************************************************/

write_img(img, img_name, sz)
unsigned char **img;
char *img_name;
int sz;
{
   int fd ;
   register int i,j ;
   char img1[MAXSIZE] ;

      /* create a file that will contain graph */
     if((fd=creat(img_name,PMODE)) == -1)
           error("write_img: "," can't create",img_name) ;

     for(i=0;i<sz;i++) {
        for(j=0;j<sz;j++)
          img1[j] = img[i][j] ;
        if(write(fd,img1,sz) != sz)
          error(" write error in routine write_img\n","","");
     } /* end of for */
     close(fd);
}



/*********************************************************************/
/*            Save the floating point image to File                  */
/*-------------------------------------------------------------------*/
/*************   img  : floating-point image             *************/
/*************   img_name: The output file name          *************/
/*************   sz   : File Size                        *************/
/*********************************************************************/

write_img1(img, img_name, sz)
float img[][MAXSIZE];
char img_name[20];
int sz;
{
   int fd ;
   register int i,j ;
   unsigned char img1[MAXSIZE] ;

      /* create a file that will contain graph */
     if((fd=creat(img_name,PMODE)) == -1)
           error("write_img: "," can't create",img_name) ;
     for(i=0;i<sz;i++) {
        for(j=0;j<sz;j++) {
          img1[j] = (unsigned char) img[i][j] ;
         }
        if(write(fd,img1,sz) != sz)
          error(" write error in routine write_img\n","","");
     } /* end of for */
     close(fd);
}

/************************************************************************/   
/* Purpose : allocates a picture matrix with range [nrl..nrh][ncl..nch] */
/* Calling Sequence : alloc_img(fsz,)                                   */
/* Algorithm :                                                          */
/* 	from Numerical Recepes in C, Cambridge University Press,        */
/*      pp 706                                                          */
/************************************************************************/
unsigned char **
alloc_img(fsz)
/*alloc_img(nrl, nrh, ncl, nch) */
	int  fsz ; 
{
	register int	i;
	int nrl,nrh,ncl,nch ;
	unsigned char **m;

	nrl = ncl = 0;
	nrh = nch = fsz-1;
	/* allocate pointers to rows. */
	m = (unsigned char **) calloc((nrh - nrl + 1),sizeof(unsigned char *));
	if (!m)
		error("alloc_img()","allocation failure","1");
	m -= nrl;

	/* allocate rows and set pointers to them */
	for (i = nrl; i <= nrh; i++) {
	      m[i] = (unsigned char *)
                      calloc((unsigned)(nch - ncl + 1),sizeof(unsigned char));
	      if (!m[i])
	       	      error("alloc_img()","allocation failure","4");
	      m[i] -= ncl;
	}

	/*** Return pointer to array of pointers to rows. */
	return m;
}


/************************************************************************/
/* Function : free_img                                                  */
/* Purpose : frees a picture matrix allocated with pmatrixi             */
/* Calling Sequence : free_img(m,fsz)                                   */
/* Algorithm :                                                          */
/*           from Numerical Recipes in C, Cambridge University Press,   */
/*           pp 708                                                     */
/************************************************************************/
free_img(m,fsz)
	unsigned char **m;
	int        fsz;
{
	int             nrl, nrh, ncl, nch;  
	register int	i;
	nrl = ncl = 0;
	nrh = nch = fsz-1;
	for (i = nrh; i >= nrl; i--)
		free((char *) (m[i] + ncl));
	free((char *) (m + nrl));
}

/*********************************************************************/
/***************   print error messages and exit    ******************/
/*********************************************************************/
  error(s1, s2, s3)

   char *s1, *s2 , *s3;
   {
      fprintf(stderr,"%s %s %s\07\n", s1, s2, s3) ;
      exit(1) ;
   }


/********************************************************************/
/* display binary image                                             */
/* Gray level images can be displyed on a color monotor with        */ 
/* proper look up table setting                                     */ 
/********************************************************************/
img_dsp1(img,fsz, vloc, hloc)
unsigned char **img ;
int fsz;
{
  Pixrect *screen;
  register int i,j ;

  screen = pr_open("/dev/fb");
     for(i=0;i<fsz;i++)
     for(j=0;j<fsz;j++)
        pr_put(screen,j+hloc,i+vloc,(int)img[i][j]);

     pr_close(screen);
  }


/********************************************************************/
/*               Display the floating point image on CRT            */
/*------------------------------------------------------------------*/
/*************    img  : image array                   **************/
/*************    fsz  : file size                     **************/
/*************    vloc : vertical location             **************/
/*************    hloc : horizontal location           **************/
/********************************************************************/
img_dsp2(img,fsz, vloc, hloc)
float img[][MAXSIZE] ;
int fsz;
{
  Pixrect *screen;
  register int i,j ;

  screen = pr_open("/dev/fb");
     for(i=0;i<fsz;i++)
     for(j=0;j<fsz;j++)
        pr_put(screen,j+hloc,i+vloc,(int)img[i][j]);

     pr_close(screen);
  }


/********************************************************************/
/***********    Read the image from file which is sz X sz    ********/
/********************************************************************/

unsigned char **
read_img(f_name, sz)
 char *f_name ;         /* input file name */
 int *sz ;              /* size of img to be written */

{  register int i,j ;
   unsigned char **img, **alloc_img();
   char img1[MAXSIZE] ;
   int fd , ssz,
   mode = 0;  /* mode 0 : read 
		    1 : write 
		    2 : read & write */

  if((fd = open(f_name,mode )) == -1 )
       error(" file :", f_name," open error ");

   i=lseek(fd,0L,2) ;
     lseek(fd,0L,0) ;
   ssz = (int)sqrt((double)i) ;
   *sz = ssz ;

   img = alloc_img(ssz);
   for(i=0;i<ssz;i++) {
      if(read(fd,img1,ssz) != ssz)
	   error("read error from file %s \n", f_name,"");
     for(j=0;j<ssz;j++)
       img[i][j] = img1[j] ;
   } /* end of for */
   close(fd);
   return(img);
}



/*******************************************************/
/*                 Function Menu                       */
/*-----------------------------------------------------*/
/************  fnum : function number        ***********/
/************  fname: Input image File Name  ***********/
/*******************************************************/
void menu(fnum,fname)
char *fnum, *fname ;
{
       system("clear");
       printf("\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n") ;
       printf("        ==========================================\n"); 
       printf("                      Wavelet Transform           \n");  
       printf("        ==========================================\n");       
       printf("          1. Decomposition & Composition\n");       
       printf("          2. Edge Detection          \n");        
       printf("          3. End                     \n");        
       printf("                                 Choose one : ");          
       scanf("%c",fnum) ;  getchar() ;
       printf("\n") ;
       if (*fnum > '2' || *fnum < '1' ) exit(1) ; 
       printf("          # Input file name :");                    
       scanf("%s",fname) ;  getchar() ;
       printf("          # Save the results on file (y/n)? "); 
       scanf("%c",&save_img) ;   getchar() ;
}


/**********************************************************************/
/*******************                                *******************/
/*******************        M   A   I   N           *******************/
/*******************                                *******************/
/**********************************************************************/

int h_loc = 0 , v_loc = 0 ;

void main(argc,argv)
int argc ;
char *argv[] ;
{
    int  fsz ;
    float THRESH;
    char fname[20], fnum, flag = 'y';
    unsigned char **img1, 
                  **read_img(), **alloc_img() ;
    float img2[MAXSIZE][MAXSIZE] ;
    void Composition(),Decomposition(),Detect_edge() ;
    float Compress_ratio() ;

    while (flag == 'y') {
       h_loc = 0 ; 
    
       menu(&fnum,fname) ;
       switch(fnum) {
        case '1' :   
          img1 = read_img(fname,&fsz) ;
          img_dsp1(img1,fsz, v_loc , h_loc) ; 
         
          Decomposition  (img1,img2,fsz,fname) ;          
          h_loc += fsz+20 ;                /* display the decomposed image */
          img_dsp2(img2,fsz, v_loc , h_loc) ; 
          printf("          # Threshod value :");
          scanf("%f",&THRESH); getchar() ;
          printf("           Compress ratio : %4.2f %\n",
                                     Compress_ratio(img2,fsz,THRESH)*100);
          Composition(img2,fsz,fname,THRESH) ;
          h_loc += fsz+20 ;                /* display the reconstructed image */
          img_dsp2(img2,fsz, v_loc , h_loc) ;
          break ;

       case '2' :
          img1 = read_img(fname,&fsz) ;
          img_dsp1(img1,fsz, v_loc , h_loc) ; 
          printf("          # Thresh value ? ") ;
          scanf("%f",&THRESH) ;  getchar() ;
          h_loc += fsz+20 ;                /* display the edge image */
          Detect_edge(img1,fsz,THRESH,fname) ;
          break ;
       default : exit(1) ;
      } /* switch */

      printf("\n          # Again ? ");
      scanf("%c",&flag) ; getchar() ;
      free_img(img1,fsz);

   }  /* while  */

}  /*  M A I N */



/*******************************************************/
/*******  Copy the pixel with floating-point type ******/
/*******************************************************/
void char2float (img1,img2,fsz) 
unsigned char **img1 ;
float img2[][MAXSIZE] ;
{
  int i,j ;

    for ( i= 0 ; i < fsz ; i ++)
      for ( j= 0 ; j < fsz ; j ++)
        img2[i][j] = (float)img1[i][j] ;
}




/***********************************************************/
/*                    Decomposition                        */ 
/*---------------------------------------------------------*/ 
/***********                                     ***********/ 
/***********   img1 : Oringinal Image            ***********/ 
/***********   img2 : Edge frame                 ***********/ 
/***********   fsz  : size of image              ***********/ 
/***********   fname: file name                  ***********/
/***********************************************************/

void Decomposition(img1,img2,fsz,fname)
unsigned char **img1 ;
float    img2[][MAXSIZE];
int fsz ;
char *fname ;
{
  int i, j, size ;
  char ch,outname[20], tail[20] ;
  float     tmp[MAXSIZE] ;

   char2float(img1,img2,fsz) ;

   for ( size = fsz  ; size > 2 ; size/=2 ) {
         for ( i = 0 ; i < size ; i ++)  {
            for ( j = 0 ; j < size ; j +=2 )  {
              tmp[j/2+size/2] = (img2[i][j] - img2[i][j+1]) / 2.0 ; /* H */
              tmp[j/2] = (img2[i][j]+img2[i][j+1]) / 2.0 ;          /* L */
            }               
            for ( j = 0 ; j < size ; j ++ )  
               img2[i][j] = tmp[j] ;
         }
             

         for ( j = 0 ; j < size ; j ++)  {
            for ( i = 0 ; i < size ; i +=2 )  {
              tmp[i/2+size/2] = (img2[i][j] - img2[i+1][j]) / 2.0 ; /* H */
              tmp[i/2] = (img2[i][j]+img2[i+1][j]) / 2.0 ;          /* L */
            }               
            for ( i = 0 ; i < size ; i ++ )  
               img2[i][j] = tmp[i] ;
         }
         /* ------- save the final result on  file ------- */
         if (save_img == 'y') {  
             sprintf(outname,"%s.w%d",fname,size) ;
             write_img1(img2,outname,fsz) ;
         } 
   } /* for size */
 } /*  Decomposition  */                                                                        
 
float Compress_ratio(img,fsz,THRESH)
float img[][MAXSIZE] ;
int fsz ;
float THRESH ;
{
  int i,j ;
  long count ; 
 
   count = 0 ;
   for ( i = 0 ; i < fsz ; i++)
      for ( j = 0 ; j < fsz ; j++)
        if ( fabs((double)img[i][j]) < (double)THRESH) {
                    count ++ ; 
                    img[i][j] = 0 ;
         }
   return( count/(float)(fsz*fsz) ) ;
}


/***********************************************************/
/*              Composition with Decomposed Image          */ 
/*---------------------------------------------------------*/ 
/***********   img1 : Oringinal Image            ***********/ 
/***********   img2 : Edge frame                 ***********/ 
/***********   fsz  : size of image              ***********/ 
/***********************************************************/
void Composition(img2,fsz,fname,THRESH)
float  img2[][MAXSIZE];
int    fsz ;
char *fname ;
float THRESH ;
{
  int i, j, size ;
  char outname[20] ;
  float tmp[MAXSIZE] ;
  long count ;


   for ( size = 2   ; size <= fsz/2 ; size*=2 ) {
         for ( j = 0 ; j < size*2 ; j ++ )  {
            for ( i = 0 ; i < size ; i ++) {
              tmp[i*2] = img2[i][j] + img2[i+size][j] ;   
              tmp[i*2+1] = img2[i][j] - img2[i+size][j] ;
            }               
            for ( i = 0 ; i < size*2 ; i ++ ) 
              img2[i][j] = tmp[i] ; 
         }               
         for ( i = 0 ; i < size*2 ; i ++) {
            for ( j = 0 ; j < size ; j ++ )  {
              tmp[j*2] = img2[i][j] + img2[i][j+size] ;   
              tmp[j*2+1] = img2[i][j] - img2[i][j+size] ;
            }               
            for ( j = 0 ; j < size*2 ; j ++ ) 
              img2[i][j] = tmp[j] ; 
         }               
   }
   if (save_img == 'y') {
               sprintf(outname,"%s.t%d",fname,(int)THRESH);
               write_img1(img2,outname,fsz) ;
   }
 }                                                                        


/***********************************************************/
/*                 Detect Edge with Haar's wavelet         */ 
/*---------------------------------------------------------*/ 
/**********    img1 : Oringinal Image             **********/ 
/**********    fsz  : size of image               **********/ 
/**********    THRESH : Threshold value           **********/ 
/**********    fname: File name                   **********/ 
/***********************************************************/
int pi1[][2] = { 1, -1,  1, -1 } ;  /*   Windows for edge detection */
int pi2[][2] = { 1,  1, -1, -1 } ;

void Detect_edge(img1,fsz,THRESH,fname)
unsigned char **img1 ;
int fsz ;
float THRESH ;
char *fname ;
{
   int i, j ;
   long tmp1,tmp2 ;
   unsigned char **edge ;
   char outname[20] ;
   int M[MAXSIZE][MAXSIZE] ; 

   edge = alloc_img(fsz);
   /*---------------------------------------------------------------*/
   /*                    Make Modulus Image                         */
   /*---------------------------------------------------------------*/
   for ( i = 0 ; i < fsz-1 ; i++)
     for ( j = 0 ; j < fsz-1 ; j++) {
        tmp1 = img1[i][j]*pi1[0][0]+img1[i][j+1]*pi1[0][1] + 
               img1[i+1][j]*pi1[1][0] + img1[i+1][j+1]*pi1[1][1] ; 
        tmp1 *= tmp1 ;

        tmp2 = img1[i][j]*pi2[0][0]+img1[i][j+1]*pi2[0][1] + 
               img1[i+1][j]*pi2[1][0] + img1[i+1][j+1]*pi2[1][1] ; 
        tmp2 *= tmp2 ;

        M[i][j] = (int)sqrt((double)(tmp1 + tmp2)) ;
     } 

   /*--------------------------------------------------------------*/
   /*          Find the local maxima and threshold it              */
   /*--------------------------------------------------------------*/
   for ( i = 0 ; i < fsz-1 ; i++)
     for ( j = 1 ; j < fsz-1 ; j++) 
        if ((M[i][j]-M[i+1][j]) > THRESH || (M[i][j]-M[i][j+1])>THRESH ||
            (M[i][j]-M[i+1][j-1]) > THRESH || (M[i][j]-M[i+1][j+1]) > THRESH )
                edge[i][j] = BLACK ;
           else edge[i][j] = WHITE ;

   img_dsp1(edge,fsz,v_loc,h_loc) ;
   if (save_img == 'y') {
               sprintf(outname,"%s.eg%4.1f",fname,THRESH);
               write_img(edge,outname,fsz) ;
   }

}