This is a cached copy of http://cfa-www.harvard.edu/~jhora/mirac/xcor.tar.Z
	


/*************************************************************************/
/* xcor.c - Determine offsets for image based on cross correlation       */
/*************************************************************************/
/*     Febuary 22, 1993                                                  */
/*     Joseph L Hora    adapted from the CROSSR fortran program written  */
/*                      by R. Tresch-Feinberg and J. Hora                */
/*     requires the files xfitsio.c and xfiocom.h to compile             */
/*************************************************************************/
/*                                                                       */
/* xcor                                                                  */
/*                                                                       */
/*                                                                       */
/* This program reads in image files in FITS format (8, 16, 32, or -32   */
/* BITPIX) and calculates the cross correlation between two images.  At  */
/* least two images must be entered, the reference image (whose offsets  */
/* are by definition zero) and the search or window image.  A small      */
/* window out of the search image is compared to part or all of the      */
/* reference image.  The offset result given by the program is the offset*/
/* in pixels necessary to align the search image with the reference      */
/* image.  The program is currently configured for a maximum image size  */
/* of 256x256, but that can be easily changed using the definitions in   */
/* the file "fiocom.h" for the paramters XPIX and YPIX.                  */
/*                                                                       */
/* The program writes the output to the screen, and optionally to an     */
/* output file.  The following information is in the output file for each*/
/* image:  the image name, X and Y coord. of the "center" structure, the */
/* X and Y offsets necessary to align the image with the reference image,*/
/* and the max, min, and average correlation values.                     */
/*                                                                       */
/*                                                                       */
/* The program can be compiled on a SUN workstation using the following  */
/* command:                                                              */
/*                                                                       */
/* cc xcor.c xfitsio.c -lm -o xcor                                       */
/*                                                                       */
/* The executable file "xcor" will be created.                           */
/*                                                                       */
/*  see message in program below for info on command line parameters     */
/*************************************************************************/
/*                                                                       */
/*  MODIFICATION LOG                                                     */
/*                                                                       */
/*   1/25/95 - JLH                                                       */
/*    Changed the smoothing routines so that they use global arrays,     */
/*    and now only necessary part of window image is smoothed.           */
/*                                                                       */
/*  10/20/93 - JLH                                                       */
/*    Changed program in xfitsio.c and xfiocom.h - the data files are    */
/*    read in as floating point, then scaled if necessary and converted  */
/*    to long integer.  The scaling now also handles small numbers and   */
/*    large, so that rounding errors don't adversely affect the cross-   */
/*    correlation.                                                       */
/*                                                                       */
/*  10/14/93 - JLH                                                       */
/*    Corrected error in -q option: when specifying coordinates, it was  */
/*    getting them switched.                                             */
/*                                                                       */
/*************************************************************************/

#include <stdio.h>
#include <string.h>
#include <stdlib.h>
#include <math.h>
#include "xfiocom.h"

#define NULL 0
#define MAXES 6
#define TRUE 1
#define FALSE 0


static long refimg[XPIX][YPIX];   /* this is where final image is stored */

static long winimg[XPIX][YPIX];     /* mask file for a single image */
static long tmpimg[XPIX][YPIX];     /* mask file for a single image */
static float rdimg[XPIX][YPIX];     /* mask file for a single image */

static long dmax;                 /* max. good data value */

static float maxvalue;
static int mwin,lwin;                 /* number of frames used in final image*/
static char speedy, xcorrelate, do_centroid;

void save_final_image();
void do_xcorr();
void winsmooth();
void refsmooth();
int max();
int min();


