/*---------------------------------------------------------------------------*/
/* Program:  unsharp.c                                                       */
/*                                                                           */
/* Purpose:  This program performs unsharp masking using binomial            */
/*           smoothing to obtain the "unsharp" image for subtraction.        */
/*                                                                           */
/* Author:   John Gauch, Edu Metz                                            */
/*                                                                           */
/* Date:     April 20, 1994 - Original program.                              */
/*           July 23, 1996  - Color support added.                           */
/*           July 31, 1997  - Simplified color, vector, scalar code.         */
/*                                                                           */
/* Note:     Copyright (C) The University of Kansas, 1994-1996               */
/*---------------------------------------------------------------------------*/
#include <IM.h>

/* Program variables */
int Debug = FALSE;
int Sharp = 1;
int Smooth = 5;
float Fraction = 0.8;

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

/*---------------------------------------------------------------------------*/
/* Purpose:  This routine performs unsharp masking of a scalar image.        */
/*---------------------------------------------------------------------------*/
void UnsharpScalar()
{
   FLOAT_TYPE **Data1;
   FLOAT_TYPE **Data2;
   FLOAT_TYPE **Data3;
   FLOAT_TYPE **Temp;
   int i, x, y;
   float Min, Max, Scale, Value;

   /* Read input image */
   Data1 = (FLOAT_TYPE **) im_alloc2D(Image1, FLOAT);
   im_read(Image1, FLOAT, (char *) &(Data1[0][0]));

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

   /* Calculate min and max */
   Min = Max = Data1[0][0];
   for (y = 0; y < Ydim; y++)
      for (x = 0; x < Xdim; x++)
      {
	 if (Min > Data1[y][x])
	    Min = Data1[y][x];
	 if (Max < Data1[y][x])
	    Max = Data1[y][x];
      }

   /* Scale the image to 0..255 range */
   Scale = 255.0 / (Max - Min + 0.01);
   for (y = 0; y < Ydim; y++)
      for (x = 0; x < Xdim; x++)
	 Data1[y][x] = (Data1[y][x] - Min) * Scale;

   /* Perform binomial smoothing */
   for (i = 0; i < Smooth; i++)
   {
      /* Save sharp image */
      if (i == Sharp)
	 for (y = 0; y < Ydim; y++)
	    for (x = 0; x < Xdim; x++)
	       Data3[y][x] = Data1[y][x];

      /* Average pixels */
      for (y = 1; y < (Ydim - 1); y++)
	 for (x = 1; x < (Xdim - 1); x++)
	    Data2[y][x] = (Data1[y - 1][x - 1] + Data1[y - 1][x + 1]
			   + Data1[y + 1][x - 1] + Data1[y + 1][x + 1]
			   + (Data1[y - 1][x] + Data1[y + 1][x]
			      + Data1[y][x - 1] + Data1[y][x + 1]) * 2
			   + Data1[y][x] * 4) / 16;

      /* Fix boundary */
      for (y = 0; y < Ydim; y++)
      {
	 Data2[y][0] = Data1[y][1];
	 Data2[y][Xdim - 1] = Data1[y][Xdim - 2];
      }
      for (x = 0; x < Xdim; x++)
      {
	 Data2[0][x] = Data1[1][x];
	 Data2[Ydim - 1][x] = Data1[Ydim - 2][x];
      }

      /* Swap image buffers */
      Temp = Data1;
      Data1 = Data2;
      Data2 = Temp;
   }

   /* Save sharp-smooth difference image */
   for (y = 0; y < Ydim; y++)
      for (x = 0; x < Xdim; x++)
      {
	 Value = (Data3[y][x] - Data1[y][x] * Fraction) / (1 - Fraction);
	 if (Value < 0)
	    Value = 0;
	 if (Value > 255)
	    Value = 255;
	 Data3[y][x] = Value;
      }

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

/*---------------------------------------------------------------------------*/
/* Purpose:  This routine performs unsharp masking of a vector image.        */
/*---------------------------------------------------------------------------*/
void UnsharpVector()
{
   FLOAT_TYPE ***Data1;
   FLOAT_TYPE ***Data2;
   FLOAT_TYPE ***Data3;
   FLOAT_TYPE ***Temp;
   int i, x, y, z;
   float Min, Max, Scale, Value;

   /* Read input image */
   Data1 = (FLOAT_TYPE ***) im_alloc3D(Image1, FLOAT);
   im_read(Image1, FLOAT, (char *) &(Data1[0][0][0]));

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

   /* Calculate min and max */
   Min = Max = Data1[0][0][0];
   for (z = 0; z < Zdim; z++)
      for (y = 0; y < Ydim; y++)
	 for (x = 0; x < Xdim; x++)
	 {
	    if (Min > Data1[z][y][x])
	       Min = Data1[z][y][x];
	    if (Max < Data1[z][y][x])
	       Max = Data1[z][y][x];
	 }

   /* Scale the image to 0..255 range */
   Scale = 255.0 / (Max - Min + 0.01);
   for (z = 0; z < Zdim; z++)
      for (y = 0; y < Ydim; y++)
	 for (x = 0; x < Xdim; x++)
	    Data1[z][y][x] = (Data1[z][y][x] - Min) * Scale;

   /* Perform binomial smoothing */
   for (i = 0; i < Smooth; i++)
   {
      /* Save sharp image */
      if (i == Sharp)
	 for (z = 0; z < Zdim; z++)
	    for (y = 0; y < Ydim; y++)
	       for (x = 0; x < Xdim; x++)
		  Data3[z][y][x] = Data1[z][y][x];

      /* Average pixels */
      for (z = 0; z < Zdim; z++)
	 for (y = 1; y < (Ydim - 1); y++)
	    for (x = 1; x < (Xdim - 1); x++)
	       Data2[z][y][x] = (Data1[z][y - 1][x - 1] + Data1[z][y - 1][x + 1]
			   + Data1[z][y + 1][x - 1] + Data1[z][y + 1][x + 1]
				 + (Data1[z][y - 1][x] + Data1[z][y + 1][x]
			      + Data1[z][y][x - 1] + Data1[z][y][x + 1]) * 2
				 + Data1[z][y][x] * 4) / 16;

      /* Fix boundary */
      for (z = 0; z < Zdim; z++)
	 for (y = 0; y < Ydim; y++)
	 {
	    Data2[z][y][0] = Data1[z][y][1];
	    Data2[z][y][Xdim - 1] = Data1[z][y][Xdim - 2];
	 }
      for (z = 0; z < Zdim; z++)
	 for (x = 0; x < Xdim; x++)
	 {
	    Data2[z][0][x] = Data1[z][1][x];
	    Data2[z][Ydim - 1][x] = Data1[z][Ydim - 2][x];
	 }

      /* Swap image buffers */
      Temp = Data1;
      Data1 = Data2;
      Data2 = Temp;
   }

   /* Save sharp-smooth difference image */
   for (z = 0; z < Zdim; z++)
      for (y = 0; y < Ydim; y++)
	 for (x = 0; x < Xdim; x++)
	 {
	    Value = (Data3[z][y][x] - Data1[z][y][x] * Fraction)
	       / (1 - Fraction);
	    if (Value < 0)
	       Value = 0;
	    if (Value > 255)
	       Value = 255;
	    Data3[z][y][x] = Value;
	 }

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

/*---------------------------------------------------------------------------*/
/* Purpose:  This routine performs unsharp masking of a color image.         */
/*---------------------------------------------------------------------------*/
void UnsharpColor()
{
   COLOR_TYPE **Data1;
   COLOR_TYPE **Data2;
   COLOR_TYPE **Data3;
   COLOR_TYPE **Temp;
   int i, x, y, temp;

   /* Read input image */
   Data1 = (COLOR_TYPE **) im_alloc2D(Image1, COLOR);
   im_read(Image1, COLOR, (char *) &(Data1[0][0]));

   /* Create output image */
   Image2 = im_create(Name2, PixType, Xdim, Ydim, Zdim);
   Data2 = (COLOR_TYPE **) im_alloc2D(Image2, COLOR);
   Data3 = (COLOR_TYPE **) im_alloc2D(Image2, COLOR);

   /* Perform binomial smoothing */
   for (i = 0; i < Smooth; i++)
   {
      /* Save sharp image */
      if (i == Sharp)
	 for (y = 0; y < Ydim; y++)
	    for (x = 0; x < Xdim; x++)
	    {
	       Data3[y][x].r = Data1[y][x].r;
	       Data3[y][x].g = Data1[y][x].g;
	       Data3[y][x].b = Data1[y][x].b;
	    }

      /* Average pixels */
      for (y = 1; y < (Ydim - 1); y++)
	 for (x = 1; x < (Xdim - 1); x++)
	 {
	    temp = ((int) Data1[y - 1][x - 1].r + (int) Data1[y - 1][x + 1].r
		 + (int) Data1[y + 1][x - 1].r + (int) Data1[y + 1][x + 1].r
		    + ((int) Data1[y - 1][x].r + (int) Data1[y + 1][x].r
		    + (int) Data1[y][x - 1].r + (int) Data1[y][x + 1].r) * 2
		    + (int) Data1[y][x].r * 4) / 16;
	    if (temp < 0)
	       temp = 0;
	    else if (temp > 255)
	       temp = 255;
	    Data2[y][x].r = temp;

	    temp = ((int) Data1[y - 1][x - 1].g + (int) Data1[y - 1][x + 1].g
		 + (int) Data1[y + 1][x - 1].g + (int) Data1[y + 1][x + 1].g
		    + ((int) Data1[y - 1][x].g + (int) Data1[y + 1][x].g
		    + (int) Data1[y][x - 1].g + (int) Data1[y][x + 1].g) * 2
		    + (int) Data1[y][x].g * 4) / 16;
	    if (temp < 0)
	       temp = 0;
	    else if (temp > 255)
	       temp = 255;
	    Data2[y][x].g = temp;

	    temp = ((int) Data1[y - 1][x - 1].b + (int) Data1[y - 1][x + 1].b
		 + (int) Data1[y + 1][x - 1].b + (int) Data1[y + 1][x + 1].b
		    + ((int) Data1[y - 1][x].b + (int) Data1[y + 1][x].b
		    + (int) Data1[y][x - 1].b + (int) Data1[y][x + 1].b) * 2
		    + (int) Data1[y][x].b * 4) / 16;
	    if (temp < 0)
	       temp = 0;
	    else if (temp > 255)
	       temp = 255;
	    Data2[y][x].b = temp;
	 }

      /* Fix boundary */
      for (y = 0; y < Ydim; y++)
      {
	 Data2[y][0].r = Data1[y][1].r;
	 Data2[y][Xdim - 1].r = Data1[y][Xdim - 2].r;
	 Data2[y][0].g = Data1[y][1].g;
	 Data2[y][Xdim - 1].g = Data1[y][Xdim - 2].g;
	 Data2[y][0].b = Data1[y][1].b;
	 Data2[y][Xdim - 1].b = Data1[y][Xdim - 2].b;
      }

      for (x = 0; x < Xdim; x++)
      {
	 Data2[0][x].r = Data1[1][x].r;
	 Data2[Ydim - 1][x].r = Data1[Ydim - 2][x].r;
	 Data2[0][x].g = Data1[1][x].g;
	 Data2[Ydim - 1][x].g = Data1[Ydim - 2][x].g;
	 Data2[0][x].b = Data1[1][x].b;
	 Data2[Ydim - 1][x].b = Data1[Ydim - 2][x].b;
      }

      /* Swap image buffers */
      Temp = Data1;
      Data1 = Data2;
      Data2 = Temp;
   }

   /* Save sharp-smooth difference image */
   for (y = 0; y < Ydim; y++)
      for (x = 0; x < Xdim; x++)
      {
	 temp = (Data3[y][x].r - Data1[y][x].r * Fraction) / (1 - Fraction);
	 if (temp < 0)
	    temp = 0;
	 else if (temp > 255)
	    temp = 255;
	 Data3[y][x].r = temp;
	 temp = (Data3[y][x].g - Data1[y][x].g * Fraction) / (1 - Fraction);
	 if (temp < 0)
	    temp = 0;
	 else if (temp > 255)
	    temp = 255;
	 Data3[y][x].g = temp;
	 temp = (Data3[y][x].b - Data1[y][x].b * Fraction) / (1 - Fraction);
	 if (temp < 0)
	    temp = 0;
	 else if (temp > 255)
	    temp = 255;
	 Data3[y][x].b = temp;
      }

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

/*---------------------------------------------------------------------------*/
/* Purpose:  This is the main program.                                       */
/*---------------------------------------------------------------------------*/
int main(int argc, char *argv[])
{
   /* Program variables */
   int i = 0;

   /* Interpret program options */
   printf("UNSHARP Program - KUIM Version 2.1\n\n");
   while ((++i < argc) && (argv[i][0] == '-'))
      switch (argv[i][1])
      {
      case 'k':
	 if (sscanf(argv[++i], "%d", &Sharp) == 0)
	    Error("Could not get number of sharp iterations");
	 if (sscanf(argv[++i], "%d", &Smooth) == 0)
	    Error("Could not get number of smooth iterations");
	 if (Sharp >= Smooth)
	    Error("Must have sharp iterations < smooth iterations");
	 break;
      case 'f':
	 if (sscanf(argv[++i], "%f", &Fraction) == 0)
	    Error("Could not get fraction parameter");
	 if ((Fraction <= 0) || (Fraction >= 1))
	    Error("Must have fraction between 0 and 1");
      case 'd':
	 Debug = TRUE;
	 break;
      default:
	 Error("Invalid option encountered");
	 break;
      }

   /* Check smoothing parameters */
   if (Sharp >= Smooth)
      Error("Number of sharp iterations must be less than smooth iterations");

   /* Check number of file names */
   if (argc - i != 2)
   {
      fprintf(stderr, "Usage: unsharp [options] infile outfile\n");
      fprintf(stderr, "\t[-d]      Print debugging information\n");
      fprintf(stderr, "\t[-k # #]  Number of smoothing iterations (1,5)\n");
      fprintf(stderr, "\t[-f #]    Image subtraction fraction (0.8)\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");

   /* Determine PixType */
   Image1 = im_open(Name1, &PixType, &Xdim, &Ydim, &Zdim, &DimCnt);

   /* Handle color images */
   if (PixType == COLOR || PixType == PSEUDO || PixType == JPEG_RGB)
      UnsharpColor();

   /* Handle scalar images */
   else if (DimCnt == 2)
      UnsharpScalar();

   /* Handle vector images */
   else if (DimCnt == 3)
      UnsharpVector();

   /* Unsupported image type */
   else
      Error("Can not process image");
   return (0);
}