#ifdef VMS
#include <file.h>
#else
#include <fcntl.h>
#endif
#include <math.h>
/*
********************************
**     Program: gnomconv      **
********************************

Purpose: Convert a film or image-plate image to its gnomonic projection
         using the maximum pixel method.

Author:  John Campbell, February 1992

Use: 

         gnomconv [-film] [-ip] [-size xrasts yrasts] 
                  [-cent xcen ycen] [-dist ctof] [-rmin rmin]
                  [-echo] [file_in] [file_out] 

              The user will be prompted for parameters not supplied
*/
main (argc, argv)
int argc;
char ** argv;
{
   char str[80],str2[80];
   char filin[80],filout[80];
   char *p;
   int i,j,k;
   int irast;
   int type;
   int ifil;
   int echo;
   int cent;
   int dist;
   int irmin;
   int err;
   int fin, fout;
   int iarg;
   int xrasts;
   int yrasts;
   float x;
   float y;
   float xcen;
   float ycen;
   float ctof;
   float rmin;
   unsigned char *byte_img;
   unsigned short *i2_img;
   int nbytes;

   /* Initialisations */

   filin[0] = 0;
   filout[0] = 0;
   irast = 0;
   type = 0;
   ifil = 0;
   echo = 0;
   cent = 0;
   dist = 0;
   irmin = 0;

   /* Examine parameters passed to the program */

   iarg = 1;
   while (iarg<argc)
   {
       p = argv[iarg];
       if (p[0]=='-')
       {
          switch (p[1])
          {
             case 'f':  /*film*/
             case 'F':
                if (type==2)
                {
                  printf("\n**Error** Cannot be both film and image plate**\n");
                  exit(1);
                }
                type = 1;
                break;
             case 'i': /*image plate*/
             case 'I':
                if (type==1)
                {
                  printf("\n**Error** Cannot be both film and image plate**\n");
                  exit(1);
                }
                type = 2;
                break;
             case 's': /*Raster sizes*/
             case 'S':
                if (irast==1)
                {
                   printf("\n**Error** Option -s occurs more than once\n");
                   exit(1);
                }
                if (iarg+2>=argc)
                {
                   printf(
                   "\n**Error** Two raster numbers must follow -s option\n");
                   exit(1);
                }
                iarg++;
                i = sscanf(argv[iarg],"%d",&xrasts);
                if (i==0)
                {
                   printf("**Error** Invalid number follows -s option\n");
                   exit(1);
                }
                iarg++;
                i = sscanf(argv[iarg],"%d",&yrasts);
                if (i==0)
                {
                   printf("**Error** Invalid number follows -s option\n");
                   exit(1);
                }
                irast = 1;
                break;
             case 'c': /*Centre position*/
             case 'C':
                if (cent==1)
                {
                   printf("\n**Error** Option -c occurs more than once\n");
                   exit(1);
                }
                if (iarg+2>=argc)
                {
                   printf(
                   "\n**Error** Two raster values must follow -c option\n");
                   exit(1);
                }
                iarg++;
                i = sscanf(argv[iarg],"%f",&xcen);
                if (i==0)
                {
                   printf("**Error** Invalid number follows -c option\n");
                   exit(1);
                }
                iarg++;
                i = sscanf(argv[iarg],"%f",&ycen);
                if (i==0)
                {
                   printf("**Error** Invalid number follows -c option\n");
                   exit(1);
                }
                cent = 1;
                break;
             case 'd': /*Crystal to film distance*/
             case 'D':
                if (dist==1)
                {
                   printf("\n**Error** Option -d occurs more than once\n");
                   exit(1);
                }
                if (iarg+1>=argc)
                {
                   printf(
                   "\n**Error** Distance value must follow -d option\n");
                   exit(1);
                }
                iarg++;
                i = sscanf(argv[iarg],"%f",&ctof);
                if (i==0)
                {
                   printf("**Error** Invalid number follows -d option\n");
                   exit(1);
                }
                dist = 1;
                break;
             case 'r': /*Rmin value*/
             case 'R':
                if (irmin==1)
                {
                   printf("\n**Error** Option -r occurs more than once\n");
                   exit(1);
                }
                if (iarg+1>=argc)
                {
                   printf(
                   "\n**Error** Radius value must follow -r option\n");
                   exit(1);
                }
                iarg++;
                i = sscanf(argv[iarg],"%f",&rmin);
                if (i==0)
                {
                   printf("**Error** Invalid number follows -r option\n");
                   exit(1);
                }
                irmin = 1;
                break;
             case 'e': /*Echo details of transformation*/
             case 'E':
                echo = 1;
                break;
          }
      }
      else
      {
         ifil = ifil + 1;
         if (ifil>2)
         {
            printf("\n**Error** More than two file names given\n");
            exit(1);
         }
         if (ifil==1)
         {
            strcpy(filin,p);
            fin = open (filin, O_RDONLY, 0);
            if (fin<0)
            {
               printf("\n**Error** Unable to open input file %s\n",filin);
               exit(1);
            }
         }
         else
         {
            strcpy(filout,p);
            fout = creat (filout, 0664);
            if (fout<0)
            {
               printf("\n**Error** Unable to open output file %s\n",filout);
               exit(1);
            }
         }
      }
      iarg++;
   }

   /* Input details not supplied via the parameters */
   /* --------------------------------------------- */

   /* Film/image-plate */

   if (type==0)
   {
      echo = 1;
      labfip:
      printf("Film (f) or image-plate (i) file [f]: ");
      gets(str);
      sscanf(str,"%s",str2);
      switch (str2[0])
      {
         case '\0':
            type = 1;
            break;
         case 'f':
         case 'F':
            type = 1;
            break;
         case 'i':
         case 'I':
            type = 2;
            break;
         default:
            printf("**Error** Invalid file type\n");
            goto labfip;
      }      
   }

   /* Input file name */

   if (ifil==0)
   {
      echo = 1;
      labinf:
      printf("Name of input file: ");
      gets(str);
      i = sscanf(str,"%s",filin);
      if (filin[0]==0) goto labinf;
      fin = open (filin, O_RDONLY, 0);
      if (fin<0)
      {
         printf("**Error** Unable to open input file\n");
         goto labinf;
      }
   }

   /* Output file name */

   if (ifil<2)
   {
      echo = 1;
      laboutf:
      printf("Name of output file: ");
      gets(str);
      i = sscanf(str,"%s",filout);
      if (filout[0]==0) goto laboutf;
      fout = creat (filout, 0664);
      if (fout<0)
      {
         printf("**Error** Unable to open output file\n");
         goto laboutf;
      }
   }
 
   /* Raster sizes */

   if (irast==0)
   {
      echo = 1;
      if (type==1)
      {
         xrasts = 2400;
         yrasts = 2400;
      }
      else
      {
         xrasts = 1200;
         yrasts = 1200;
      }
      labrast:
      printf("No. of x (slow) rasters, no. of y (fast) rasters [%d %d]: ",
             xrasts,yrasts);
      gets(str);
      i = sscanf(str,"%d%d",&j,&k);
      if (i!=-1&&i!=2)
      {
         printf("**Error** Invalid reply\n");
         goto labrast;
      }
      if (i>0) xrasts = j;
      if (i>1) yrasts = k;
   }

   /* Centre position */

   if (cent==0)
   {
      echo = 1;
      xcen = xrasts/2.0;
      ycen = yrasts/2.0;
      labcent:
      printf("Centre x, y in raster units [%.1f %.1f]: ",xcen ,ycen);
      gets(str);
      i = sscanf(str,"%f%f",&x,&y);
      if (i!=-1&&i!=2)
      {
         printf("**Error** Invalid reply\n");
         goto labcent;
      }
      if (i>0) xcen = x;
      if (i>1) ycen = y;
   }

   /* Crystal to film distance */

   if (dist==0)
   {
      echo = 1;
      ctof = xrasts/2.0;
      labdist:
      printf("Crystal to film distance in raster units [%.1f]: ",ctof);
      gets(str);
      i = sscanf(str,"%f",&x);
      if (i!=-1&&i!=1)
      {
         printf("**Error** Invalid reply\n");
         goto labdist;
      }
      if (i>0) ctof = x;
   }

   /* Rmin value */

   if (irmin==0)
   {
      echo = 1;
      rmin = xrasts/20.0;
      labrmin:
      printf("Rmin value in raster units [%.1f]: ",rmin);
      gets(str);
      i = sscanf(str,"%f",&x);
      if (i!=-1&&i!=1)
      {
         printf("**Error** Invalid reply\n");
         goto labrmin;
      }
      if (i>0) rmin = x;
   }

   /* Print details of transformation */       

   if (echo)
   {
      if (type==1)
      {
         printf("\nTransforming film data\n");
      }
      else
      {
         printf("\nTransforming image plate data\n");
      }
      printf("Name of input  file: %s\n",filin);
      printf("Name of output file: %s\n",filout);
      printf("No. x rasters = %d, no. y rasters = %d\n",xrasts,yrasts);
   }

   /* Transform the data */

   xcen = xrasts - xcen;

   if (type==1)
   {
      if (echo) printf("*Reading file*\n");
      nbytes = xrasts*yrasts;
      byte_img = (unsigned char *) malloc(nbytes);
      if (byte_img==0)
      {
         printf("\n**Error** Unable to allocate memory to hold film image\n");
         exit(1);
      }
      i = read_full_film (fin,1,xrasts,yrasts,byte_img);
      if (i<0)
      {
         printf("\n**Error** Unable to allocate temporary buffer\n");
         exit(1);
      }
      if (i>0)
      {
         printf("\n**Error** Error reading record %d\n",i);
         exit(1);
      }
      if (echo) printf("*Gnomonic transformation*\n");
      err = img_maxpix_gnom(type, byte_img, yrasts, xrasts, 1, 1, 1, 1,
                ycen, xcen, ctof, rmin);
      if (echo) printf("*Writing file*\n");
      i = write (fout,byte_img,nbytes);
      if (i!=nbytes)
      {
         printf("\n**Error** Error writing file\n");
         exit(1);
      }
      if (echo) printf("*Tranformation successfully completed*");
      close (fin);
      close (fout);
      free (byte_img);
      exit(0);
   }
   else
   {
      if (echo) printf("*Reading file*\n");
      nbytes = 2*xrasts*yrasts;
      i2_img = (unsigned short *) malloc(nbytes);
      if (i2_img==0)
      {
         printf("\n**Error** Unable to allocate memory to hold image\n");
         exit(1);
      }
      i = read_full_i2 (fin,1,xrasts,yrasts,i2_img);
      if (i<0)
      {
         printf("\n**Error** Unable to allocate temporary buffer\n");
         exit(1);
      }
      if (i>0)
      {
         printf("\n**Error** Error reading record %d\n",i);
         exit(1);
      }
      if (echo) printf("*Gnomonic transformation*\n");
      err = img_maxpix_gnom(type, i2_img, yrasts, xrasts, 1, 1, 1, 1,
                ycen, xcen, ctof, rmin);
      if (echo) printf("*Writing file*\n");
      i = write (fout,i2_img,nbytes);
      if (i!=nbytes)
      {
         printf("\n**Error** Error writing file\n");
         exit(1);
      }
      if (echo) printf("*Transformation successfully completed*");
      close (fin);
      close (fout);
      free (i2_img);
      exit(0);
   }
}
/*
*****************************
**     read_full_film      **
*****************************

Purpose: Read a film image file and return the image data in an array. 
         The film_file may have its data in any of the eight possible
         orders (the default is assumed to  have yf moving fastest and xf 
         slowest) and is in unsigned byte format without any header. The 
         number of x rasters and number of y rasters must be supplied.

Return:  =0 OK; >0 Record number where a read error occurred; -1 unable
         to allocate temporary buffer space.

*/

