/*---------------------------------------------------------------------------*/
/* Program:  histo.c                                                         */
/*                                                                           */
/* Purpose:  This program calculates a histogram of intensity values         */
/*           for an image in KUIM format.                                    */
/*                                                                           */
/* Author:   John Gauch, Bona Nasution                                       */
/*                                                                           */
/* Date:     March 30, 1994   - Original program.                            */
/*           August 3, 1994   - Version 2.                                   */
/*           November 1, 1994 - Version 3.                                   */
/*           July 16, 1997    - Output 2D image of histogram.                */
/*                                                                           */
/* Note:     Copyright (C) The University of Kansas, 1994                    */
/*---------------------------------------------------------------------------*/
#include <IM.h>

/*---------------------------------------------------------------------------*/
/* 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;
   BYTE_TYPE **Data2;
   int Histo[256];
   int Cumulative[256];
   int PixType, Xdim, Ydim, Zdim, DimCnt;

   /* Program variables */
   int Debug = FALSE;
   int Ascii = FALSE;
   int i = 0;
   int Count;
   float Min, Max;
   float Scale;
   int MAXBIN = 40;
   float Ave, Ave2, Var;
   FILE *Fd;
   int x, y;

   /* Interpret program options */
   fprintf(stderr, "HISTO Program - KUIM Version 3.0\n\n");
   while ((++i < argc) && (argv[i][0] == '-'))
      switch (argv[i][1])
      {
      case 'd':
	 Debug = TRUE;
	 break;
      case 'a':
	 Ascii = TRUE;
	 break;
      default:
	 Error("Invalid option encountered");
	 break;
      }

   /* Check number of file names */
   if (argc - i < 1)
   {
      fprintf(stderr, "Usage: histo [options] infile [outfile]\n");
      fprintf(stderr, "       [-d]    Print debugging information\n");
      fprintf(stderr, "       [-a]    Output ascii for xgraph\n");
      exit(1);
   }

   /* Get image file names from argument list */
   Name1[0] = Name2[0] = '\0';
   if (sscanf(argv[i++], "%s", Name1) == 0)
      Error("Could not get input file name");
   if ((i < argc) && (sscanf(argv[i++], "%s", Name2) == 0))
      Error("Could not get output file name");
   if (Debug == TRUE)
      printf("%s %s\n", Name1, Name2);

   /* Read input image */
   Image1 = im_open(Name1, &PixType, &Xdim, &Ydim, &Zdim, &DimCnt);
   Data1 = (FLOAT_TYPE *) im_alloc1D(Image1, FLOAT);
   im_read(Image1, FLOAT, (char *) Data1);

   /* Get image range */
   Ave = Ave2 = Var = 0.0;
   Min = Max = Data1[0];
   Count = Xdim * Ydim * Zdim;
   for (i = 0; i < Count; i++)
   {
      if (Data1[i] < Min)
	 Min = Data1[i];
      else if (Data1[i] > Max)
	 Max = Data1[i];
      Ave += Data1[i];
      Ave2 += Data1[i] * Data1[i];
   }
   Ave /= Count;
   Ave2 /= Count;
   Var = Ave2 - Ave * Ave;
   if (Max == Min)
      Error("Pixels all have same the value");

   /* Create output image */
   if ((strcmp(Name2, "") != 0) || (Ascii == TRUE))
      MAXBIN = 256;

   /* Get intensity histogram */
   Scale = (MAXBIN - 1) / (Max - Min);
   for (i = 0; i < MAXBIN; i++)
      Histo[i] = Cumulative[i] = 0;
   for (i = 0; i < Count; i++)
      Histo[(int) ((Data1[i] - Min) * Scale)]++;

   /* Get cumulative histogram */
   Cumulative[0] = Histo[0];
   for (i = 1; i < MAXBIN; i++)
      Cumulative[i] = Cumulative[i - 1] + Histo[i];

   /* Write histogram to file */
   if ((strcmp(Name2, "") != 0) && (Ascii == FALSE))
   {
      /* Scale Histo between 0..200 */
      Min = Max = (float) Histo[0];
      for (i = 0; i < MAXBIN; i++)
      {
	 if (Min > Histo[i])
	    Min = (float) Histo[i];
	 else if (Max < Histo[i])
	    Max = (float) Histo[i];
      }
      for (i = 0; i < MAXBIN; i++)
	 if (Max > Min)
	    Histo[i] = (int) (((Histo[i] - Min) * 200) / (Max - Min));
	 else
	    Histo[i] = 0;

      /* Create output image */
      Image2 = im_create(Name2, BYTE, MAXBIN, MAXBIN, 1);
      Data2 = (BYTE_TYPE **) im_alloc2D(Image2, BYTE);
      for (y = 0; y < MAXBIN; y++)
	 for (x = 0; x < MAXBIN; x++)
	    if ((MAXBIN - y) < (Histo[x]))
	       Data2[y][x] = 200;
	    else
	       Data2[y][x] = 0;
      im_write(Image2, BYTE, &(Data2[0][0]));
   }

   /* Print histogram in xgraph format */
   else if ((strcmp(Name2, "") != 0) && (Ascii == TRUE))
   {
      Fd = fopen(Name2, "w");
      fprintf(Fd, "TitleText: %s\n", Name1);
      for (i = 0; i < MAXBIN; i++)
	 fprintf(Fd, "%d %d\n", i, Histo[i]);
      fclose(Fd);
   }

   /* Print histogram in xgraph format */
   else if (Ascii == TRUE)
   {
      printf("TitleText: %s\n", Name1);
      for (i = 0; i < MAXBIN; i++)
	 printf("%d %d\n", i, Histo[i]);
   }

   /* Print histogram and cumulative histogram */
   else
   {
      printf("   Intensity Range       Histogram      Cumulative\n");
      for (i = 0; i < MAXBIN; i++)
	 printf("[%6.2f .. %6.2f) %8d %6.2f %8d %6.2f\n",
		Min + i / Scale, Min + (i + 1) / Scale, Histo[i], Histo[i] * 100.0 / Count,
		Cumulative[i], Cumulative[i] * 100.0 / Count);

      printf("\nImage Statistics:\n");
      printf("Mean     = %7.2f\n", Ave);
      printf("Variance = %7.2f\n", Var);
      printf("StdDev   = %7.2f\n", (float) sqrt((double) Var));
   }
   return (0);
}