/*---------------------------------------------------------------------------*/
/* Program:  butter.c                                                        */
/*                                                                           */
/* Purpose:  This program performs Butterworth filtering in the frequency    */
/*           domain.  Images are automatically converted to COMPLEX type     */
/*           filtered, and converted back to FLOAT type.                     */
/*                                                                           */
/* Author:   John Gauch                                                      */
/*                                                                           */
/* Date:     October 7, 1994                                                 */
/*                                                                           */
/* Note:     Copyright (C) The University of Kansas, 1994                    */
/*---------------------------------------------------------------------------*/
#include <IM.h>

/* Global variables */
int LowPass = TRUE;
int Homomorph = FALSE;
float Power = 2.0;
float CutOff = 20.0;
float Boost = 0.4;

/*---------------------------------------------------------------------------*/
/* Purpose:  This routine performs Butterworth filtering.                    */
/*---------------------------------------------------------------------------*/
void Filter(COMPLEX_TYPE ** Data, int Xdim, int Ydim)
{
   int x, y, x1, y1, halfx, halfy;
   float MinVal, Filter;

   /* Take log of input image */
   if (Homomorph == TRUE)
   {
      MinVal = Data[0][0].re;
      for (y = 0; y < Ydim; y++)
	 for (x = 0; x < Xdim; x++)
	    if (MinVal > Data[y][x].re)
	       MinVal = Data[y][x].re;
      for (y = 0; y < Ydim; y++)
	 for (x = 0; x < Xdim; x++)
	    Data[y][x].re = (float) log((double) (1.0 + Data[y][x].re - MinVal));
   }

   /* Convert to frequency domain */
   if (fft_2d(&(Data[0][0]), Xdim, Ydim, -1) == INVALID)
      Error("Could not perform FFT");

   /* Prepare to filter */
   halfx = Xdim / 2;
   halfy = Ydim / 2;

   /* Loop over Y axis */
   for (y = 0; y < Ydim; y++)
   {
      if (y < halfy)
	 y1 = y;
      else
	 y1 = y - Ydim;

      /* Loop over X axis */
      for (x = 0; x < Xdim; x++)
      {
	 if (x < halfx)
	    x1 = x;
	 else
	    x1 = x - Xdim;

	 /* Calculate value of Butterworth filter */
	 if (LowPass == TRUE)
	    Filter = (float) (1 / (1 + pow((x1 * x1 + y1 * y1) /
					   (CutOff * CutOff), Power)));
	 else if ((x1 != 0) || (y1 != 0))
	    Filter = (float) (1 / (1 + pow((CutOff * CutOff) /
					   (x1 * x1 + y1 * y1), Power)));
	 else
	    Filter = (float) 0.0;
	 if (Homomorph == TRUE)
	    Filter = Boost + (1 - Boost) * Filter;

	 /* Do pointwise multiplication */
	 Data[y][x].re *= Filter;
	 Data[y][x].im *= Filter;
      }
   }

   /* Convert to spatial domain */
   if (fft_2d(&(Data[0][0]), Xdim, Ydim, 1) == INVALID)
      Error("Could not perform inverse FFT");

   /* Take exp of input image */
   if (Homomorph == TRUE)
      for (y = 0; y < Ydim; y++)
	 for (x = 0; x < Xdim; x++)
	    Data[y][x].re = (float) (exp((double) Data[y][x].re) - 1.0 + MinVal);
}

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

   /* Program variables */
   int Debug = FALSE;
   int i = 0;

   /* Interpret program options */
   printf("BUTTER Program - KUIM Version 1.0\n\n");
   while ((++i < argc) && (argv[i][0] == '-'))
      switch (argv[i][1])
      {
      case 'l':
	 LowPass = TRUE;
	 Homomorph = FALSE;
	 break;
      case 'h':
	 LowPass = FALSE;
	 Homomorph = FALSE;
	 break;
      case 'L':
	 LowPass = TRUE;
	 Homomorph = TRUE;
	 break;
      case 'H':
	 LowPass = FALSE;
	 Homomorph = TRUE;
	 break;
      case 'f':
	 if (sscanf(argv[++i], "%f", &CutOff) == 0)
	    Error("Could not get CutOff argument");
	 if (CutOff < 1)
	    Error("Cut off frequency must be >= 1");
	 break;
      case 'n':
	 if (sscanf(argv[++i], "%f", &Power) == 0)
	    Error("Could not get Power argument");
	 if (Power < 0)
	    Error("Power must be >= 0");
	 break;
      case 'b':
	 if (sscanf(argv[++i], "%f", &Boost) == 0)
	    Error("Could not get Boost argument");
	 if ((Boost < 0) || (Boost > 1))
	    Error("Boost must be within [0..1]");
	 break;
      case 'd':
	 Debug = TRUE;
	 break;
      default:
	 Error("Invalid option encountered");
	 break;
      }

   /* Check number of file names */
   if (argc - i != 2)
   {
      fprintf(stderr, "Usage: butter [options] infile outfile\n");
      fprintf(stderr, "       [-l]    Lowpass filtering\n");
      fprintf(stderr, "       [-h]    Highpass filtering\n");
      fprintf(stderr, "       [-L]    Homomorphic lowpass filtering\n");
      fprintf(stderr, "       [-H]    Homomorphic highpass filtering\n");
      fprintf(stderr, "       [-f #]  Frequency with filter = 1/2 (20.0)\n");
      fprintf(stderr, "       [-n #]  Denominator exponent (2.0)\n");
      fprintf(stderr, "       [-b #]  Filter boost (Homomorphic only) (0.4)\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 and create output image */
   Image1 = im_open(Name1, &PixType, &Xdim, &Ydim, &Zdim, &DimCnt);
   if (DimCnt != 2)
      Error("Can not process 1D or 3D images");
   Image2 = im_create(Name2, FLOAT, Xdim, Ydim, Zdim);
   Data = (COMPLEX_TYPE **) im_alloc2D(Image1, COMPLEX);
   im_read(Image1, COMPLEX, (char *) &(Data[0][0]));

   /* Perform filtering */
   Filter(Data, Xdim, Ydim);

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