main(argc, argv)           /***********Main XCOR program start**********/
int argc;
char *argv[];
{
   char  fname[50];
   float xtmp,ytmp,cmax,cmin,cavg,usemax;
   char  fits_input,expand,discard_file, do_centroid, write_individual_files,
          outfileopen,interpolate,peak_coords,wsmooth,rsmooth,peakset;
   int   ii,
         imin,
         imax,
         jmin,
         jmax,
         refpkx,refpky,
         nbuf,
         itmin,
         itmax,
         bitpix,
         jtmin,
         jtmax,
         tnum,
         rsmfac,wsmfac,
         xst,yst,ipk,jpk,xreg,yreg,xsim,ysim,
         hfsize,  rxoff,ryoff,rxend,ryend,wxoff,wyoff,wxend,wyend,wxsiz,wysiz,
         rxsiz,rysiz,
         fmin,fmax;
   int numexclude;
   struct BUFFER searchead, windhead;

   char  object[8];
   char wname[60],excludefile[60], peakfile[60], fits_name[60],refname[60],
        winname[60],outfile[60],infile[60];

   char lbuf[80];

   int   framenum,
         nump,
         st_arg;
   int   n_rawimg;
   int   i,j, k, l;
   char tstr[20];
   FILE *fp1,*fp2;

/*  Initialize parameter values  */

   strcpy(peakfile,"NONE");
   strcpy(fits_name,"NONE");
   strcpy(infile,"NONE");
   strcpy(outfile,"NONE");
   imin=-1;
   imax=XPIX-1;
   jmin=0;
   fmin = fmax = -1;
   rsmfac=4;
   wsmfac=4;
   dmax=216000000;
   mwin=20;
   lwin=20;
   do_centroid=FALSE;
   xcorrelate=TRUE;
   jmax=YPIX-1;
   hfsize = 15;
   refpkx = refpky = -1;
   peakset = outfileopen=FALSE;
   st_arg = 1;
   do_centroid=FALSE;
   write_individual_files = FALSE;
   ipk = -999;
   fits_input = FALSE;
   wsmooth = FALSE;
   rsmooth = FALSE;
   peak_coords = FALSE;
   expand = FALSE;
   interpolate = FALSE;
   speedy = FALSE;

   printf("Cross correlation program ( xcor.c )\n\n");

/*    If enough command line parameters have been entered, read them in    */
   if (argc >=3)
   {
      /* Fits filename  of reference image is 1st argument */

      sscanf(argv[1],"%s",refname);

      sscanf(argv[2],"%s",fname);
      if (strcspn(fname,"-")==0)
      {
         st_arg = 2;
         strcpy(fname,"NONE");
      }
      else
      {
         st_arg = 3;
         strcpy(wname,fname);
      }

      if (argc>st_arg)
      {
         for (i= st_arg; i<argc; i++)
         {
            sscanf(argv[i],"%s",tstr);
            if (strcspn(tstr,"-")==0)
            {
               switch (tstr[1])
               {
                  case 's':switch (tstr[2])
                           {
                              case 'r':
                                 rsmooth = TRUE;
                                 if (i<argc-1)
                                 {
                                   sscanf(argv[i+1],"%s",tstr);
                                   if (strcspn(tstr,"-")!=0)
                                      nump = sscanf(argv[++i],"%i",&rsmfac);
                                 }
                                 break;

                              case 'w':
                                 wsmooth = TRUE;
                                 if (i<argc-1)
                                 {
                                   sscanf(argv[i+1],"%s",tstr);
                                   if (strcspn(tstr,"-")!=0)
                                      nump = sscanf(argv[++i],"%i",&wsmfac);
                                 }
                                 break;
                           }
                           break;

                  case 'c':do_centroid = TRUE;  /* Turn centroiding on */
                           break;

                  case 'q':speedy = TRUE;
                           if (i<argc-1)
                           {
                              sscanf(argv[i+1],"%s",tstr);
                              if (strcspn(tstr,"-")!=0)
                              {
                                 nump = sscanf(argv[++i],"%i",&refpkx);
                                 nump = sscanf(argv[++i],"%i",&refpky);
                                 if (nump==0)
                                 {
                                    printf("Error specifying reference peak\n");
                                    exit(1);
                                 }
                                 sscanf(argv[i+1],"%s",tstr);
                                 if (strcspn(tstr,"-")!=0)
                                 {
                                    nump = sscanf(argv[++i],"%i",&ipk);
                                    nump = sscanf(argv[++i],"%i",&jpk);
                                    if (nump==0)
                                    {
                                       printf("Error specifying reference peak\n");
                                       exit(1);
                                    }
                                    peakset=TRUE;
                                 }
                              }
                           }
                           break;

                  case 'r':sscanf(argv[++i],"%i",&fmin);
                           nump= sscanf(argv[++i],"%i",&fmax);
                           if (nump==0)
                           {
                              printf("Error specifying file limits\n");
                              exit(1);
                           }
                           break;


                                     /* user has specified search limits */
                  case 'l':sscanf(argv[++i],"%i",&imin);
                           sscanf(argv[++i],"%i",&imax);
                           sscanf(argv[++i],"%i",&jmin);
                           nump = sscanf(argv[++i],"%i",&jmax);
                           if (nump==0)
                           {
                              printf("Error specifying limits\n");
                              exit(1);
                           }
                           imin=max(0,min(XPIX-1,imin));
                           imax=max(0,min(XPIX-1,imax));
                           jmin=max(0,min(YPIX-1,jmin));
                           jmax=max(0,min(YPIX-1,jmax));
                           break;
                                        /* enter the search window size */
                  case 'b':if (i<argc+2)
                           {
                              sscanf(argv[++i],"%i",&mwin);
                              nump= sscanf(argv[++i],"%i",&lwin);
                              if (nump==0)
                              {
                                 printf("Error entering search window size\n");
                                 exit(1);
                              }
                           }
                           else
                           {
                              printf("Error entering search window size\n");
                              exit(1);
                           }
                           break;

                  case 'd':sscanf(argv[++i],"%d",&dmax);
                           break;


                  case 'o':sscanf(argv[++i],"%s",outfile);
                           break;

                  case 'i':sscanf(argv[++i],"%s",infile);
                           break;

                  case 'p':sscanf(argv[++i],"%s",peakfile);
                           if (strcspn(peakfile,"-")==0)
                           {
                              printf("\nError specifying peak coordinate file\n");
                              exit(1);
                           }
                           peak_coords=TRUE;
                           break;

               }
            }
         }
      }
   }
   else
   {
            /* Display general info on the program */
printf("\nusage: xcor referenceimg [windowimg] [-p peakfile][-sr (sbox)][-sw (sbox)]\n");
printf("            [-b xsize ysize] [-d datamax] [-q (xpk ypk)] [-i inputlist]\n");
printf("            [-o outputlist] [-r min# max#] [-l imin imax jmin jmax]\n");
printf("  where\n");
printf("        referenceimg is the FITS file of the image to compare all the search\n");
printf("              images to.\n");
printf("        windowimg is the FITS image to cut the window from to compare it to\n");
printf("              the reference image.  The window image does not need to be\n");
printf("              entered if an input list is entered with the -i or -p option.\n");
printf("       -l  specifies the searchbox pixel limits in the search image. If\n");
printf("              this switch is not set, then the searchbox pixel limits are \n");
printf("              automatically calculated, using the peak value in the image.\n");
printf("       -b  sets the size of the searchbox.  Default is 20 pixels in X & Y.\n");
printf("       -d  specifies the maximum Data value to accept as a valid peak.\n");
printf("              This helps in eliminating bad pixels.  Default is 2^32/2.\n");
printf("              (the correlation calculation is done using long integers)\n");
printf("       -i  specify input file list; files are listed one per line, full name\n");
printf("              If this option is selected, it overrides any windowimg name\n");
printf("              entered in the command line.\n");
printf("       -o  specify output file for results.  The following information is\n");
printf("           printed: file name (or #), X and Y of the matching structure,\n");
printf("           X and Y shift distance, correlation max, min, and avg.\n");
printf("\n");
printf(" Type return to continue...");
interpolate = getchar();
printf("\n\n");
printf("       -sr smooth the reference image.  Smooth factor can be set, default=4\n");
printf("       -sw smooth the window image(s).  Smooth factor can be set, default=4\n");
printf("       -q  Quick mode: only search region within searchbox of reference \n");
printf("           position.  Optional reference position may be entered, otherwise\n");
printf("           peak value is used\n");
printf("       -p  Peak coordinate list option.  The search window will be centered\n");
printf("           on the position given for each file in the input list.  This option\n");
printf("           superseeds the windowimg specification and the -i option.  The\n");
printf("           format of the peak file is one line per file, with the file name\n");
printf("           followed by the X and Y coordinate of the peak position.  This is\n");
printf("           the format of the file output by the IRAF script peak_eval \n");
printf("       -r  range of files: the min and max of range is entered, then a three\n");
printf("           number extension is added on to the windowimg root name to read files\n");
printf("\n");
printf("  Example 1:  Single file\n");
printf("     xcor ref1.fts data1.fts -sr -sw -q -o xcorpeaks.txt\n");
printf("  Example 2:  Use list of object peak positions\n");
printf("     xcor f921210a.210 -p pk200 -sr -sw -q -o objoffsets \n");
         exit(0);
   }

   for (i=0; i<XPIX; i++)             /* initialize image buffers */
      for (j=0; j<YPIX; j++)
         refimg[i][j] = winimg[i][j] = 0;
   printf("Now reading reference image...\n");
   if (speedy && (refpkx>0))
   {
      rxoff = max(refpkx-2*lwin-rsmfac,0);
      rxend = min(refpkx+2*lwin+rsmfac,XPIX+rxoff-1);
      ryoff = max(refpky-2*mwin-wsmfac,0);
      ryend = min(refpky+2*mwin+wsmfac,YPIX+ryoff-1);
      if (readfits(refname,&searchead,rdimg,&rxoff,&rxend,&ryoff,&ryend)!=0)
      {
         printf("Program terminating.\n");
         exit(1);
      }
      refpkx -= rxoff;
      refpky -= ryoff;
   }
   else
   {
      rxoff = ryoff = 0;
      rxend = XPIX;    ryend = YPIX;
      if (readfits(refname,&searchead,rdimg,&rxoff,&rxend,&ryoff,&ryend)!=0)
      {
         printf("Program terminating.\n");
         exit(1);
      }
   }
   rxsiz = rxend - rxoff + 1;
   rysiz = ryend - ryoff + 1;
   bitpix = searchead.bitpix;

   if ((rxsiz<refpkx) || (rysiz<refpky))
   {
      printf("reference point is outside map!\n");
      exit(1);
   }
   printf("File %s header read, parameters: %i %i %i\n",refname,rxsiz,rysiz,
           searchead.bitpix);
   usemax = fabs(max( searchead.max, -searchead.min));
/*   Image may have to be scaled so that the I*4 calculations do not overflow */
   if ((usemax>100000)||(usemax<1000))
   {
      printf("Now scaling reference image...\n");
      for (i=0; i<rxsiz; i++)
         for (j=0; j<rysiz; j++)
            refimg[i][j] = (long) (rdimg[i][j]/usemax*100000.0);

   }
   else
   {
      printf("Now loading ref image into array...\n");
      for (i=0; i<rxsiz; i++)
         for (j=0; j<rysiz; j++)
            refimg[i][j] = (long)(rdimg[i][j]);
   }

   if (rsmooth)
      refsmooth(rsmfac,rxsiz,rysiz);

/*   If search window specified for reference image, find center and load into */
/*   the proper variables                                                      */

   tnum = -21600000;

/*   Find peak position in reference image if necessary                        */
   if ((speedy) && (refpkx==-1))
   {
      for (i=0; i<rxsiz; i++)
         for (j=0; j<rysiz; j++)
            if ((refimg[i][j]>tnum)&&(refimg[i][j]<dmax))
             {
                tnum = refimg[i][j];
                refpkx = i;
                refpky = j;
             }
      printf("Reference peak found: %i at %i,%i\n",tnum,refpkx+1,refpky+1);
   }
   else if (refpkx!=-1)
      printf("Using input reference point: %i,%i\n",refpkx+rxoff,refpky+ryoff);

/*   If results are being written to an output file, open the file here        */
   if (strcspn(outfile,"NONE")!=0)
   {
      if ((fp2 = fopen(outfile,"w"))==NULL)
      {
         printf("Error opening output file: %s\n",outfile);
      }
      else
         outfileopen=TRUE;
   }

/*   Processing loops - for each of the subsections below, files are read in,  */
/*   scaled if necessary and smoothed if that option is chosen.  Then the      */
/*   offsets are calculated and the results displayed to the screen and written*/
/*   to the output file if open.                                               */

   if (fmin>-1)   /*  The user has entered a range of files with the -r option */
   {
      for (i=fmin; i<=fmax; i++)
      {
         strcpy(wname,fname);   /* construct file names by adding extension to */
         if (i==0)              /* the root name entered on the command line   */
            strcat(wname,"00");
         else if (i<10)
            strcat(wname,"00");
         else if (i<100)
            strcat(wname,"0");
         sprintf(tstr,"%i",i);
         strcat(wname,tstr);
         wxoff = 0;
         wxend = XPIX;
         wyoff = 0;
         wyend = YPIX;
         if(readfits(wname,&windhead,rdimg,&wxoff,&wxend,&wyoff,&wyend)!=0)
         {
            printf("Program terminating.\n");
            exit(1);
         }
         wxsiz = wxend - wxoff + 1;
         wysiz = wyend - wyoff + 1;
         printf("File %s header read, parameters: %i %i %i\n",wname,wxsiz,wysiz,
           windhead.bitpix);
         if ((windhead.max>100000)||(windhead.max<100))
         {
            printf("Now rescaling search image...\n");
            for (ii=0; ii<wxsiz; ii++)
               for (j=0; j<wysiz; j++)
                  winimg[ii][j] = (long) (((double)rdimg[ii][j])/
                                        ((double)windhead.max)*100000.0);
         }
         else
         {
            for (ii=0; ii<wxsiz; ii++)
               for (j=0; j<wysiz; j++)
                  winimg[ii][j] = (long)(rdimg[ii][j]);
         }
         if (!peakset)
         {
            ipk = -999;
            jpk = -999;
         }

         do_xcorr(rxsiz,rysiz,wxsiz,wysiz,lwin,mwin,imin, imax, jmin, jmax,
              ipk,jpk,refpkx,refpky,hfsize,&xreg,&yreg,&cmax,&cmin,&cavg,
              wsmooth,wsmfac);

         xreg += rxoff - wxoff;
         yreg += ryoff - wyoff;

         printf(" Results:  %5i %5i %10.5f %10.5f %10.5f\n\n",xreg,yreg,
            cmax,cmin,cavg);
         if (outfileopen)
           fprintf(fp2,"%s %3i %3i %4i %4i %3i %3i %10.5f %10.5f %10.5f\n",
            wname,refpkx-xreg+1,refpky-yreg+1,windhead.naxis1,windhead.naxis2,
            xreg,yreg,cmax,cmin,cavg);

      }
   }
   else if (peak_coords)  /*  The user has entered a file that is a list of    */
                          /*  image file names, with the peak positions given  */
   {
      printf("Now reading file from peak coordinate file %s\n",peakfile);
      if ((fp1 = fopen(peakfile,"r"))==NULL)
      {
         printf("Error opening input file: %s\n",peakfile);
         exit(1);
      }
      while (!feof(fp1))     /*  Read in the image file name and peak coords   */
      {
         if (fgets(lbuf,80,fp1)!=NULL)
         {
            sscanf(lbuf,"%s%f%f",wname,&xtmp,&ytmp);
            ipk = (int)(xtmp);
            jpk = (int)(ytmp);
            printf("Peak read: %s %i %i\n",wname,ipk,jpk);

            if (speedy)
            {
               wxoff = ipk - lwin;
               wxend = ipk + lwin + 1;
               wyoff = jpk - mwin;
               wyend = jpk + mwin + 1;
               if(readfits(wname,&windhead,rdimg,&wxoff,&wxend,&wyoff,&wyend)!=0)
               {
                  printf("Program terminating.\n");
                  exit(1);
               }
            }
            else
            {
               wxoff = 0;
               wxend = XPIX;
               wyoff = 0;
               wyend = YPIX;
               if(readfits(wname,&windhead,rdimg,&wxoff,&wxend,&wyoff,&wyend)!=0)
               {
                  printf("Program terminating.\n");
                  exit(1);
               }
            }
            ipk -= wxoff;
            jpk -= wyoff;
            wxsiz = wxend - wxoff + 1;
            wysiz = wyend - wyoff + 1;
            printf("File %s header read, parameters: %i %i %i\n",wname,
              wxsiz,wysiz,windhead.bitpix);

            if ((windhead.max>100000)||(windhead.max<100))
            {
               printf("Now rescaling search image...\n");
               for (i=0; i<wxsiz; i++)
                  for (j=0; j<wysiz; j++)
                     winimg[i][j] = (long) (((double)rdimg[i][j])/
                                        ((double)windhead.max)*100000.0);
            }
            else
            {
               for (i=0; i<wxsiz; i++)
                  for (j=0; j<wysiz; j++)
                     winimg[i][j] = (long)(rdimg[i][j]);
            }

            do_xcorr(rxsiz,rysiz,wxsiz,wysiz,lwin,mwin,imin, imax, jmin, jmax,
                 ipk,jpk,refpkx,refpky,hfsize,&xreg,&yreg,&cmax,&cmin,&cavg,
                 wsmooth,wsmfac);
            xreg += rxoff - wxoff;
            yreg += ryoff - wyoff;
            printf(" %5i %5i %10.5f %10.5f %10.5f\n",xreg,yreg,cmax,cmin,cavg);

            if (outfileopen)
           fprintf(fp2,"%s %3i %3i %4i %4i %3i %3i %10.5f %10.5f %10.5f\n",
            wname,refpkx-xreg+1,refpky-yreg+1,windhead.naxis1,windhead.naxis2,
            xreg,yreg,cmax,cmin,cavg);

         }
      }
      fclose(fp1);
   }
   else if (strcspn(infile,"NONE")!=0)  /*  The user has input a file that has */
                                        /*  a list of image file names         */
   {
      printf("Now reading file from input file %s\n",infile);;
      if ((fp1 = fopen(infile,"r"))==NULL)
      {
         printf("Error opening input file: %s\n",infile);
         exit(1);
      }
      while (!feof(fp1))
      {
         if (fgets(lbuf,80,fp1)!=NULL)
         {
            sscanf(lbuf,"%s",wname);
            wxoff = 0;
            wxend = XPIX;
            wyoff = 0;
            wyend = YPIX;
            if(readfits(wname,&windhead,rdimg,&wxoff,&wxend,&wyoff,&wyend)!=0)
            {
               printf("Program terminating.\n");
               exit(1);
            }
            wxsiz = wxend - wxoff + 1;
            wysiz = wyend - wyoff + 1;

            printf("File %s header read, parameters: %i %i %i\n",wname,
              wxsiz,wysiz,windhead.bitpix);

            if ((windhead.max>100000)||(windhead.max<100))
            {
               printf("Now rescaling search image...\n");
               for (i=0; i<wxsiz; i++)
                  for (j=0; j<wysiz; j++)
                     winimg[i][j] = (long) (((double)rdimg[i][j])/
                                        ((double)windhead.max)*100000.0);
            }
            else
            {
               for (i=0; i<wxsiz; i++)
                  for (j=0; j<wysiz; j++)
                     winimg[i][j] = (long)(rdimg[i][j]);
            }
            ipk = -999;

            do_xcorr(rxsiz,rysiz,wxsiz,wysiz,lwin,mwin,imin, imax, jmin, jmax,
                 ipk,jpk,refpkx,refpky,hfsize,&xreg,&yreg,&cmax,&cmin,&cavg,
                 wsmooth,wsmfac);
            xreg += rxoff - wxoff;
            yreg += ryoff - wyoff;
            printf(" %5i %5i %10.5f %10.5f %10.5f\n",xreg,yreg,cmax,cmin,cavg);

            if (outfileopen)
           fprintf(fp2,"%s %3i %3i %4i %4i %3i %3i %10.5f %10.5f %10.5f\n",
            wname,refpkx-xreg+1,refpky-yreg+1,windhead.naxis1,windhead.naxis2,
            xreg,yreg,cmax,cmin,cavg);

         }
      }
      fclose(fp1);
   }
   else       /*  The user is cross correlating a single search image with the */
   {          /*  reference image                                              */
         wxoff = 0;
         wxend = XPIX;
         wyoff = 0;
         wyend = YPIX;
         if(readfits(wname,&windhead,rdimg,&wxoff,&wxend,&wyoff,&wyend)!=0)
         {
            printf("Program terminating.\n");
            exit(1);
         }
         wxsiz = wxend - wxoff + 1;
         wysiz = wyend - wyoff + 1;
         printf("File %s header read, parameters: %i %i %i\n",fname,wxsiz,wysiz,
           windhead.bitpix);

         if ((windhead.max>100000)||(windhead.max<100))
         {
            printf("Now rescaling search image...\n");
            for (i=0; i<wxsiz; i++)
               for (j=0; j<wysiz; j++)
                  winimg[i][j] = (long) (((double)rdimg[i][j])/
                                        ((double)windhead.max)*100000.0);
         }
         else
         {
            for (i=0; i<wxsiz; i++)
               for (j=0; j<wysiz; j++)
                  winimg[i][j] = (long)(rdimg[i][j]);
         }

      do_xcorr(rxsiz,rysiz,wxsiz,wysiz,lwin,mwin,imin, imax, jmin, jmax,
              ipk,jpk,refpkx,refpky,hfsize,&xreg,&yreg,&cmax,&cmin,&cavg,
              wsmooth,wsmfac);
      xreg += rxoff - wxoff;
      yreg += ryoff - wyoff;
      printf(" %5i %5i %10.5f %10.5f %10.5f\n",xreg,yreg,cmax,cmin,cavg);

      if (outfileopen)
         fprintf(fp2,"%s %3i %3i %4i %4i %3i %3i %10.5f %10.5f %10.5f\n",
           wname,refpkx-xreg+1,refpky-yreg+1,windhead.naxis1,windhead.naxis2,
           xreg,yreg,cmax,cmin,cavg);
   }
   if (outfileopen)
      fclose(fp2);
}     /************************* End main *************************************/

