/* @(#)fold_image.c	17.1.1.1 (ESO-IPG) 01/25/02 17:51:59 */
/* 
   
   Otmar Stahl

   fold_image.c

   simple folding routine 
   strong alias peaks with image slicer!

*/

/* system includes */

#include <math.h>

/* FEROS specific includes */

#include <proto_nrutil.h>
#include <misc.h>

int fold_image 
#ifdef __STDC__
(
 float imageo[], float imagef[], float template[],
 int offset, int npix[], int npixi
 )
#else
     (
      imageo, imagef, template, offset, npix, npixi
      )
     float imageo[], imagef[], template[];
     int offset, npix[], npixi;
     
#endif

{
  int width, first, last, center;
  float sum, xcorr, mini;
  float *templ;
  
  int i, jj;

  templ = vector(0, npixi);

  width = (npixi - 1) / 2;
  first = width;
  last = npix[0] - width;

  /* normalize template */

  mini = 3.0e34;  
  for (jj = -width; jj <= width; jj++)
    {
      if (mini > template[jj + width])
	{
	  mini = template[jj + width];
	}
    }

  /* subtract minumum  */

  xcorr = 0.0;
  for (jj = -width; jj <= width; jj++)
    {
      xcorr += SQUARE(template[jj + width] - mini);
    }
  xcorr = (float) sqrt((double) xcorr);

  /* and normalize auto-correlation of template */

  for (jj = -width; jj <= width; jj++)
    {
      templ[jj + width] = (template[jj + width] - mini) / xcorr;
    }

  /*  printf("template normalized\n"); */

  for (i = first ; i <= last; i++)
    {
      center = offset + i;

      /* subtract minimum */

      mini = 3.0e34;
      for (jj = -width; jj <= width; jj++)
	{
	  if ( mini > imageo[center + jj])
	    {
	      mini = imageo[center + jj];
	    }
	}

      /* compute auto-correlation */

      xcorr = 0.0;
      for (jj = -width; jj <= width; jj++)
        {
          xcorr += SQUARE(imageo[center + jj] - mini);
        }
      xcorr = (float) sqrt((double) xcorr);

      /* and compute cross-correlation */

      sum = 0.0;
      for (jj = -width; jj <= width; jj++)
	{
	  sum += (imageo[center + jj] - mini) / xcorr * (templ[jj + width]);
	}
      imagef[i] = sum;
    }
  return 0;
}