 /*
  * Khoros: $Id: lvhxray.c,v 1.1 1991/05/10 15:41:54 khoros Exp $
  */

#if !defined(lint) && !defined(SABER)
static char rcsid[] = "Khoros: $Id: lvhxray.c,v 1.1 1991/05/10 15:41:54 khoros Exp $";
#endif

 /*
  * $Log: lvhxray.c,v $
 * Revision 1.1  1991/05/10  15:41:54  khoros
 * Initial revision
 *
  */ 

/*
 *----------------------------------------------------------------------
 *
 * 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: lvhxray.c
 >>>>
 >>>>      Program Name: vhxray
 >>>>
 >>>> Date Last Updated: Tue Mar  5 22:33:27 1991 
 >>>>
 >>>>          Routines: lvhxray - the library call for vhxray
 >>>>
 >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>   <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<*/


#include "vinclude.h"


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


/****************************************************************
*
* Routine Name: lvhxray - library call for vhxray
*
* Purpose:
*    
*    Enhance images using histogram stretch or equalization  with
*    overlapping windows.
*    
*    
* Input:
*    
*         image - struct xvimage -- the input image
*                                   and the output image
*         type -- int -- 1 = do histogram equalization
*                        2 = do histogram stretch
*         xsize -- int -- window width size
*         ysize -- int -- window heigth size
*    
*    
* Output:
*    
*         image - struct xvimage -- after the enhancement
*    
*         image is used for both the input xvimage structure  and  the
*         output  result  xvimage  structure.   This  is  done to save
*         space, but you must be careful not  to  overwrite  important
*         data.
*    
*         lvhxray returns 1 upon success and 0 on failure.
*    
*    
*
* Written By: Tom Sauer, Charlie Gage
*    
*    
****************************************************************/


/* -library_def */
int
lvhxray(image, type, xsize, ysize)

struct xvimage *image;
int type, xsize, ysize;
/* -library_def_end */

/* -library_code */
{

    int row,                            /* points image row */ 
        col,                            /* points image col */ 
        g1,                             /* auxiliar 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 */
        i, j, k, m;                     /* general counters */

    unsigned char med;                  /* this is the median value */

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

    unsigned char *ptr,                 /* pointer to the image data char */
                  *result;              /* pointer to resulting image char */

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


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

      /* can only have one data band */

    if (! (proper_num_images(program,image,1,FALSE))) {
         (void) fprintf(stderr,
                 "%s: ERROR All inputs must contain only one image\n"
                 , program);
          return(0);
     }

      /* can only have one data band */

     if (! (proper_num_bands(program,image,1,FALSE))) {
          (void) fprintf(stderr,
                  "%s: ERROR All inputs must contain only one data band\n"
                  , program);
           return(0);
      }

      /* can only have data type of BYTE */

     if (! (propertype(program,image,VFF_TYP_1_BYTE,FALSE))) {
          (void) fprintf(stderr,
                  "%s: ERROR input must have data type of BYTE\n"
                  , program);
           return(0);
      }

      /* must have map enable set to OPTIONAL */

     if (! (proper_map_enable(program,image,VFF_MAP_OPTIONAL,FALSE))) {
          (void) fprintf(stderr,
                  "%s: ERROR input must have map enable type OPTIONAL\n"
                  ,program);
           return(0);
      }


    nr = image->col_size;
    nc = image->row_size;

         /* compute the border size of pixels not updated */
    xoff = xsize/2;
    yoff = ysize/2;

        /* compute where the window should stop on the right and bottom */
    xend = nc - xsize;
    yend = nr - ysize + 1;
    m = 0; 

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

   /* Allocate space for resulting image */ 
    result = (unsigned char *)malloc((unsigned int)nc*nr*sizeof(unsigned char));
    if(result == NULL) {
        (void) fprintf(stderr,"%s: insufficient memory available\n", program);
        return(0);
    }


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

      /* scan each row of the image, stepping column by column */
      /* with the window */

    for (row = 0 ; row < yend; row++){

       /* setup the histogram for the current window */

        for(i = 0; i < ysize; i++){
            for(j = 0; j < xsize; j++){
                g1 = ptr[nc*(row + i) + j];
                hist[g1]++;
            }
        }

          /* compute the equalization or stretch of the window */

        switch(type) {
          case 0:
              (void) hist_equalization(hist, hist2);
              break;
          case 1:
              (void) hist_stretch(hist, hist2, xsize, ysize);
              break;
          default:
             fprintf(stderr,"%s: invalid histogram enhancement type\n",
                     program);
             return(0);
             break;
        }

         /* get the middle pixel value and assign that image location
            with the middle pixel value mapped through the histogram */

        med = image->imagedata[nc *(row + yoff) + xoff + m];
        result[nc *(row + yoff) + xoff + m++] = (unsigned char) hist2[med];

          

           /* step across the current row with the window */

        for (col = 0 ; col < xend ; col++){

              /* adjust the histogram for the next window */

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

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


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

            }
            
               /* compute the equalization or stretch of the window */

            switch(type) {
              case 0:
                  (void) hist_equalization(hist, hist2);
                  break;
              case 1:
                  (void) hist_stretch(hist, hist2, xsize, ysize);
                  break;
              default:
                 fprintf(stderr,"%s: invalid histogram enhancement type\n",
                         program);
                 return(0);
                 break;
            }

             /* get the middle pixel value and assign that image location
                with the middle pixel value mapped through the histogram */

            med = image->imagedata[nc *(row + yoff) + xoff + m];
            result[nc *(row + yoff) + xoff + m++] = (unsigned char) hist2[med];

        }

       /* Reset histogram and variables */

        bzero((char *) hist, 256 * sizeof(unsigned int));
        m = 0; 
    }

      /* assign the new image to the old image pointer */
    image->imagedata = (char *) result;

    return(1);
}

