/*---------------------------------------------------------------------------*/
/* Program:  interp.c                                                        */
/*                                                                           */
/* Purpose:  This program performs scaling of KUIM format images.            */
/*           Floating point interpolation is used for accuracy and not       */
/*           speed.  Fast interpolation methods are used for scale           */
/*           factors of 0.5 and 2.0.  It also works for color images.        */
/*                                                                           */
/* Author:   John Gauch / Chunyen Liu                                        */
/*                                                                           */
/* Date:     March 1, 1996                                                   */
/*                                                                           */
/* Note:     Copyright (C) The University of Kansas, 1996                    */
/*---------------------------------------------------------------------------*/
#include <IM.h>

/* Global variables */
FLOAT_TYPE **Data1;
FLOAT_TYPE **Data2;
COLOR_TYPE **Data3;
COLOR_TYPE **Data4;
int Xdim, Ydim, Zdim;
int NewXdim = 0;
int NewYdim = 0;

/*---------------------------------------------------------------------------*/
/* Purpose:  This routine quickly handles expanding by factors of 2.         */
/*---------------------------------------------------------------------------*/
void Expand()
{
   int x, y;

   /* Linearly interpolate between pixels */
   for (y = 0; y < (Ydim - 1); y++)
      for (x = 0; x < (Xdim - 1); x++)
      {
	 Data2[y + y][x + x] = Data1[y][x];
	 Data2[y + y][x + x + 1] = (Data1[y][x] + Data1[y][x + 1]) / 2;
	 Data2[y + y + 1][x + x] = (Data1[y + 1][x] + Data1[y][x]) / 2;
	 Data2[y + y + 1][x + x + 1] = (Data1[y + 1][x] + Data1[y][x + 1]) / 2;
      }

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

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

   /* Handle last pixel */
   Data2[NewYdim - 2][NewXdim - 2] = Data1[Ydim - 1][Xdim - 1];
   Data2[NewYdim - 1][NewXdim - 2] = Data1[Ydim - 1][Xdim - 1];
   Data2[NewYdim - 2][NewXdim - 1] = Data1[Ydim - 1][Xdim - 1];
   Data2[NewYdim - 1][NewXdim - 1] = Data1[Ydim - 1][Xdim - 1];
}