int read_full_film (fd, order,nxrasts, nyrasts, img_full)

int fd;           /* The file descriptor of the opened film format 
                     file (R)*/
int order;        /* Order of the data in the input film image file
                     in terms of the xf, yf coordinate system as
                     a number from 1 to 8.
                     1   +xf slow   +yf fast
                     2   +xf slow   -yf fast
                     3   -xf slow   +yf fast
                     4   -xf slow   -yf fast
                     5   +yf slow   +xf fast
                     6   +yf slow   -xf fast
                     7   -yf slow   +xf fast
                     8   -yf slow   -xf fast
                     (R)*/
int nxrasts;      /* The number of x rasters (no. of records) (R)*/
int nyrasts;      /* The number of y rasters (pixels in a record) (R)*/
unsigned char img_full[];
                  /* Array returning the read image data; 
                     The slower moving index in
                     the returned data is x going from low x to
                     high x and the faster moving index is y going 
                     from low y to high y. (W)*/
{
   int i, j, k;
   int irec;
   int step;
   int nrecs, nrasts;
   unsigned char *buf;
   int ir, ix;

   /* Initialisations */

   if(order<1||order>8) order=1;
   nrecs = nxrasts;
   nrasts = nyrasts;
   if (order>4)
   {
      nrecs = nyrasts;
      nrasts = nxrasts;
   }

   /* Read data */

   buf = (unsigned char *) malloc(nrasts);
   if (buf==0) return -1;
   irec = 1;
   for (ir=0;ir<nrecs;++ir)
   {
      i = read (fd, buf, nrasts);
      if (i!=nrasts) goto tidy;
      switch (order)
      {
         case 1:
            j = ir*nrasts;
            step = 1;
            break;
         case 2:
            j = (ir+1)*nrasts - 1;
            step = -1;
            break;
         case 3:
            j = (nrecs-ir-1)*nrasts;
            step = 1;
            break;
         case 4:
            j = (nrecs-ir)*nrasts - 1;
            step = -1;
            break;
         case 5:
            j = ir;
            step = nrecs;
            break;
         case 6:
            j = (nrasts-1)*nrecs + ir;
            step = -nrecs;
            break;
         case 7:
            j = nrecs - 1 - ir;
            step = nrecs;
            break;
         case 8:
            j = nrasts*nrecs - 1 - ir;
            step = -nrecs;
            break;
      }
      for (k=0;k<nrasts;++k)
      {
          img_full[j] = buf[k];
          j += step;
      }
      irec++;
   }
   irec = 0;   
   
   /* Tidy up */

   tidy:
   free (buf);
   return irec;
}
/*
*****************************
**     read_full_i2        **
*****************************

Purpose: Read an image-plate 'i2' image file and return the image data in 
         an array. The file may have its data in any of the eight
         possible orders and may be byte swapped if needed (the default
         is assumed to  have yf moving fastest and xf 
         slowest) and is in unsigned two byte integer format without any 
         header. The number of x rasters and number of y rasters must be 
         supplied.

Return:  =0 OK; >0 Record number where a read error occurred.

*/
int read_full_i2 (fd, order, nxrasts, nyrasts, img_full)