/**************************************************************************/
/*   DO_XCORR - find peak if needed and do the cross correlation based    */
/*              on the selected parameters                                */
/**************************************************************************/

void do_xcorr(xsize, ysize, wxd,wyd,lwin,mwin, imin, imax, jmin, jmax,
        ipk, jpk, refpkx, refpky, hfsize, ireg, jreg, cmax,cmin,cavg,
        wsmooth,wsmfac)
int xsize;              /*  x size of reference image                      */
int ysize;              /*  y size of   "          "                       */
int wxd;                /*  x size of search image                         */
int wyd;                /*  y size of   "      "                           */
int lwin;               /*  x size of search window in search image        */
int mwin;               /*  y size of   "      "     "   "      "          */
int imin;               /*  x min coordinate of optional restricted search */
int imax;               /*  x max      "      "    "         "        "    */
int jmin;               /*  y min      "      "    "         "        "    */
int jmax;               /*  y max      "      "    "         "        "    */
int ipk;                /*  x coord of peak position; if =-999, must find  */
int jpk;                /*  y  "    "   "      "                           */
int refpkx;             /*  x coord of reference peak position             */
int refpky;             /*  y coord of    "        "     "                 */
int hfsize;             /*  size of centroid region, used only if centroid */
int *ireg;              /*  x coord of maximum cross correlation value     */
int *jreg;              /*  y coord of    "      "        "        "       */
float *cmax;            /*  maximum cross correlation value in search region */
float *cmin;            /*  minimum cross correlation value in search region */
float *cavg;            /*  average cross correlation value in search region */
char wsmooth;
int wsmfac;
{

   int  numpix, i, j, k, l,m,n, k_new, l_new, k_off, l_off,
        kmin, kmax, lmin, lmax, num, ko, lo,npix , iwin, jwin;

   double nearsum, xcentrd, ycentrd, dsum;

   double rsqr, lstmax,sum, rmn,ww,ss,ws,rmx;
   char nobadneighbors,
      result;

   int maxpt,            /* list of maximum values in map */
      maxi,                    /* "i" coordinate of maxima */
      maxj;                    /* "j" coordinate of maxima */
   int itmp;

   *ireg=*jreg=-1;
   *cmin=*cmax=*cavg=0.0;

      for (i=0; i<XPIX; i++)
         for (j=0; j<YPIX; j++)
            rdimg[i][j]=0.0;

      if (ipk==-999)
      {                               /* find maxima in image            */
         maxpt = -32765;
         maxi = 0;
         maxj = 0;

printf("%i %i %i %i\n",imin,imax,jmin,jmax);
         for (i = 0; i < wxd; i++)
         {
            for (j = 0; j < wyd; j++)
            {                            /* if pixel is good and in limits,
                                              test for maximum */
               if (((i>=imin) && (i<=imax) &&
                   (j>=jmin) && (j<=jmax)) && (winimg[i][j]<dmax)
                 &&(winimg[i][j]>maxpt))
               {
                  maxpt= winimg[i][j];
                  maxi = i;
                  maxj = j;
               }
            }
         }
         ipk=maxi;
         jpk = maxj;
      }

      if (ipk!=-999)       /*  If peak was found,  */
      {
            /*  Good peak found, now use offsets to add map to others */
         if (do_centroid!=0)
         {
            xcentrd = ycentrd = dsum = 0;
            kmin = max(ipk-hfsize,0);
            kmax = min(ipk+hfsize,XPIX-1);
            lmin = max(jpk-hfsize,0);
            lmax = min(jpk+hfsize,YPIX-1);

            for (k=kmin; k<=kmax; k++)
               for (l=lmin; l<=lmax; l++)
               {
                     xcentrd += k*winimg[k][l];
                     ycentrd += l*winimg[k][l];
                     dsum += winimg[k][l];
                }
            if (dsum!=0)
            {
               ipk = (int)(xcentrd/dsum+0.5);
               if (ipk>=xsize)
                  ipk = -1;
               jpk = (int)(ycentrd/dsum+0.5);
               if (jpk>=ysize)
                  jpk = -1;
            }
            if ((ipk<0) || (jpk<0))
            {
               ipk = -999;
               printf("xcor:  Warning- no valid peak found \n");
            }
         }
      }

      iwin = min(wxd-lwin,max(0,ipk - lwin/2));
      jwin = min(wyd-mwin,max(0,jpk - mwin/2));
      printf(" Peak used is (%i) at  %i , %i\n",winimg[ipk][jpk],ipk+1,jpk+1);

      printf("  Window corner at %i %i\n",iwin+1,jwin+1);

      /*  We have a peak position now, do xcorrelation if selected */

      if (ipk!=-999)
      {
/*  Now do cross correlation                     */
         if (wsmooth)
            winsmooth(wsmfac, xsize, ysize, iwin, jwin, lwin, mwin);

         if (xcorrelate)
         {
            for (j = 0; j< ysize-mwin; j++)
               for (i = 0; i < xsize - lwin; i++)
               {
                  ww = 0.0;
                  ss = 0.0;
                  ws = 0.0;
                  if ( (!speedy) || ((i>refpkx-lwin) && (i<refpkx) &&
                                     (j>refpky-mwin) && (j<refpky)))
                  {
                     for (m=0; m<mwin; m++)
                        for (l=0; l<lwin; l++)
                           if ((winimg[l+iwin][m+jwin]!=0) &&
                               (refimg[i+l][j+m]!=0))
                           {
                              ww += ((double)winimg[l+iwin][m+jwin])
                                    *((double)winimg[l+iwin][m+jwin]);
                              ss += ((double)refimg[i+l][j+m])
                                    *((double)refimg[i+l][j+m]);
                              ws += ((double)winimg[l+iwin][m+jwin])
                                    *((double)refimg[i+l][j+m]);
                           }
                  }
                  else
                  {
                     ww = 1.0;
                     ss = 1.0;
                     ws = 0.01;
                  }
                  if (ww*ss!=0)
                  {
                     rsqr = (ws*ws)/(ww*ss);
                     rdimg[i][j] = sqrt((double)rsqr);
                  }
                  else
                  {
                     rdimg[i][j] = 0;
                  }
               }
            rmn = 2;
            *ireg = -1;
            sum = numpix = rmx = 0;
            for (i = 0; i<xsize-lwin; i++)
               for (j = 0; j<ysize-mwin; j++)
               {
                  if ((rdimg[i][j]<rmn) && (rdimg[i][j]>0.01))
                     rmn = rdimg[i][j];
                  if (rdimg[i][j]>=rmx)
                  {
                     lstmax = rmx;
                     rmx = rdimg[i][j];
                     *ireg = i-iwin;
                     *jreg = j-jwin;
                  }
                  if (rdimg[i][j]>0.01)
                  {
                     sum += rdimg[i][j];
                     numpix++;
                  }
               }
            if  (lstmax == rmx)
               printf("xcor: Warning - no unique correlation max in image\n");
            if (numpix>0)
               *cavg = sum/numpix;
            *cmax =rmx;
            *cmin = rmn;
         }
      }
      else
      {
            /*  No good peak found  in image   */
         printf("xcor: Warning - no good peak found in image\n");
      }
}      /* End do_xcorr */


