//**********************************************************************************
//
//...Modulation Transfer Function (edge trace)
//
//...C.T. Yeung
//
//...4-29-94
//
//...This program is created to determine the resolution of our current Kodak 4000.
//   However, this program will be used to determine the resolution for all of our
//   imaging equipment.
//
//...Input requirement for this program is a 1 or 2-D image of an edge.
//   Also known as step function or half of a rect
//
//...Output of this program is an array to be opened in xcel for plotting.
//                                                                                   
//...version 1
//
//***********************************************************************************
                
#include <fstream.h>
#include <math.h>    
#include <string.h>
  
const int MAXPIX=128;
const double pi=3.1416;
      
   //variables in readimage()
   char imagename[20];
   int imagewidth=0;   
   int imageheight=0;

   //variables in kernel()
   int kernelwidth=0;  
   int kernel[9];
   double coef=0.0;

   //variables in derivative()
   unsigned char image[MAXPIX];          		        //input image
   double der[MAXPIX];     	                            //output image    
   char dername[20]= "derivative";

   //variables in fourier()
   double real[MAXPIX];
   double imag[MAXPIX];                                
   double mtf[MAXPIX];        
   char mtfname[20]= "mtf"; 
   char realname[20]= "real";                        


readimage()
{
   cout << "**********************INPUT IMAGE**********************\n\n";
   cout << "\nenter the name of the file you want to mess with\n\n";
   cin >> imagename;
   
   cout << "\n\n What is the width of " << imagename << "\n\n";
   cin >> imagewidth;
         
   cout << "hello!!  I am about to read in the image...\n\n";	// test
                                        
   ifstream in (imagename);
   in.read((unsigned char*) &image, imagewidth);
   in.close();           
             
   for (int a=0; a<20; a++)                    	// test parameter
   {
      cout << "\n (in readimage) image " << a << " is " << (int)image[a];
   }             
   return 0;
}

kern()
{
   cout << "\n\n***************************KERNEL****************************\n\n";
   cout << "\n\n What is the width of the kernel? (1 to 9) \n\n";
   cin >> kernelwidth; 
   cout << "kernelwidth is "<< kernelwidth << "\n\n";           
   
   for (int y=0; y<kernelwidth; y++)
   {
      cout << "\n\n what is the coeficient for bin #   " << y+1 << "\n\n";
      cin >> kernel[y];
      
      coef += kernel[y];
   }
   cout << "\n the coefficient is \t  " << coef << "\n\n";
   cout << "\n\n the kernel you have input is as follows: \n\n";
   
   for (int x=0; x<kernelwidth; x++)
   {
      cout << "\t";
      cout << kernel [x] << "\t";
   }
   cout << "\n\n";
   
   return 0;
}  
                               
                           
derivative()
{ 
   double temp;
                                                        
   for (int a=0; a<(imagewidth-(kernelwidth-1));a++)
   {
      for (int b=0; b<kernelwidth; b++)
      {
         temp += image[a+b] * kernel[b];
      } 
      der[a] = temp;
      temp = 0;           
   }                     
   
   for (int c=0; c<20; c++)                     // test parameter
   {
      cout << "\n (in derivative) der " << c << " is " << (int)der[c];  
   }
   
   return 0;
}            
                           
                           
fourier()
{
   double gral = 0.0;

//*******************************REAL*******************************

   for (int k=-imagewidth/2; k<imagewidth/2; k++)
   {
      for (int x=-imagewidth/2; x<imagewidth/2; x++)
      {
         gral += der[x+imagewidth/2] * cos((2.0 * pi * k * x)/imagewidth);
      }
      real[k+imagewidth/2] = gral;
      gral = 0.0;
   }
   
//*******************************IMAGINARY**************************

   for (int j=-imagewidth/2; j<imagewidth/2; j++)
   {
      gral = 0.0;
      for (int y=-imagewidth/2; y<imagewidth/2; y++)
      {
         gral += der[y+imagewidth/2] * -sin((2.0 * pi * j * y)/imagewidth);   
      }
      imag [j+imagewidth/2] = gral;
   }        
   return 0;
}

magnitude()
{
   for (int x=0; x<imagewidth; x++)
   {
      mtf[x] = sqrt (real[x] * real[x] + imag[x] * imag[x]);
   }
   return 0;
}                             
                           
histout(double array[], char name[])
{
   char outfile[20];
   cout << "\n\n enter the histogram file name for the ";
   
   for (int x=0; x<20; x++)
       cout << name[x];
       
   cout << "\n\n";
   cin >> outfile;             

   ofstream hist_out(outfile);
   for (int a=0; a<imagewidth; a++)
   {
      hist_out << (double)array[a] << "\n";
   }
   return 0;                             
}   

readout(double array[], char name[])
{
   char outfile[20];
   cout << "\n\n enter the image file name for ";
   for (int i=0; i<20; i++)
      cout << name[i];
      
   cout << "\n\n";  
      
   cin >> outfile;
   ofstream out(outfile);
   out.write((unsigned char*)&array, imagewidth);
   out.close();

   return 0;
}           


main()
{            
   readimage();
   kern();
   derivative();
   fourier();
   magnitude();
   histout(mtf, mtfname);      
   /*readout(der, dername); 
   readout(mtf, mtfname);   */
   
   return 0;
}