/*---------------------------------------------------------------------------*/
/* Purpose:  This routine handles expanding by factors of 2 in color.        */
/*---------------------------------------------------------------------------*/
void ExpandColor()
{
   int x, y;

   /* Linearly interpolate between color pixels */
   for (y = 0; y < (Ydim - 1); y++)
      for (x = 0; x < (Xdim - 1); x++)
      {
	 Data4[y + y][x + x].r = Data3[y][x].r;
	 Data4[y + y][x + x].g = Data3[y][x].g;
	 Data4[y + y][x + x].b = Data3[y][x].b;
	 Data4[y + y][x + x + 1].r = (Data3[y][x].r + Data3[y][x + 1].r) / 2;
	 Data4[y + y][x + x + 1].g = (Data3[y][x].g + Data3[y][x + 1].g) / 2;
	 Data4[y + y][x + x + 1].b = (Data3[y][x].b + Data3[y][x + 1].b) / 2;
	 Data4[y + y + 1][x + x].r = (Data3[y + 1][x].r + Data3[y][x].r) / 2;
	 Data4[y + y + 1][x + x].g = (Data3[y + 1][x].g + Data3[y][x].g) / 2;
	 Data4[y + y + 1][x + x].b = (Data3[y + 1][x].b + Data3[y][x].b) / 2;
	 Data4[y + y + 1][x + x + 1].r = (Data3[y + 1][x].r + Data3[y][x + 1].r) / 2;
	 Data4[y + y + 1][x + x + 1].g = (Data3[y + 1][x].g + Data3[y][x + 1].g) / 2;
	 Data4[y + y + 1][x + x + 1].b = (Data3[y + 1][x].b + Data3[y][x + 1].b) / 2;
      }

   /* Handle last column */
   x = Xdim - 1;
   for (y = 0; y < (Ydim - 1); y++)
   {
      Data4[y + y][x + x].r = Data4[y + y][x + x + 1].r = Data3[y][x].r;
      Data4[y + y][x + x].g = Data4[y + y][x + x + 1].g = Data3[y][x].g;
      Data4[y + y][x + x].b = Data4[y + y][x + x + 1].b = Data3[y][x].b;
      Data4[y + y + 1][x + x].r = Data4[y + y + 1][x + x + 1].r
	 = (Data3[y + 1][x].r + Data3[y][x].r) / 2;
      Data4[y + y + 1][x + x].g = Data4[y + y + 1][x + x + 1].g
	 = (Data3[y + 1][x].g + Data3[y][x].g) / 2;
      Data4[y + y + 1][x + x].b = Data4[y + y + 1][x + x + 1].b
	 = (Data3[y + 1][x].b + Data3[y][x].b) / 2;
   }

   /* Handle last row */
   y = Ydim - 1;
   for (x = 0; x < (Xdim - 1); x++)
   {
      Data4[y + y][x + x].r = Data4[y + y + 1][x + x].r = Data3[y][x].r;
      Data4[y + y][x + x].g = Data4[y + y + 1][x + x].g = Data3[y][x].g;
      Data4[y + y][x + x].b = Data4[y + y + 1][x + x].b = Data3[y][x].b;
      Data4[y + y][x + x + 1].r = Data4[y + y + 1][x + x + 1].r
	 = (Data3[y][x + 1].r + Data3[y][x].r) / 2;
      Data4[y + y][x + x + 1].g = Data4[y + y + 1][x + x + 1].g
	 = (Data3[y][x + 1].g + Data3[y][x].g) / 2;
      Data4[y + y][x + x + 1].b = Data4[y + y + 1][x + x + 1].b
	 = (Data3[y][x + 1].b + Data3[y][x].b) / 2;
   }

   /* Handle last pixel */
   Data4[NewYdim - 2][NewXdim - 2].r = Data3[Ydim - 1][Xdim - 1].r;
   Data4[NewYdim - 2][NewXdim - 2].g = Data3[Ydim - 1][Xdim - 1].g;
   Data4[NewYdim - 2][NewXdim - 2].b = Data3[Ydim - 1][Xdim - 1].b;
   Data4[NewYdim - 1][NewXdim - 2].r = Data3[Ydim - 1][Xdim - 1].r;
   Data4[NewYdim - 1][NewXdim - 2].g = Data3[Ydim - 1][Xdim - 1].g;
   Data4[NewYdim - 1][NewXdim - 2].b = Data3[Ydim - 1][Xdim - 1].b;
   Data4[NewYdim - 2][NewXdim - 1].r = Data3[Ydim - 1][Xdim - 1].r;
   Data4[NewYdim - 2][NewXdim - 1].g = Data3[Ydim - 1][Xdim - 1].g;
   Data4[NewYdim - 2][NewXdim - 1].b = Data3[Ydim - 1][Xdim - 1].b;
   Data4[NewYdim - 1][NewXdim - 1].r = Data3[Ydim - 1][Xdim - 1].r;
   Data4[NewYdim - 1][NewXdim - 1].g = Data3[Ydim - 1][Xdim - 1].g;
   Data4[NewYdim - 1][NewXdim - 1].b = Data3[Ydim - 1][Xdim - 1].b;
}

/*---------------------------------------------------------------------------*/
/* Purpose:  This routine quickly handles contracting by factors of 2.       */
/*---------------------------------------------------------------------------*/
void Contract()
{
   int x, y;

   /* Perform 2x2 averaging */
   for (y = 0; y < NewYdim; y++)
      for (x = 0; x < NewXdim; x++)
	 Data2[y][x] = (Data1[y + y][x + x] + Data1[y + y][x + x + 1]
	       + Data1[y + y + 1][x + x] + Data1[y + y + 1][x + x + 1]) / 4;

   /* Perform subsampling */
   for (y = 0; y < NewYdim; y++)
      for (x = 0; x < NewXdim; x++)
	 Data2[y][x] = Data1[y + y][x + x];
}

/*---------------------------------------------------------------------------*/
/* Purpose:  This routine handles contracting by factors of 2 in color.      */
/*---------------------------------------------------------------------------*/
void ContractColor()
{
   int x, y;

   /* Perform 2x2 averaging */
   for (y = 0; y < NewYdim; y++)
      for (x = 0; x < NewXdim; x++)
      {
	 Data4[y][x].r = (Data3[y + y][x + x].r + Data3[y + y][x + x + 1].r
	   + Data3[y + y + 1][x + x].r + Data3[y + y + 1][x + x + 1].r) / 4;
	 Data4[y][x].g = (Data3[y + y][x + x].g + Data3[y + y][x + x + 1].g
	   + Data3[y + y + 1][x + x].g + Data3[y + y + 1][x + x + 1].g) / 4;
	 Data4[y][x].b = (Data3[y + y][x + x].b + Data3[y + y][x + x + 1].b
	   + Data3[y + y + 1][x + x].b + Data3[y + y + 1][x + x + 1].b) / 4;
      }

   /* Perform subsampling */
   for (y = 0; y < NewYdim; y++)
      for (x = 0; x < NewXdim; x++)
      {
	 Data4[y][x].r = Data3[y + y][x + x].r;
	 Data4[y][x].g = Data3[y + y][x + x].g;
	 Data4[y][x].b = Data3[y + y][x + x].b;
      }
}