/*******************************************************************************/
/*   Smooth - boxcar smooth an image with the chosen box width                 */
/*******************************************************************************/
void refsmooth(sfactor,xsize,ysize)
int sfactor;           /*  smooth box size         */
int xsize;             /*  x size of data image    */
int ysize;             /*  y size of data image    */
{
   double newpix,numpix;
   int i,j,k,l,sfd2;

   sfd2 = sfactor/2;
   for (j=0; j<ysize; j++)
      for (i=0; i<xsize; i++)
         winimg[i][j] = refimg[i][j];

   for (j = 0; j<ysize-sfactor; j++)
      for (i = 0; i<xsize-sfactor; i++)
      {
         newpix = 0;
         numpix = 0;
         for (l = j; l<j+sfactor; l++)
            for (k = i; k< i+sfactor; k++)
               if (refimg[k][l]!=0)
               {
                  newpix+=(double)(refimg[k][l]);
                  numpix++;
               }
         if(numpix>0)
            winimg[i][j] = (long)(newpix/numpix);
         else
            winimg[i][j] = 0;
      }
   for (i = 0; i< xsize-sfd2; i++)
      for (j = 0; j<ysize-sfd2; j++)
         refimg[i+sfd2][j+sfd2] = winimg[i][j];
   }  /*  End of smooth routine */

