/*---------------------------------------------------------------------------*/
/* Program:  wrr.c                                                           */
/*                                                                           */
/* Purpose:  This program performs local contrast enhancement by computing   */
/*           the weighted region ranking (WRR) for each pixel in the image.  */
/*           Two parameters control this process:                            */
/*                                                                           */
/*           1)  The size of the contextual region about each pixel is       */
/*           specified by the standard deviation of a Gaussian weighting     */
/*           function.  The values of this function control the amount       */
/*           each pixel in the contextual region contributes to the pseudo-  */
/*           histogram and hence the rank of each pixel.  Like WAHE this     */
/*           provides a smooth transition from pixels in the contextual      */
/*           region to those outside.                                        */
/*                                                                           */
/*           2)  The amount of contrast enhancement of each pixel is         */
/*           specified by the standard deviation of a Gaussian smoothing     */
/*           function.  By blurring the pseudo-histogram the first and       */
/*           second derivatives of the cumulative histogram are reduced.     */
/*           Like CLAHE this reduces the contrast enhancement at each        */
/*           pixel.  Unlike CLAHE this method reduces the rate of change     */
/*           of contrast enhancement.                                        */
/*                                                                           */
/* Author:   John Gauch                                                      */
/*                                                                           */
/* Date:     April 22, 1994                                                  */
/*                                                                           */
/* Note:     Copyright (C) The University of Kansas, 1994                    */
/*---------------------------------------------------------------------------*/
#include <IM.h>

/* Global variables */
float BlurStdDev = 1.0;
float WeightStdDev = 8.0;
int Radius = 12;
FLOAT_TYPE **Data1;
FLOAT_TYPE **Data2;
int Xdim, Ydim, Zdim;
int Debug = FALSE;
double exp();
void WRR();

/*---------------------------------------------------------------------------*/
/* Purpose:  This is the main program.                                       */
/*---------------------------------------------------------------------------*/
int main(int argc, char *argv[])
{

   /* Image variables */
   char Name1[50];
   char Name2[50];
   IM_TYPE *Image1;
   IM_TYPE *Image2;
   int PixType, DimCnt;

   /* Program variables */
   int i = 0;

   /* Interpret program options */
   printf("WRR (Weighted Region Ranking) Program - KUIM Version 1.0\n\n");
   while ((++i < argc) && (argv[i][0] == '-'))
      switch (argv[i][1])
      {
      case 'b':
	 if (sscanf(argv[++i], "%f", &BlurStdDev) == 0)
	    Error("Could not get BlurStdDev argument");
	 if (BlurStdDev < 0.1)
	    BlurStdDev = 0.1;
	 break;
      case 's':
	 if (sscanf(argv[++i], "%f", &WeightStdDev) == 0)
	    Error("Could not get WeightStdDev argument");
	 if (WeightStdDev < 0.1)
	    WeightStdDev = 0.1;
	 break;
      case 'r':
	 if (sscanf(argv[++i], "%d", &Radius) == 0)
	    Error("Could not get Radius argument");
	 if (Radius < 3)
	    Radius = 3;
	 break;
      case 'd':
	 Debug = TRUE;
	 break;
      default:
	 Error("Invalid option encountered");
	 break;
      }

   /* Check number of file names */
   if (argc - i != 2)
   {
      fprintf(stderr, "Usage: wrr [options] infile outfile\n");
      fprintf(stderr, "       [-d]    Print debugging information\n");
      fprintf(stderr, "       [-s #]  StdDev of weighting Gaussian (8)\n");
      fprintf(stderr, "       [-r #]  Radius of weighting region (12)\n");
      fprintf(stderr, "       [-b #]  StdDev of blurring Gaussian (1)\n");
      exit(1);
   }

   /* Get image file names from argument list */
   if (sscanf(argv[i++], "%s", Name1) == 0)
      Error("Could not get input file name");
   if (sscanf(argv[i++], "%s", Name2) == 0)
      Error("Could not get output file name");

   /* Read input image */
   Image1 = im_open(Name1, &PixType, &Xdim, &Ydim, &Zdim, &DimCnt);
   if (DimCnt != 2)
      Error("Can not process 1D or 3D images");
   Data1 = (FLOAT_TYPE **) im_alloc2D(Image1, FLOAT);
   im_read(Image1, FLOAT, (char *) &(Data1[0][0]));

   /* Create output image */
   Image2 = im_create(Name2, FLOAT, Xdim, Ydim, Zdim);
   Data2 = (FLOAT_TYPE **) im_alloc2D(Image2, FLOAT);

   /* Calculate weighted region ranking */
   WRR();

   /* Write information to output image */
   im_write(Image2, FLOAT, (char *) &(Data2[0][0]));
   im_free2D((char **) Data1);
   im_free2D((char **) Data2);
   return (0);
}

/*---------------------------------------------------------------------------*/
/* Purpose:  This routine computes the WRR for the image.                    */
/*---------------------------------------------------------------------------*/
void WRR()
{

   int x, y, dx, dy;
   float Variance;
   IM_TYPE *Image3;
   FLOAT_TYPE **Weight;
   float Temp, Area, Below, Center, Blur;

   /* Allocate weight image */
   Image3 = im_create("/dev/null", FLOAT, 2 * Radius + 1, 2 * Radius + 1, 1);
   Weight = (FLOAT_TYPE **) im_alloc2D(Image3, FLOAT);

   /* Initialize weight table */
   Variance = 2 * WeightStdDev * WeightStdDev;
   for (dy = -Radius; dy <= Radius; dy++)
      for (dx = -Radius; dx <= Radius; dx++)
	 Weight[dy + Radius][dx + Radius] = 
            (float) exp(-(dx * dx + dy * dy) / Variance);

   /* Handle case without histogram blurring */
   if (BlurStdDev <= 1)
   {
      for (y = 0; y < Ydim; y++)
	 for (x = 0; x < Xdim; x++)
	 {
	    /* Consider neighborhood radius R about point */
            Area = 0.0;
	    Below = 0.0;
	    Center = Data1[y][x];

	    for (dy = -Radius; dy <= Radius; dy++)
	       for (dx = -Radius; dx <= Radius; dx++)
	       {
		  if ((y + dy >= 0) && (y + dy < Ydim - 1) &&
		      (x + dx >= 0) && (x + dx < Xdim - 1))
                     {
                     Temp = Weight[dy + Radius][dx + Radius];
                     Area += Temp;
		     if (Data1[y + dy][x + dx] < Center)
			Below += Temp;
                     }
	       }
	    Data2[y][x] = 255 * Below / Area;
	 }
   }

   /* Handle slower case with histogram blurring */
   else
   {
      for (y = 0; y < Ydim; y++)
	 for (x = 0; x < Xdim; x++)
	 {
	    /* Consider neighborhood radius R about point */
	    Area = 0.0;
	    Below = 0.0;
	    Center = Data1[y][x];

	    for (dy = -Radius; dy <= Radius; dy++)
	       for (dx = -Radius; dx <= Radius; dx++)
	       {
		  if ((y + dy >= 0) && (y + dy < Ydim - 1) &&
		      (x + dx >= 0) && (x + dx < Xdim - 1))
		  {
		     Blur = (Center - Data1[y + dy][x + dx]) / BlurStdDev;
		     if (Blur < 0)
			Blur = 0;
		     else if (Blur > 1)
			Blur = 1;
                     Temp = Weight[dy + Radius][dx + Radius];
                     Area += Temp;
		     Below += Temp * Blur;
		  }
	       }
	    Data2[y][x] = 255 * Below / Area;
	 }
   }
}