/*---------------------------------------------------------------------------*/
/* Purpose:  This routine handles general interpolation.                     */
/*---------------------------------------------------------------------------*/
void Interp()
{
   int x, y;
   float RealX, RealY, StepX, StepY;
   int IntX, IntY;
   float Low, High;
   IM_TYPE *Image6;
   FLOAT_TYPE **Data6;

   /* Allocate temporary image */
   Image6 = im_create("/dev/null", FLOAT, NewXdim, Ydim, Zdim);
   Data6 = (FLOAT_TYPE **) im_alloc2D(Image6, FLOAT);

   /* First pass - expand in X direction */
   StepX = (float) Xdim / (float) NewXdim;
   for (y = 0; y < Ydim; y++)
      for (x = 0; x < NewXdim; x++)
      {
	 RealX = x * StepX;
	 IntX = (int) RealX;
	 if (IntX + 1 < Xdim)
	 {
	    Low = Data1[y][IntX];
	    High = Data1[y][IntX + 1];
	    Data6[y][x] = Low + (High - Low) * (RealX - IntX);
	 }
	 else
	    Data6[y][x] = Data1[y][IntX];
      }

   /* Second pass - expand in Y direction */
   StepY = (float) Ydim / (float) NewYdim;
   for (y = 0; y < NewYdim; y++)
      for (x = 0; x < NewXdim; x++)
      {
	 RealY = y * StepY;
	 IntY = (int) RealY;
	 if (IntY + 1 < Ydim)
	 {
	    Low = Data6[IntY][x];
	    High = Data6[IntY + 1][x];
	    Data2[y][x] = Low + (High - Low) * (RealY - IntY);
	 }
	 else
	    Data2[y][x] = Data6[IntY][x];
      }
   im_free2D((char **) Data6);
}

/*---------------------------------------------------------------------------*/
/* Purpose:  This routine handles general interpolation in color.            */
/*---------------------------------------------------------------------------*/
void InterpColor()
{
   int x, y;
   float RealX, RealY, StepX, StepY;
   int IntX, IntY;
   float LowR, HighR, LowG, HighG, LowB, HighB;
   IM_TYPE *Image5;
   COLOR_TYPE **Data5;

   /* Allocate temporary image */
   Image5 = im_create("/dev/null", COLOR, NewXdim, Ydim, Zdim);
   Data5 = (COLOR_TYPE **) im_alloc2D(Image5, COLOR);

   /* First pass - expand in X direction */
   StepX = (float) Xdim / (float) NewXdim;
   for (y = 0; y < Ydim; y++)
      for (x = 0; x < NewXdim; x++)
      {
	 RealX = x * StepX;
	 IntX = (int) RealX;
	 if (IntX + 1 < Xdim)
	 {
	    LowR = Data3[y][IntX].r;
	    LowG = Data3[y][IntX].g;
	    LowB = Data3[y][IntX].b;
	    HighR = Data3[y][IntX + 1].r;
	    HighG = Data3[y][IntX + 1].g;
	    HighB = Data3[y][IntX + 1].b;
	    Data5[y][x].r = (unsigned char) (LowR + (HighR - LowR) * (RealX - IntX));
	    Data5[y][x].g = (unsigned char) (LowG + (HighG - LowG) * (RealX - IntX));
	    Data5[y][x].b = (unsigned char) (LowB + (HighB - LowB) * (RealX - IntX));
	 }
	 else
	 {
	    Data5[y][x].r = Data3[y][IntX].r;
	    Data5[y][x].g = Data3[y][IntX].g;
	    Data5[y][x].b = Data3[y][IntX].b;
	 }
      }

   /* Second pass - expand in Y direction */
   StepY = (float) Ydim / (float) NewYdim;
   for (y = 0; y < NewYdim; y++)
      for (x = 0; x < NewXdim; x++)
      {
	 RealY = y * StepY;
	 IntY = (int) RealY;
	 if (IntY + 1 < Ydim)
	 {
	    LowR = Data5[IntY][x].r;
	    LowG = Data5[IntY][x].g;
	    LowB = Data5[IntY][x].b;
	    HighR = Data5[IntY + 1][x].r;
	    HighG = Data5[IntY + 1][x].g;
	    HighB = Data5[IntY + 1][x].b;
	    Data4[y][x].r = (unsigned char) (LowR + (HighR - LowR) * (RealY - IntY));
	    Data4[y][x].g = (unsigned char) (LowG + (HighG - LowG) * (RealY - IntY));
	    Data4[y][x].b = (unsigned char) (LowB + (HighB - LowB) * (RealY - IntY));
	 }
	 else
	 {
	    Data4[y][x].r = Data5[IntY][x].r;
	    Data4[y][x].g = Data5[IntY][x].g;
	    Data4[y][x].b = Data5[IntY][x].b;
	 }
      }
   im_free2D((char **) Data5);
}