int fd;           /* The file descriptor of the opened film format 
                     file (R)*/
int order;        /* Order of the data in the input film image file
                     in terms of the xf, yf coordinate system as
                     a number from 1 to 8. (-order can be given to
                     byte swap the data as well as change the order)
                     1   +xf slow   +yf fast
                     2   +xf slow   -yf fast
                     3   -xf slow   +yf fast
                     4   -xf slow   -yf fast
                     5   +yf slow   +xf fast
                     6   +yf slow   -xf fast
                     7   -yf slow   +xf fast
                     8   -yf slow   -xf fast
                     (R)*/
int nxrasts;      /* The number of x rasters (no. of records) (R)*/
int nyrasts;      /* The number of y rasters (pixels in a record) (R)*/
unsigned short img_full[];
                  /* Array returning the read image data; 
                     The slower moving index in
                     the returned data is x going from low x to
                     high x and the faster moving index is y going 
                     from low y to high y i.e. the same order as in the
                     file. (W)*/
{
   int i, j, k, l;
   int swap;
   int irec;
   int step;
   int nrecs, nrasts;
   unsigned short *buf;
   unsigned char *cbuf;
   unsigned char tempc;
   int ir, ix;

   /* Initialisations */

   swap = 0;
   if (order<0)
   {
      swap = 1;
      order = -order;
   }
   if(order<1||order>8) order=1;
   nrecs = nxrasts;
   nrasts = nyrasts;
   if (order>4)
   {
      nrecs = nyrasts;
      nrasts = nxrasts;
   }

   /* Read data */

   buf = (unsigned short *) malloc(nrasts*sizeof(unsigned short));
   if (buf==0) return -1;
   cbuf = (unsigned char *) buf;
   irec = 1;
   for (ir=0;ir<nrecs;++ir)
   {
      i = read (fd, buf, 2*nrasts);
      if (i!=2*nrasts) goto tidy;
      switch (order)
      {
         case 1:
            j = ir*nrasts;
            step = 1;
            break;
         case 2:
            j = (ir+1)*nrasts - 1;
            step = -1;
            break;
         case 3:
            j = (nrecs-ir-1)*nrasts;
            step = 1;
            break;
         case 4:
            j = (nrecs-ir)*nrasts - 1;
            step = -1;
            break;
         case 5:
            j = ir;
            step = nrecs;
            break;
         case 6:
            j = (nrasts-1)*nrecs + ir;
            step = -nrecs;
            break;
         case 7:
            j = nrecs - 1 - ir;
            step = nrecs;
            break;
         case 8:
            j = nrasts*nrecs - 1 - ir;
            step = -nrecs;
            break;
      }
      if (swap)
      {
         l=0;
         for (k=0;k<nrasts;++k)
         {
             tempc = cbuf[l];
             cbuf[l] = cbuf[l+1];
             cbuf[l+1] = tempc;
             l += 2;
         }
      }
      for (k=0;k<nrasts;++k)
      {
          img_full[j] = buf[k];
          j += step;
      }
      irec++;
   }
   irec = 0;   
   
   /* Tidy up */

   tidy:
   free (buf);
   return irec;
}
/*
*****************************
**     img_maxpix_gnom     **
*****************************

PURPOSE: Convert an image to its gnomonic projection using 'maximum pixel'
         method

RETURN:  Error flag =0 no error
                    =1 Unable to allocate memory for image work array
                    =2 unable to allocate memory for ratios work array

AUTHOR:  John W. Campbell, July 1991

NOTE:    The fast moving axis in the input image is designated x and
         the slower moving index as y.

*/

