 /*
  * Khoros: $Id: lvspatial.c,v 1.2 1992/03/20 23:07:40 dkhoros Exp $
  */

#if !defined(lint) && !defined(SABER)
static char rcsid[] = "Khoros: $Id: lvspatial.c,v 1.2 1992/03/20 23:07:40 dkhoros Exp $";
#endif

 /*
  * $Log: lvspatial.c,v $
 * Revision 1.2  1992/03/20  23:07:40  dkhoros
 * VirtualPatch5
 *
  */ 

/*
 *----------------------------------------------------------------------
 *
 * Copyright 1991, University of New Mexico.  All rights reserved.
 * Permission to copy and modify this software and its documen-
 * tation only for internal use in your organization is hereby
 * granted, provided that this notice is retained thereon and
 * on all copies.  UNM makes no representations as too the sui-
 * tability and operability of this software for any purpose.
 * It is provided "as is" without express or implied warranty.
 * 
 * UNM DISCLAIMS ALL WARRANTIES WITH REGARD TO THIS SOFTWARE,
 * INCLUDING ALL IMPLIED WARRANTIES OF MERCHANTABILITY AND FIT-
 * NESS.  IN NO EVENT SHALL UNM BE LIABLE FOR ANY SPECIAL,
 * INDIRECT OR CONSEQUENTIAL DAMAGES OR ANY OTHER DAMAGES WHAT-
 * SOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER
 * IN AN ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS
 * ACTION, ARISING OUT OF OR IN CONNECTION WITH THE USE OR PER-
 * FORMANCE OF THIS SOFTWARE.
 * 
 * No other rights, including for example, the right to redis-
 * tribute this software and its documentation or the right to
 * prepare derivative works, are granted unless specifically
 * provided in a separate license agreement.
 *---------------------------------------------------------------------
 */

#include "unmcopyright.h"        /* Copyright 1991 by UNM */

/*>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>  <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
 >>>>
 >>>>         File Name: lvspatial.c
 >>>>
 >>>>      Program Name: vspatial
 >>>>
 >>>> Date Last Updated: Tue Mar  5 22:27:15 1991 
 >>>>
 >>>>          Routines: lvspatial - the library call for vspatial
 >>>>
 >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>   <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<*/


#include "vinclude.h"


/* -library_includes */
/* -library_includes_end */


/****************************************************************
*
* Routine Name: lvspatial - library call for vspatial
*
* Purpose:
*    
*    Perform  spatial  feature  extraction  of   an   image   via
*    statistical analysis of window pixels.
*    
*    
* Input:
*    
*    img            xvimage structure
*    
*    wsize          window width argument
*    
*    hsize          window height argument
*    
*    m              mean operator argument
*    
*    v              variance operator argument
*    
*    c              contrast operator argument
*    
*    s              angular 2nd moment argument
*    
*    e              entropy operator argument
*    
*    d              dispersion operator argument
*    
*    
* Output:
*    
*    img            holds the resulting feature of the image.  img  is
*                   used  for both the input xvimage structure and the
*                   output result xvimage structure.
*    
*    
*
* Written By: Charlie Gage
*    
*    
****************************************************************/


/* -library_def */
int
lvspatial(img, xsize, ysize, m, v, c, s, e, d)
struct xvimage *img;
int xsize, ysize, m, v, c, s, e, d;
/* -library_def_end */