/*******************************************************************************/
/*   Smooth - boxcar smooth an image with the chosen box width                 */
/*******************************************************************************/
void winsmooth(sfactor,xsize,ysize, iwin, jwin, lwin, mwin)
int sfactor;           /*  smooth box size         */
int xsize;             /*  x size of data image    */
int ysize;             /*  y size of data image    */
int iwin;              /*  x window corner */
int jwin;              /*  y window corner */
int lwin;              /*  x window size */
int mwin;              /*  y window size */
{
   double newpix,numpix;
   int i,j,k,l,sfd2,sxst,syst;

   sfd2 = sfactor/2;
   for (j=max(jwin-mwin-1,0); j<min(ysize,jwin+mwin+sfactor); j++)
      for (i=max(iwin-lwin-1,0); i<min(iwin+lwin+sfactor,xsize); i++)
         tmpimg[i][j] = winimg[i][j];

   sxst = max(iwin-sfactor,0);
   syst = max(jwin-sfactor,0);
   for (j = syst; j<min(ysize-sfactor,jwin+mwin+sfd2); j++)
      for (i = sxst; i<min(iwin+lwin+sfd2,xsize-sfactor); i++)
      {
         newpix = 0;
         numpix = 0;
         for (l = j; l<j+sfactor; l++)
            for (k = i; k< i+sfactor; k++)
               if (winimg[k][l]!=0)
               {
                  newpix+=(double)(winimg[k][l]);
                  numpix++;
               }
         if(numpix>0)
            tmpimg[i][j] = (long)(newpix/numpix);
         else
            tmpimg[i][j] = 0;
      }
   for (i=sxst; i<min(iwin+lwin,xsize); i++)
      for (j=syst; j<min(ysize,jwin+mwin); j++)
         winimg[i+sfd2][j+sfd2] = tmpimg[i][j];
   }  /*  End of smooth routine */

/*******************END OF FILE************************************************/