int img_maxpix_gnom (type, imgdata, img_width, img_height, xo_pix, yo_pix, 
                nxp, nyp, x_cen, y_cen, cfd, rmin)

int type;                /*Type 1=film data, 2=image-plate (i2) data*/
unsigned char *imgdata;  /*Image data array (x fastest) (M)*/
int img_width;           /*Width of image in pixels in 'imgdata' (R)*/
int img_height;          /*Height of image in pixels im 'imgdata' (R)*/
int xo_pix;              /*The x origin pixel number (pixels numbered from 1 up)
                           for the section of the film image passed to the 
                           routine in 'imgdata' (R)*/
int yo_pix;              /*The y origin pixel number for the section of the
                           film image passed to the routine in 'imgdata' (R)*/
int nxp;                 /*The no. of x pixels from the original image
                           compressed into one in 'imgdata' (R)*/
int nyp;                 /*The no. of y pixels from the original image
                           compressed into one in 'imgdata' (R)*/
float x_cen;             /*x centre pixel position (uncompressed image) (R)*/
float y_cen;             /*y centre pixel position (uncompressed image) (R)*/
float cfd;               /*Crystal to film distance in (y) raster units of the
                           uncompressed image(R)*/
float rmin;              /*Minimum radius in (y) raster units of the 
                           uncompressed image (R)*/