/* -library_code */
{

    int row,                            /* points image row */ 
        col,                            /* points image col */ 
        x,                              /* pixel variable */
        xoff,                           /* horiz. offset in the img */
        yoff,                           /* vertic. offset in the img */
        xend,                           /* step in img up to here */
        yend,                           /* step in img up to here */
        nc,                             /* column size of image */
        nr,                             /* row size of image */
        n,                              /* window size --> (xsize*ysize) */
        sum,                            /* sum of pixel values in window */
        l,                              /* number of gray levels */
        nb,                             /* number of bands in image */
        bandsize,                       /* size of one image band */
        i, j, k, p;                     /* general counters */

    float  *result,                     /* pointer to resulting image char */
           z,                           /* temporary pixel value */
           mean,                        /* mean value of a window */
           variance,                    /* variance value of a window */
           contrast,                    /* contrast value of a window */
           ang_sec_mom,                 /* Angular Second Moment of a window */
           dispersion;                  /* dispersion value of a window */

    double entropy,                     /* entropy value of a window */
           log_base2,                   /* logarithm to the base 2 */
           log10_two;                   /* constant for log_base2 calculation */

    float  *mptr,                       /* pointer to mean feature band */
           *vptr,                       /* pointer to variance feature band */
           *cptr,                       /* pointer to contrast feature band */
           *sptr,                       /* pointer to ang_sec_mom feat. band*/
           *eptr,                       /* pointer to entropy feature band */
           *dptr;                       /* pointer to dispersion feature band */

    float *bptr[6];                     /* array of pointers to feature bands */
    unsigned int flag[6];               /* array of flags for feature bands */

    unsigned int hist[256];             /* holds histogram of window */

    unsigned char *ptr;                 /* pointer to the image data char */

    char  *malloc(),                    /* function definition */
          *program = "lvspatial";       /* contains the library name */


   /* check for odd windows */
    if((xsize & 1) == 0 || (ysize & 1) == 0){
        (void) fprintf(stderr,"lvspatial: Error window must be odd\n");
        return(0);
    }

    if (!propertype (program, img, VFF_TYP_1_BYTE, FALSE)) {
        (void) fprintf (stderr,
        "lvspatial: ERROR Input image must be of data type BYTE\n");
         return (0);
    }

    n = xsize * ysize;          /* size of window */
    nr = img->col_size;         /* number of rows in image */
    nc = img->row_size;         /* number of cols in image */

    bandsize = ((unsigned int)nc * nr);   /* image bandsize */

    xoff = xsize/2;             /* offset from center of window (x-dir) */
    yoff = ysize/2;             /* offset from center of window (y-dir) */

    xend = nc - xsize;          /* step in image (x-dir) up to here */
    yend = nr - ysize + 1;      /* step in image (y-dir) up to here */

    l = 256;                    /* number of gray levels */
    p = 0;                      /* general counter */
    j = 0;
    nb = 0;                     /* initialize number of bands to zero */
    log10_two = log10(2.0);     /* used to calculate log_base2 */
    sum = 0;
    ang_sec_mom = 0.0;          /* Angular Second Moment */
    contrast = 0.0;
    dispersion = 0.0;
    entropy = 0.0;
    variance = 0.0;

   /* Zero flags */
    bzero(flag, 6 * sizeof(unsigned int));

   /* Zero out histogram */
    bzero((char *) hist, 256 * sizeof(unsigned int));

   /* determine number of output bands and which bands are requested */
   if (m == 1) {            /* set mean flag */
      nb++;
      flag[0] = 1;
   }
   if (v == 1) {            /* set variance flag */
      nb++;
      flag[1] = 1;
   }
   if (c == 1) {            /* set contrast flag */
      nb++;
      flag[2] = 1;
   }
   if (s == 1) {            /* set ang_sec_mom flag */
      nb++;
      flag[3] = 1;
   }
   if (e == 1) {            /* set entropy flag */
      nb++;
      flag[4] = 1;
   }
   if (d == 1) {            /* set dispersion flag */
      nb++;
      flag[5] = 1;
   }

   /* Allocate space for resulting image -- of type FLOAT */ 
    result = (float *) malloc((unsigned int)nb*nc*nr*sizeof(float));
    if(result == NULL) {
        (void) fprintf(stderr,"lvspatial: insufficient memory available\n");
        return(0);
    }

   /* Assign input image data address to ptr */
    ptr = (unsigned char *)img->imagedata;

   /* Assign pointers to appropriate bands in output image (result) */
   for (i=0; i < 6; i++) {
       if (flag[i]) {
          bptr[i] = &result[j * bandsize];
          j++;
       }
   }
   
   /* Assign meaningful pointer names to appropriate bands */
   mptr = bptr[0];           /* pointer to mean band */
   vptr = bptr[1];           /* pointer to variance band */
   cptr = bptr[2];           /* pointer to contrast band */
   sptr = bptr[3];           /* pointer to ang_sec_mom band */
   eptr = bptr[4];           /* pointer to entropy band */
   dptr = bptr[5];           /* pointer to dispersion band */

   /* perform spatial textural analysis of an image */
        /* loop through rows -- (y-direction) */
    for (row = 0 ; row < yend; row++){

       /* setup the histogram for the current row */
        for(i = 0; i < ysize; i++){
            for(j = 0; j < xsize; j++){
                x = ptr[nc*(row + i) + j];
                sum += x;
                hist[x]++;
            }
        }

       /* calculate the mean value in a window */
        mean = (float) sum / n;

       /* perform selected feature extraction on window */
         if (m == 1) {      /* update center pixel using mean */
              mptr[nc * (row + yoff) + xoff + p] = mean;
         }
         if (v == 1) {      /* calculate variance and update center pixel */
              for (i = 0; i < l; i++) {
                  if (hist[i] != 0) {
                      z = (float) i;
                      variance += (z - mean)*(z - mean)*((float) hist[i]/n );
                  }
              }
              vptr[nc * (row + yoff) + xoff + p] = variance;
         }            
         if (c == 1) {      /* calculate contrast and update center pixel */
              for (i = 0; i < l; i++) {
                  if (hist[i] != 0) {
                      z = (float) i;
                      contrast += z * z * ((float) hist[i]/n );
                  }
              }
              cptr[nc * (row + yoff) + xoff + p] = contrast;
         }            
         if (s == 1) {      /* calculate ASM and update center pixel */
              for (i = 0; i < l; i++) {
                  if (hist[i] != 0) {
                    ang_sec_mom += ((float) hist[i]/n) * ((float) hist[i]/n);
                  }
              }
              sptr[nc * (row + yoff) + xoff + p] = ang_sec_mom;
         }            
         if (e == 1) {      /* calculate entropy and update center pixel */
              for (i = 0; i < l; i++) {
                  if (hist[i] != 0) {
                      log_base2 = log10((double) hist[i]/n) / log10_two;
                      entropy += ((double) hist[i]/n) * log_base2;
                  }
              }
              eptr[nc * (row + yoff) + xoff + p] = (float) -entropy;
         }            
         if (d == 1) {      /* calculate dispersion and update center pixel */
              for (i = 0; i < l; i++) {
                  if (hist[i] != 0) {
                      z = (float) i;
                      dispersion += fabs(z - mean) * ((float) hist[i]/n);
                  }
              }
              dptr[nc * (row + yoff) + xoff + p] = dispersion;
         }    
       /* reset variables to zero for next window */
       variance = 0.0;
       contrast = 0.0;
       ang_sec_mom = 0.0;
       entropy = 0.0;
       dispersion = 0.0;
       p++;

       /* step through the current row of image */
                /* loop through the cols -- (x-direction) */
        for (col = 0 ; col < xend ; col++){

            for(k = 0; k < ysize; k++){

               /* delete the lefmost column of the previous window*/
                x =  ptr[col + nc*(row + k)];
                sum = sum - x;
                hist[x] = hist[x] - 1;

               /* add the rightmost column of the current window*/
                x = ptr[xsize + nc*(row + k) + col];
                sum = sum + x;
                hist[x] = hist[x] + 1;
            }

       /* calculate the mean value in a window */
        mean = (float) sum / n;

       /* perform selected feature extraction on window */
         if (m == 1) {      /* update center pixel using mean */
              mptr[nc * (row + yoff) + xoff + p] = mean;
         }
         if (v == 1) {      /* calculate variance and update center pixel */
              for (i = 0; i < l; i++) {
                  if (hist[i] != 0) {
                      z = (float) i;
                      variance += (z - mean)*(z - mean)*((float) hist[i]/n);
                  }
              }
              vptr[nc * (row + yoff) + xoff + p] = variance;
         }            
         if (c == 1) {      /* calculate contrast and update center pixel */
              for (i = 0; i < l; i++) {
                  if (hist[i] != 0) {
                      z = (float) i;
                      contrast += z * z * ((float) hist[i]/n);
                  }
              }
              cptr[nc * (row + yoff) + xoff + p] = contrast;
         }            
         if (s == 1) {      /* calculate ASM and update center pixel */
              for (i = 0; i < l; i++) {
                  if (hist[i] != 0) {
                    ang_sec_mom += ((float) hist[i]/n) * ((float) hist[i]/n);
                  }
              }
              sptr[nc * (row + yoff) + xoff + p] = ang_sec_mom;
         }            
         if (e == 1) {      /* calculate entropy and update center pixel */
              for (i = 0; i < l; i++) {
                  if (hist[i] != 0) {
                      log_base2 = log10((double) hist[i]/n) / log10_two;
                      entropy += ((double) hist[i]/n) * log_base2;
                  }
              }
              eptr[nc * (row + yoff) + xoff + p] = (float) -entropy;
         }            
         if (d == 1) {      /* calculate dispersion and update center pixel */
              for (i = 0; i < l; i++) {
                  if (hist[i] != 0) {
                      z = (float) i;
                      dispersion += fabs(z - mean) * ((float) hist[i]/n);
                  }
              }
              dptr[nc * (row + yoff) + xoff + p] = dispersion;
         }    

        /* reset variables to zero for next window */
        ang_sec_mom = 0.0;
        contrast = 0.0;
        dispersion = 0.0;
        entropy = 0.0;
        variance = 0.0;
        p++;

      } /* End of col loop */
       /* Reset histogram and variables for next pass */
        bzero((char *) hist, 256 * sizeof(unsigned int));
        p = 0;
        sum = 0;
    }   /* End of row loop */

    (void) free (img->imagedata);
    img->data_storage_type = VFF_TYP_FLOAT;
    img->imagedata = (char *) result;
    img->num_data_bands = nb;

    return(1);
}

/* -library_code_end */