/*---------------------------------------------------------------------------*/
/* 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 Debug = FALSE;
   float Scale = 1.0;
   int i = 0;

   /* Interpret program options */
   printf("INTERP Program - KUIM Version 2.0\n\n");
   while ((++i < argc) && (argv[i][0] == '-'))
      switch (argv[i][1])
      {
      case 's':
	 if (sscanf(argv[++i], "%f", &Scale) == 0)
	    Error("Could not get Scale argument");
	 break;
      case 'x':
	 if (sscanf(argv[++i], "%d", &NewXdim) == 0)
	    Error("Could not get new X dimension");
	 break;
      case 'y':
	 if (sscanf(argv[++i], "%d", &NewYdim) == 0)
	    Error("Could not get new Y dimension");
	 break;
      case 'd':
	 Debug = TRUE;
	 break;
      default:
	 Error("Invalid option encountered");
	 break;
      }

   /* Check number of file names */
   if (argc - i != 2)
   {
      fprintf(stderr, "Usage: interp [options] infile outfile\n");
      fprintf(stderr, "              [-d]      Print debugging information\n");
      fprintf(stderr, "              [-s #]    Image scaling factor\n");
      fprintf(stderr, "              [-x #]    New size in X dimension\n");
      fprintf(stderr, "              [-y #]    New size in Y dimension\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 == 1)
      Error("Can not process 1D images");
   if ((PixType == PSEUDO) || (PixType == COLOR) || (PixType == JPEG_RGB))
   {
      Data3 = (COLOR_TYPE **) im_alloc2D(Image1, COLOR);
      im_read(Image1, COLOR, (char *) &(Data3[0][0]));
   }
   else
   {
      Data1 = (FLOAT_TYPE **) im_alloc2D(Image1, FLOAT);
      im_read(Image1, FLOAT, (char *) &(Data1[0][0]));
   }

   /* Check new image dimensions */
   if (NewXdim == 0 && NewYdim == 0 && Scale == 1.0)
      Error("Must specify new X and Y dimensions");
   if (NewXdim == 0)
      NewXdim = (int) (Xdim * Scale);
   if (NewYdim == 0)
      NewYdim = (int) (Ydim * Scale);
   if (Debug == TRUE)
      printf("Old image dimensions %d %d\n", Xdim, Ydim);
   if (Debug == TRUE)
      printf("New image dimensions %d %d\n", NewXdim, NewYdim);

   /* Create output image */
   Image2 = im_create(Name2, PixType, NewXdim, NewYdim, Zdim);
   if ((PixType == PSEUDO) || (PixType == COLOR) || (PixType == JPEG_RGB))
      Data4 = (COLOR_TYPE **) im_alloc2D(Image2, COLOR);
   else
      Data2 = (FLOAT_TYPE **) im_alloc2D(Image2, FLOAT);

   /* Handle 3D images by unfolding in Y dimension */
   Ydim *= Zdim;
   NewYdim *= Zdim;

   /* Perform interpolation */
   if ((NewXdim == 2 * Xdim) && (NewYdim == 2 * Ydim))
   {
      if ((PixType == PSEUDO) || (PixType == COLOR) || (PixType == JPEG_RGB))
	 ExpandColor();
      else
	 Expand();
   }
   else if ((Xdim == 2 * NewXdim) && (Ydim == 2 * NewYdim))
   {
      if ((PixType == PSEUDO) || (PixType == COLOR) || (PixType == JPEG_RGB))
	 ContractColor();
      else
	 Contract();
   }
   else
   {
      if ((PixType == PSEUDO) || (PixType == COLOR) || (PixType == JPEG_RGB))
	 InterpColor();
      else
	 Interp();
   }

   /* Write information to output image */
   if ((PixType == PSEUDO) || (PixType == COLOR) || (PixType == JPEG_RGB))
   {
      im_write(Image2, COLOR, (char *) &(Data4[0][0]));
      im_free2D((char **) Data3);
      im_free2D((char **) Data4);
   }
   else
   {
      im_write(Image2, FLOAT, (char *) &(Data2[0][0]));
      im_free2D((char **) Data1);
      im_free2D((char **) Data2);
   }
   return (0);
}