{
   float cgd;        /*Crystal to gnomonic projection plane distance in
                       raster units*/
   unsigned short *i2data;    /*Pointer to image data array for use with i2
                                data*/
   unsigned char *img_temp;   /*Array to hold data during transformation*/
   unsigned short *i2_temp;   /*As above but for i2 image*/
   float *ratio;     /*Array to hold pixel expansion ratios (approx) as a 
                       function of two theta*/
   float xc;         /*x centre in raster units in compressed image*/
   float yc;         /*y centre in raster units in compressed image*/
   float xmax;       /*width of compressed image in raster units*/
   float ymax;       /*height of compressed image in raster units*/
   float rmax;       /*Maximum radius from input image size and centre*/
   float piov2;      /*Pi divided by 2.0*/
   float th2min;     /*Minimum two theta value*/
   float th2max;     /*Maximum two theta value*/
   float thg;        /*Theta for gnomonic projection*/
   float dth;        /*Theta interval for expansion ratios table*/
   float dthov2;     /*Half of dth*/
   float dx;         /*Small distance used in calculating expansion ratios*/
   float x1;         /*Temp x coordinate in raster units*/
   float x2;         /*Temp x coordinate in raster units*/
   float y1;         /*Temp y coordinate in raster units*/
   float y2;         /*Temp y coordinate in raster units*/
   float xp;         /*Untransformed x of pixel centre in raster units for
                       current pixel*/
   float yp;         /*Untransformed y of pixel centre in raster units for
                       current pixel*/
   float r1;         /*Distance of untransformed pixel from centre in raster
                       units*/
   float r2;         /*Distance of transformed pixel from centre in raster
                       units*/
   float th2;        /*Two theta value (untransformed pixel)*/
   float halfsiz;    /*Half size (each dimension) of transformed pixel area*/
   float area;       /*Transformed pixel area*/
   float xt;         /*Transformed pixel area centre x in raster units*/
   float yt;         /*Transformed pixel area centre y in raster units*/
   float xtmin;      /*x start of transformed pixel area in raster units*/
   float xtmax;      /*x end of transformed pixel area in raster units*/
   float ytmin;      /*y start of transformed pixel area in raster units*/
   float ytmax;      /*y end of transformed pixel area in raster units*/
   float xpmin;      /*x start of current pixel in transformed area*/
   float xpmax;      /*x end of current pixel in transformed area*/
   float ypmin;      /*y start of current pixel in transformed area*/
   float ypmax;      /*y end of current pixel in transformed area*/
   int i;            /*Temp/loop variable*/
   int ip;           /*Offset into image data arrays*/
   int ih;           /*Loop variable over y pixels*/
   int ihh;          /*Inner loop variable over y pixels*/
   int iw;           /*Loop variable over x pixels*/
   int iww;          /*Inner loop variable over x pixels*/
   int xpixmin;      /*Minimum x pixel which transformed pixel area may 
                       overlap*/
   int xpixmax;      /*Maximum x pixel which transformed pixel area may 
                       overlap*/
   int ypixmin;      /*Minimum y pixel which transformed pixel area may 
                       overlap*/
   int ypixmax;      /*Maximum y pixel which transformed pixel area may 
                       overlap*/

   /* Set Pi divided by 2 */

   piov2 = 1.5707963;

   /* Other initialisations */

   i2data = (unsigned short *) imgdata;

   /*Allocate memory for work arrays*/

   if (type==1)
   {
      img_temp = (unsigned char *) malloc (img_width*img_height*sizeof(char));
   }
   else
   {
      i2_temp = (unsigned short *) 
                malloc (img_width*img_height*sizeof(unsigned short));
   }
   if (img_temp==0) return 1;

   ratio = (float *) malloc (201*sizeof(float));
   if (ratio==0) return 2;

   /* Centre and radius calculation */

   cfd /= (float) nyp;
   rmin /= (float) nyp;
   xc = (float) (x_cen-xo_pix+1) / (float) nxp - 0.5 / (float) nxp;
   yc = ((float) (img_height*nyp) + (float) yo_pix - (float) y_cen) 
        / (float) nyp - 0.5 / (float) nyp;
   xmax = (float) img_width;
   ymax = (float) img_height;
   rmax = xc;
   if ((xmax-xc)<rmax) rmax = xmax - xc;
   if (yc<rmax) rmax = yc;
   if ((ymax-yc)<rmax) rmax = ymax - yc;

   /* Calculate two theta limits */

   th2max = (float) atan ((double)(rmax/cfd));
   th2min = (float) atan ((double)(rmin/cfd));

   /* Calculate gnomonic projection distance*/

   thg = piov2 - th2min / 2.0;
   cgd = rmax / (float) tan ((double)thg);

   /* Set up pixel expansion ratios (as a function of two theta )table */

   dth = (th2max - th2min) / 200.0;
   dthov2 = dth / 2.0;
   for (i=0;i<=200;++i)
   {
      th2 = th2min + (float) i * dth;
      dx = cfd * ((float) tan((double)(th2+dthov2)) 
           - (float) tan((double)(th2-dthov2)));
      thg = piov2 - (th2-dthov2) / 2.0;
      x1 = cgd * (float) tan ((double)thg);
      thg = piov2 - (th2+dthov2) / 2.0;
      x2 = cgd * (float) tan ((double)thg);
      ratio[i] = (x1 - x2) / dx;
   }

   /* Initialise array to hold transformed data */

   ip=0;
   for (ih=0;ih<img_height;++ih)
   {
      for (iw=0;iw<img_width;++iw)
      {
          if (type==1){ img_temp[ip] = 0; }else{ i2_temp[ip] = 0;}
          ip++;
      }
   }

   /* Do gnomonic transformation now into temporary array */

   ip=0;
   for (ih=0;ih<img_height;++ih)
   {
      yp = (float) ih + 0.5;
      for (iw=0;iw<img_width;++iw)
      {
         xp = (float) iw + 0.5;
         r1 = (float) sqrt((double) ((xp-xc)*(xp-xc) + (yp-yc)*(yp-yc)));
         th2 = (float) atan ((double) (r1 / cfd));
         if (th2<th2min||th2>th2max){ip++; continue;}
         i = (int) 200.0 * ((th2 - th2min) / (th2max - th2min));
         halfsiz = ratio[i]/2.0;
         thg = piov2 - th2 / 2.0;
         r2 = cgd * (float) tan ((double) thg);
         xt = xc + (xp - xc) * (r2 / r1);
         xtmin = xt - halfsiz;
         xtmax = xt + halfsiz;
         yt = yc - (yp - yc) * (r2 / r1);
         ytmin = yt - halfsiz;
         ytmax = yt + halfsiz;
         area = 4.0 * halfsiz * halfsiz;
         xpixmin = (int) xtmin + 1;
         if (xpixmin<0) xpixmin = 0;
         xpixmax = (int) (xtmax-0.00001) + 1;
         if (xpixmax>=img_width) xpixmax = img_width - 1;
         ypixmin = (int) ytmin + 1;
         if (ypixmin<0) ypixmin = 0;
         ypixmax = (int) (ytmax-0.00001) + 1;
         if (ypixmax>=img_height) ypixmax = img_height - 1;
         if (type==1)
         {
            img_temp[ip]=0;
            for (ihh=ypixmin;ihh<=ypixmax;++ihh)
            {
               for (iww=xpixmin;iww<=xpixmax;++iww)
               {
                  i = (ihh-1)*img_width + iww -1;
                  if (imgdata[i]>img_temp[ip]) img_temp[ip]=imgdata[i];
               } 
            }
         }
         else
         {
            i2_temp[ip]=0;
            for (ihh=ypixmin;ihh<=ypixmax;++ihh)
            {
               for (iww=xpixmin;iww<=xpixmax;++iww)
               {
                  i = (ihh-1)*img_width + iww -1;
                  if (i2data[i]>i2_temp[ip]) i2_temp[ip]=i2data[i];
               } 
            }
         }
         ip++;
      }
   }

   /* Rewrite image data array with transformed data */

   if (type==1)
   {
      ip=0;
      for (ih=0;ih<img_height;++ih)
      {
         for (iw=0;iw<img_width;++iw)
         {
             imgdata[ip] = img_temp[ip];
             ip++;
         }
      }
      free (img_temp);
   }
   else
   {
      ip=0;
      for (ih=0;ih<img_height;++ih)
      {
         for (iw=0;iw<img_width;++iw)
         {
             i2data[ip] = i2_temp[ip];
             ip++;
         }
      }
      free (i2_temp);
   }
   free (ratio);
   return 0;
}