/**************************************************************
*
* MODULE NAME: hist_stretch
*
*     PURPOSE: Perform a histogram stretch on the histogram
*               passed in. This routine will perform the
*               maximum stretch possible.
*               This is the same stretch algorithm used be vhstr
*               written by Ramiro Jordan and Marcelo Teran.
*
*       INPUT: hist -- a 256 histogram 
*
*      OUTPUT: hist -- a 256 histogram
*
*
* CALLED FROM: lvhxray
*
* ROUTINES CALLED: 
*
**************************************************************/

int
hist_stretch(hist, hist2, xsize, ysize)

unsigned int *hist;
unsigned int *hist2;
int xsize, ysize;

{
    int  i, lo, hi, klo, khi, thresh;

    /* Find the histogram boundaries by locating the 1 percent levels */
    thresh = xsize * ysize * 0.01;
    lo  =  hi = -1;
    klo = khi = 0;
    for (i = 0; i < 256; i++)
      {
        klo += hist[i]; khi += hist[255-i];
        if (klo > thresh && lo == -1) lo = i;
        if (khi > thresh && hi == -1) hi = 255-i;
        if (lo != -1 && hi != -1) break;
      }

    /* Check to see if we have a reasonable contrast span. If not, re-do
       the bounds with a zero threshold. */
    if (hi == lo)
      {
        thresh = 0;
        lo  =  hi = -1;
        klo = khi = 0;
        for (i = 0; i < 256; i++)
          {
            klo += hist[i]; khi += hist[255-i];
            if (klo > thresh && lo == -1) lo = i;
            if (khi > thresh && hi == -1) hi = 255-i;
            if (lo != -1 && hi != -1) break;
          }
      } 

    /* If still have a zero span bound, quit */
    if (hi == lo) {
       bcopy(hist, hist2, 256*sizeof(int));
       return(1);
      }

    /* Re map the histogram */
    for (i = 0; i < 256; i++)
    {
        if (i < lo)
           hist2[i] = 0;
        else if (i > hi)
           hist2[i] = 255;
        else
           hist2[i] = 255*(i-lo)/(hi-lo);
    }

    return(1);
}


/**************************************************************
*
* MODULE NAME: hist_equalization
*
*     PURPOSE: Performs a histogram equalization on
*              a histogram. This is the same equalization
*               algorithm used by vheq, written by
*               Ramiro Jordan and Marcelo Teran.
*
*       INPUT: hist -- a 256 histogram
*
*      OUTPUT: hist -- a 256 histogram
*
*
* CALLED FROM: lvhxray
*
* ROUTINES CALLED: 
*
**************************************************************/

int
hist_equalization(hist, hist2)

unsigned int hist[];
unsigned int hist2[];

{
    int  i;
    unsigned int  map[256], k;
    float f;

    /* Integrate the histogram to get the equalization map */

    k = 0;
    for (i = 0; i < 256; i++)
    {
        k += hist[i];
        map[i] = k ;
    }

    /* Compute scale factor */
    if (map[255] == 0.0)  {
       bcopy(hist, hist2, 256*sizeof(int));
       return(1);
    }

    f = 255.0/map[255];
    
    /* scale the cumulative distribution */
    for (i = 0; i < 256; i++) 
       hist2[i] = (unsigned int) ((float) map[i] * f);

    return(1);
}
/* -library_code_end */
