/*---------------------------------------------------------------------------*/
/* Program:  mult.c                                                          */
/*                                                                           */
/* Purpose:  This program multiplies (or divides) two KUIM images.           */
/*                                                                           */
/* Author:   John Gauch                                                      */
/*                                                                           */
/* Date:     April 27, 1994   - Original program.                            */
/*           October 13, 1994 - Added supprort for complex.                  */
/*                                                                           */
/* Note:     Copyright (C) The University of Kansas, 1994                    */
/*---------------------------------------------------------------------------*/
#include <IM.h>
#define THRESH 0.005

/*---------------------------------------------------------------------------*/
/* Purpose:  This is the main program.                                       */
/*---------------------------------------------------------------------------*/
int main(int argc, char *argv[])
{
   /* Image variables */
   char Name1[50];
   char Name2[50];
   char Name3[50];
   IM_TYPE *Image1;
   IM_TYPE *Image2;
   IM_TYPE *Image3;
   COMPLEX_TYPE *Data1;
   COMPLEX_TYPE *Data2;
   COMPLEX_TYPE *Data3;
   int PixType, Xdim, Ydim, Zdim, DimCnt;
   int Xdim2, Ydim2, Zdim2;

   /* Program variables */
   int Quotient = FALSE;
   int Debug = FALSE;
   int i = 0;
   int Count;
   float temp;

   /* Interpret program options */
   printf("MULT Program - KUIM Version 2.0\n\n");
   while ((++i < argc) && (argv[i][0] == '-'))
      switch (argv[i][1])
      {
      case 'q':
	 Quotient = TRUE;
	 break;
      case 'p':
	 Quotient = FALSE;
	 break;
      case 'd':
	 Debug = TRUE;
	 break;
      default:
	 Error("Invalid option encountered");
	 break;
      }

   /* Check number of file names */
   if (argc - i != 3)
   {
      fprintf(stderr, "Usage: mult [options] infile1 infile2 outfile\n");
      fprintf(stderr, "       [-d]  Print debugging information\n");
      fprintf(stderr, "       [-q]  Calculate infile1 / infile2\n");
      fprintf(stderr, "       [-p]  Calculate infile1 * infile2 (default)\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 input file name");
   if (sscanf(argv[i++], "%s", Name3) == 0)
      Error("Could not get output file name");

   /* Read input image1 */
   Image1 = im_open(Name1, &PixType, &Xdim, &Ydim, &Zdim, &DimCnt);
   if (PixType == COLOR || PixType == PSEUDO)
      Error("Can not handle COLOR images");
   Data1 = (COMPLEX_TYPE *) im_alloc1D(Image1, COMPLEX);
   im_read(Image1, COMPLEX, (char *) Data1);

   /* Read input image2 */
   Image2 = im_open(Name2, &PixType, &Xdim2, &Ydim2, &Zdim2, &DimCnt);
   if ((Xdim != Xdim2) || (Ydim != Ydim2) || (Zdim != Zdim2))
      Error("Image dimensions do not match");
   Image2 = im_open(Name2, &PixType, &Xdim, &Ydim, &Zdim, &DimCnt);
   if (PixType == COLOR || PixType == PSEUDO || PixType == JPEG_RGB)
      Error("Can not handle COLOR images");
   Data2 = (COMPLEX_TYPE *) im_alloc1D(Image2, COMPLEX);
   im_read(Image2, COMPLEX, (char *) Data2);

   /* Create output image */
   Image3 = im_create(Name3, PixType, Xdim, Ydim, Zdim);
   Data3 = (COMPLEX_TYPE *) im_alloc1D(Image3, COMPLEX);
   Count = Xdim * Ydim * Zdim;

   /* Multiply two images */
   if (Quotient == FALSE)
      for (i = 0; i < Count; i++)
      {
	 Data3[i].re = Data1[i].re * Data2[i].re - Data1[i].im * Data2[i].im;
	 Data3[i].im = Data1[i].re * Data2[i].im + Data1[i].im * Data2[i].re;
      }

   /* Divide two images */
   if (Quotient == TRUE)
      for (i = 0; i < Count; i++)
      {
	 temp = Data2[i].re * Data2[i].re + Data2[i].im * Data2[i].im;
	 if (temp < THRESH)
	    temp = THRESH;
	 Data3[i].re = (Data1[i].re * Data2[i].re + Data1[i].im * Data2[i].im) / temp;
	 Data3[i].im = (Data1[i].im * Data2[i].re - Data1[i].re * Data2[i].im) / temp;
      }

   /* Write information to output image */
   im_write(Image3, COMPLEX, (char *) Data3);
   im_free1D((char *) Data1);
   im_free1D((char *) Data2);
   im_free1D((char *) Data3);
   return (0);
}