/******************************************************************************/
/*                                                                            */
/* Lim & Lee algorithm to perform multispectral edge detection                */
/*                                                                            */
/******************************************************************************/

#include "VisXV4.h"               /* VisX structure include file     */
#include "Vutil.h"                /* VisX utility header files       */

VXparam_t par[] =                 /* command line structure          */
{
{ "if=",    0,   "input file"},
{ "of=",    0,   "output file"},
{ "g=",     0,   "gradient threshold"},
{  0,       0,   0}   /* list termination */
};
#define  IVAL   par[0].val
#define  OVAL   par[1].val
#define  GVAL   par[2].val

VisXfile_t *VXin,                  /* input file structure               */
           *VXout ;                /* output file structure              */
VisXelem_t *VXlist,                /* VisionX data structure             */
           *VXptr ;                /* VisionX pixel data pointer         */

/* global variables */
VisXimage_t r, g, b,               /* RGB image structure                */
            r_grad, g_grad,        /* image structures of gradients of   */
            b_grad ;               /* individual spectra                 */

/* Function declaration */           
void sobelr (int, int) ;
void sobelg (int, int) ;
void sobelb (int, int) ;           
           
main(argc, argv)
int argc ;
char *argv[] ; {
   int i, j,                       /* loop counters                      */
       gt ;                        /* gradient threshold                 */
   VisXimage_t in,                 /* input image structure              */
               grad ;              /* gradient image structure           */
 
   VXparse(&argc, &argv, par) ;    /* parse the command line             */
   VXin  = VXopen(IVAL, 0) ;       /* open input file                    */
   VXout = VXopen(OVAL, 1) ;       /* open the output file               */
   VXlist = VXread(VXin) ;         /* read file                          */
   gt = (GVAL ? atoi(GVAL) : 50) ; /* read gt, default is gt=50          */

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

1) copy image data into array "in" which has 3 times the WIDTH to 
   accommodate R, G, B data.
2) create arrays "r_grad", "g_grad", "b_grad" with dimentions 
   WIDTH * HEIGHT (float type)
3) create arrays "r", "g", "b" with dimensions WIDTH * HEIGHT with a border 
   of zeros (float type because don't know how to embed and change type)
4) create output array "grad" with dimensions WIDTH * HEIGHT

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

   if (VXNIL != (VXptr = VXfind(VXlist, VX_PBYTE))) {  /* find image */
      VXsetimage(&in, VXptr, VXin) ;   /* initialize input structure */

      VXmakeimage(&r_grad, VX_PFLOAT, in.bbx, 1) ; // float for more accuracy
      VXmakeimage(&g_grad, VX_PFLOAT, in.bbx, 1) ;
      VXmakeimage(&b_grad, VX_PFLOAT, in.bbx, 1) ;

      VXembedimage(&r, &r_grad, 1, 1, 1, 1) ;  // image structure with border of zeros
      VXembedimage(&g, &r_grad, 1, 1, 1, 1) ;  // image structure with border of zeros
      VXembedimage(&b, &r_grad, 1, 1, 1, 1) ;  // image structure with border of zeros

      VXmakeimage(&grad, VX_PBYTE, in.bbx, 1) ;
   
      /* split data from "in" into the 3 arrays "r", "g", "b" */
      for (i = in.ylo ; i <= in.yhi ; i++) {
         for (j = in.xlo ; j <= in.xhi ; j++) {
            if (j%3 == 0)
               r.f[i][j/3] = in.u[i][j] ;
            if (j%3 == 1)
               g.f[i][j/3] = in.u[i][j] ;
            if (j%3 == 2)
               b.f[i][j/3] = in.u[i][j] ;
         }
      }

      /* loop to calculate gradients in each spectra */  
      for (i = r_grad.ylo ; i <= r_grad.yhi ; i++) {
         for (j = r_grad.xlo ; j <= r_grad.xhi ; j++) {
            sobelr (i, j) ;
            sobelg (i, j) ;
            sobelb (i, j) ;
         }
      }

      /* loop to calculate average gradient and mark out edges */  
      for (i = r_grad.ylo ; i <= r_grad.yhi ; i++) {
         for (j = r_grad.xlo ; j <= r_grad.xhi ; j++) {
            grad.u[i][j] = ((r_grad.f[i][j] + g_grad.f[i][j] + b_grad.f[i][j]) / 3) ;
            if (grad.u[i][j] > gt)
               grad.u[i][j] = 0 ; // black
            else
               grad.u[i][j] = 255 ; // white
         }
      }
   
      VXwrite(VXout, grad.list) ;        /* write data  */  
   }
   else {
      fprintf(stderr, "no byte image data in input file\n") ;
      exit(-1) ;
   }       
   VXclose(VXin);                  /* close files  */
   VXclose(VXout);
   exit(0);
}

void sobelr (int y, int x) {
   float sx, sy ;
   sx = r.f[y+1][x-1]*(-1) + r.f[y][x-1]*(-2) + r.f[y-1][x-1]*(-1) +
        r.f[y+1][x+1] + r.f[y][x+1]*2 + r.f[y-1][x+1] ;
   sy = r.f[y+1][x-1] + r.f[y+1][x]*2 + r.f[y+1][x+1] +
        r.f[y-1][x-1]*(-1) + r.f[y-1][x]*(-2) + r.f[y-1][x+1]*(-1) ;
   r_grad.f[y][x] = sqrt (sx * sx + sy * sy) ;
}

void sobelg (int y, int x) {
   float sx, sy ;
   sx = g.f[y+1][x-1]*(-1) + g.f[y][x-1]*(-2) + g.f[y-1][x-1]*(-1) +
        g.f[y+1][x+1] + g.f[y][x+1]*2 + g.f[y-1][x+1] ;
   sy = g.f[y+1][x-1] + g.f[y+1][x]*2 + g.f[y+1][x+1] +
        g.f[y-1][x-1]*(-1) + g.f[y-1][x]*(-2) + g.f[y-1][x+1]*(-1) ;
   g_grad.f[y][x] = sqrt (sx * sx + sy * sy) ;
}

void sobelb (int y, int x) {
   float sx, sy ;
   sx = b.f[y+1][x-1]*(-1) + b.f[y][x-1]*(-2) + b.f[y-1][x-1]*(-1) +
        b.f[y+1][x+1] + b.f[y][x+1]*2 + b.f[y-1][x+1] ;
   sy = b.f[y+1][x-1] + b.f[y+1][x]*2 + b.f[y+1][x+1] +
        b.f[y-1][x-1]*(-1) + b.f[y-1][x]*(-2) + b.f[y-1][x+1]*(-1) ;
   b_grad.f[y][x] = sqrt (sx * sx + sy * sy) ;
}