/*********************************************************/
/*   program 6: generate flat field correction image     */
/*********************************************************/

#include <math.h>
#define  XSIZE  640
#define  YSIZE  480

short int *in;
short int *out;

main()
  {
  long  int  x,y;
  double     sum,sum2,fsum,fsum2;
  double     mean,var,stdev;

  /* allocate input and output memory buffers */
  in  = (short int *) malloc (XSIZE*YSIZE);
  out = (short int *) malloc (XSIZE*YSIZE);

  /* read the 16-bit raw data file */
  read_image16("raw",in);

  /* calculate the arithmetic mean
      (and variance and standard deviation) */
  sum = sum2 = 0.0;
  for (y = 0; y < YSIZE; y = y++)
    {
    for (x = 0; x < XSIZE; x = x++)
      {
      sum  = sum +   in[(y*XSIZE)+x];
      sum2 = sum2 + (in[(y*XSIZE)+x] * in[(y*XSIZE)+x]);
      }
    }
  fsum   = sum  / (double)(XSIZE * YSIZE);
  fsum2  = sum2 / (double)(XSIZE * YSIZE);
  mean  = fsum;
  var   = fsum2 - (fsum * fsum);
  stdev = sqrt(var);

  /* subtract the arithmetic mean from each pixel */
  for (y = 0; y < YSIZE; y = y++)
    {
    for (x = 0; x < XSIZE; x = x++)
      {
      out[(y*XSIZE)+x] = in[(y*XSIZE)+x] - (short int)mean;
      }
    }

  /* write the output image */
  /* note: this is a signed image, 16 bits/pixel */
  write_image16("ffcorrect",out);

  /* free memory buffers */
  free (in);
  free (out);
}