/*---------------------------------------------------------------------------*/
/* Program:  greymap.c                                                       */
/*                                                                           */
/* Purpose:  This program performs intensity scaling and windowing.          */
/*                                                                           */
/* Author:   John Gauch                                                      */
/*                                                                           */
/* Date:     March 28, 1994                                                  */
/*                                                                           */
/* Note:     Copyright (C) The University of Kansas, 1994                    */
/*---------------------------------------------------------------------------*/
#include <IM.h>
#define MAXBIN 1024

/*---------------------------------------------------------------------------*/
/* 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;
   FLOAT_TYPE *Data1;
   int PixType, Xdim, Ydim, Zdim, DimCnt;

   /* Program variables */
   int Debug = FALSE;
   int Count;
   int i = 0;
   float Min, Max;
   float Smin, Smax;
   float Wmin, Wmax;
   int ScaleFlag = FALSE;
   int WindowFlag = FALSE;
   float Scale;
   int HistoFlag = FALSE;
   int TrimFlag = FALSE;
   float Map[MAXBIN];

   /* Interpret program options */
   printf("GREYMAP Program - KUIM Version 1.0\n\n");
   while ((++i < argc) && (argv[i][0] == '-'))
      switch (argv[i][1])
      {
      case 's':
	 ScaleFlag = TRUE;
	 if (sscanf(argv[++i], "%f", &Smin) == 0)
	    Error("Could not get Smin argument");
	 if (sscanf(argv[++i], "%f", &Smax) == 0)
	    Error("Could not get Smax argument");
	 break;
      case 'w':
	 WindowFlag = TRUE;
	 if (sscanf(argv[++i], "%f", &Wmin) == 0)
	    Error("Could not get Wmin argument");
	 if (sscanf(argv[++i], "%f", &Wmax) == 0)
	    Error("Could not get Wmax argument");
	 break;
      case 'h':
	 HistoFlag = TRUE;
	 break;
      case 't':
	 TrimFlag = TRUE;
	 break;
      case 'd':
	 Debug = TRUE;
	 break;
      default:
	 Error("Invalid option encountered");
	 break;
      }

   /* Check number of file names */
   if (argc - i != 2)
   {
      fprintf(stderr, "Usage: greymap [options] infile outfile\n");
      fprintf(stderr, "       [-d]      Print debugging information\n");
      fprintf(stderr, "       [-s # #]  Scale pixels to specified range\n");
      fprintf(stderr, "       [-w # #]  Window pixels to specified range\n");
      fprintf(stderr, "       [-h]      Perform histogram equalization\n");
      fprintf(stderr, "       [-t]      Trim top and bottom 1 percent of range\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);
   Count = Xdim * Ydim * Zdim;
   if (PixType == COLOR || PixType == PSEUDO || PixType == COMPLEX)
      Error("Can not greymap color or complex images");
   Data1 = (FLOAT_TYPE *) im_alloc1D(Image1, FLOAT);
   im_read(Image1, FLOAT, (char *) Data1);

   /* Window input image */
   if (WindowFlag == TRUE)
   {
      for (i = 0; i < Count; i++)
	 if (Data1[i] < Wmin)
	    Data1[i] = Wmin;
	 else if (Data1[i] > Wmax)
	    Data1[i] = Wmax;
   }

   /* Scale input image */
   if (ScaleFlag == TRUE)
   {
      Min = Max = Data1[0];
      for (i = 0; i < Count; i++)
	 if (Data1[i] < Min)
	    Min = Data1[i];
	 else if (Data1[i] > Max)
	    Max = Data1[i];
      if (Max > Min)
	 Scale = (Smax - Smin) / (Max - Min);
      else
	 Scale = 0;
      for (i = 0; i < Count; i++)
	 Data1[i] = (Data1[i] - Min) * Scale + Smin;
   }

   /* Perform histogram equalization */
   if (HistoFlag == TRUE)
   {
      /* Get image range */
      Min = Max = Data1[0];
      for (i = 0; i < Count; i++)
	 if (Data1[i] < Min)
	    Min = Data1[i];
	 else if (Data1[i] > Max)
	    Max = Data1[i];
      if (Max > Min)
	 Scale = (MAXBIN - 1) / (Max - Min);
      else
	 Scale = 0;

      /* Get cumulative histogram */
      for (i = 0; i < MAXBIN; i++)
	 Map[i] = 0;
      for (i = 0; i < Count; i++)
	 Map[(int) ((Data1[i] - Min) * Scale)]++;
      for (i = 1; i < MAXBIN; i++)
	 Map[i] = Map[i] + Map[i - 1];

      /* Remap input image */
      for (i = 0; i < MAXBIN; i++)
	 Map[i] = Min + (Max - Min) * Map[i] / Count;
      for (i = 0; i < Count; i++)
	 Data1[i] = Map[(int) ((Data1[i] - Min) * Scale)];
   }

   /* Perform histogram trimming */
   if (TrimFlag == TRUE)
   {
      /* Get image range */
      Min = Max = Data1[0];
      for (i = 0; i < Count; i++)
	 if (Data1[i] < Min)
	    Min = Data1[i];
	 else if (Data1[i] > Max)
	    Max = Data1[i];
      if (Max > Min)
	 Scale = (MAXBIN - 1) / (Max - Min);
      else
	 Scale = 0;

      /* Get cumulative histogram */
      for (i = 0; i < MAXBIN; i++)
	 Map[i] = 0;
      for (i = 0; i < Count; i++)
	 Map[(int) ((Data1[i] - Min) * Scale)]++;
      for (i = 1; i < MAXBIN; i++)
	 Map[i] = Map[i] + Map[i - 1];

      /* Find trim values */
      for (i = 0; Map[i] < Count / 100; i++);
      Wmin = i / Scale + Min;
      printf("low index = %d low trim = %f\n", i, Wmin);
      for (i = MAXBIN - 1; Map[i] > Count * 99 / 100; i--);
      Wmax = i / Scale + Min;
      printf("high index = %d high trim = %f\n", i, Wmax);

      /* Window image to trim histogram */
      for (i = 0; i < Count; i++)
	 if (Data1[i] < Wmin)
	    Data1[i] = Wmin;
	 else if (Data1[i] > Wmax)
	    Data1[i] = Wmax;
   }

   /* Write information to output image */
   Image2 = im_create(Name2, PixType, Xdim, Ydim, Zdim);
   im_write(Image2, FLOAT, (char *) Data1);
   im_free1D((char *) Data1);
   return (0);
}