/* Program to convolve an unbinned flatfield with an OT shift pattern */
/*
 * Syntax: conflat2 <infile> <otfile> <outfile>
 *
 * options: conflat2 <> <> <> [-b -v N -z 2l -r 2l -t time -yswap]
 *
 * Assumes: unbinned flatfield *with* bias strip of width 32
 *
 * Examples:

     conflat2 sumflat.fits f_ot.103 conflat.103		[normal usage]

   Suppress diagnostic output:

     conflat2 sumflat.fits f_ot.103 conflat.103	-v 0

   Bias subtraction of (short) image flat.067:

     conflat2 flat.067 f_ot.103 conflat.103

 */

#include <stdio.h>
#include <stdlib.h>
#include <math.h>

static char rcsid[] = "$Id: conflat2.c,v 1.3 2003/12/21 14:51:25 jt Exp $";

#define MIN(a,b) (((a) < (b)) ? (a) : (b))
#define MAX(a,b) (((a) > (b)) ? (a) : (b))
#define ABS(a) (((a) > 0) ? (a) : -(a))
#define NINT(x) (x<0?(int)((x)-0.5):(int)((x)+0.5))

#define NAMP 4
#define MAXOT 100000		/* Max number of OT shifts */
#define MAXKERNEL 10000		/* Max number of distinct shift locations */

/* #define TWO15    32768		/* Scale factor for time */
#define TWO15    1		/* Scale factor for time */

#define FULLX 4096		/* Unbinned x size */
#define FULLY 4104		/* Unbinned y size */
#define STDBIAS 32		/* Normal bias count */

struct shiftstruct {
      int dx[2*NAMP];
      int dy[2*NAMP];
} shift[MAXOT];

struct regionstruct {
      int sx;
      int sy;
      int nx;
      int ny;
} regbox[2*NAMP] = {
      0,      0,   2048,    516,
      0,    516,   2048,   1536,
      0,   3588,   2048,    516,
      0,   2052,   2048,   1536,
   2048,      0,   2048,    516,
   2048,    516,   2048,   1536,
   2048,   3588,   2048,    516,
   2048,   2052,   2048,   1536};

struct regionstruct biasbox[NAMP] = {
   4096,      0,     32,   2052,
   4096,   2052,     32,   2052,
   4128,      0,     32,   2052,
   4128,   2052,     32,   2052};

int bias[NAMP];

double time[MAXOT];

double exposure[MAXKERNEL]; 
double texp[MAXKERNEL]; 
int dx[MAXKERNEL], dy[MAXKERNEL]; 
int zero[2*NAMP]={0,0,0,0,0,0,0,0};

double atof();

main(int argc, char **argv)
{
   int i, j, k, nx, ny, nz, iter, minx, miny, maxx, maxy, nkernel;
   int mx, my, sx, sy;
   int ccd, upper, region=7, sig, yswap=0, verbose=2, dobias=0;
   char flatname[256], otname[256], outname[256], line[256], version[256];
   char *head;
   float *data, *sum;
   double inttime=0.0, scale, bb, t;
   FILE *fp;

   if (argc<3) {
      fprintf(stderr, "Syntax: conflat2 <infile> <otfile> <outfile> flags\n");
      fprintf(stderr, "   flags include -v {0|1|2} for verbosity (def: 2)\n");
      fprintf(stderr, "                 -b to do bias subtraction (def: no)\n");
      exit(1);
   }

   strcpy(flatname, argv[1]); argc--; argv++;
   strcpy(otname,   argv[1]); argc--; argv++;
   strcpy(outname,  argv[1]); argc--; argv++;
  
/* Basically I hope that all of these arguments are not wanted.  They
 * are mostly here to patch up bugs in the past
 */

/* Parse the arguments */
   for(i=1; i<argc; i++) {

/* verbose? (-v) */
      if (strcmp(argv[i], "-v") == 0) {
	 verbose = atoi(argv[++i]);
	 continue;
      }

/* bias subtract? (-b) */
      if (strcmp(argv[i], "-b") == 0) {
	 dobias = 1;
	 continue;
      }

/* region? (-r) */
      if (strcmp(argv[i], "-r") == 0) {
	 ccd = atoi(argv[++i]);
	 upper = argv[i][1] == 'u';
	 region = ccd * 2 + (upper?1:0);
	 continue;
      }

/* zero shift region? (-z) */
      if (strcmp(argv[i], "-z") == 0) {
	 ccd = atoi(argv[++i]);
	 upper = argv[i][1] == 'u';
	 zero[2*ccd + (upper?1:0)] = 1;	 
	 continue;
      }

/* time? (-t) */
      if (strcmp(argv[i], "-t") == 0) {
	 inttime = atof(argv[++i]);
	 continue;
      }

   }

   sscanf(rcsid, "%*s %*s %s",version);		/* Get version number */
   if(verbose >= 1) printf("conflat2 v%s\n", version);

/* Read in flatfield image */
   if(rfitsreal(&head, &nx, &ny, &data, flatname) != 0) {
      fprintf(stderr, "Error reading %s\n", flatname);
   }

/* A few sanity checks, obviously can be relaxed */
   if(verbose >= 2) {
      printf("%s read, nx=%d ny=%d \n", flatname, nx, ny);
   }
   if(nx != FULLX+2*STDBIAS || ny != FULLY) {
     fprintf(stderr, "conflat2 only works on unbinned images, size = %d x %d\n",
	     nx, ny);
     exit(1);
   }

/* Allocate memory for the destination image */
   sum = (float *)calloc(nx*ny, sizeof(float));

/* Get the bias levels */
   if(dobias) {
      for(i=0; i<NAMP; i++) {
	 bias[i] = medcount(biasbox[i].sx, biasbox[i].sy, 
			    biasbox[i].nx, biasbox[i].ny, nx, data, &sig);
	 sprintf(line, "CCD%dBIAS", i);
	 rfitshead(head, line, &bb);
	 if(ABS(bias[0]) > 200 && ABS(bias[0]-bb) > 200) {
	    fprintf(stderr, "Bias %d at %d is strange, should be %.0f", 
		    i, bias[i], bb);
	    fprintf(stderr, "Did you compile fitsio.c with the right byte swap?");
	 }
      }
      if(verbose >= 1) {
	 printf("Bias levels at %5d %5d %5d %5d\n", 
		bias[0], bias[1], bias[2], bias[3]);
      }
   } else {
      bias[0] = bias[1] = bias[2] = bias[3] = 0;
   }

/* Get the shift data */
   if( (fp = fopen(otname, "r")) == NULL) {
      fprintf(stderr, "Can't read OT file '%s'\n", otname);
      exit(1);
   }

   fgets(line, 256, fp);
   minx = miny = maxx = maxy = 0;
   nkernel = 0;
   for(nz=0; nz<MAXOT; nz++) {
      if(fgets(line, 256, fp) == NULL) break;
      sscanf(line, "%d %lf %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d", &iter, &time[nz],
	     &shift[nz].dx[0],&shift[nz].dy[0], 
	     &shift[nz].dx[1],&shift[nz].dy[1], 
	     &shift[nz].dx[2],&shift[nz].dy[2], 
	     &shift[nz].dx[3],&shift[nz].dy[3], 
	     &shift[nz].dx[4],&shift[nz].dy[4], 
	     &shift[nz].dx[5],&shift[nz].dy[5], 
	     &shift[nz].dx[6],&shift[nz].dy[6], 
	     &shift[nz].dx[7],&shift[nz].dy[7]);
      for(i=0; i<2*NAMP; i++) {
	 minx = MIN(minx, shift[nz].dx[i]);
	 miny = MIN(miny, shift[nz].dy[i]);
	 maxx = MAX(maxx, shift[nz].dx[i]);
	 maxy = MAX(maxy, shift[nz].dy[i]);
      }
   }
   if(inttime == 0.0) inttime = 2*time[nz-1] - time[nz-2];
   if(verbose >= 1) {
      printf("%s read, nz=%d, x=[%d,%d], y=[%d,%d]\n", 
	     otname, nz, minx,maxx, miny,maxy);
   }
/* For each region perform the shifts */
   for(region=0; region<8; region++) {
      if(!zero[region]) {
	 for(k=0, nkernel=0; k<nz; k++) {
	    for(j=0; j<nkernel; j++) {
	       if(-(shift[k].dx[region] - shift[nz-1].dx[region]) == dx[j] && 
		  -(shift[k].dy[region] - shift[nz-1].dy[region]) == dy[j]) break;
	    }
	    if(j == nkernel) {
	       dx[j] = -(shift[k].dx[region] - shift[nz-1].dx[region]);
	       dy[j] = -(shift[k].dy[region] - shift[nz-1].dy[region]);
	       exposure[j] = 0.0;
	       nkernel++;
	       if(nkernel >= MAXKERNEL) {
		  fprintf(stderr, "nkernel exceeds MAXKERNEL; something's very wrong\n");
		  fprintf(stderr, "proceeding, but it's probably going to be junk\n");
		  nkernel = MAXKERNEL-1;
	       }
	    }
	    exposure[j] += (time[k] - (k==0 ? 0.0 : time[k-1])) / inttime;
	    if(k == nz-1) exposure[j] += (inttime - time[nz-1]) / inttime;
	 }
      } else {
	 nkernel = 1;
	 dx[0] = dy[0] = 0;
	 exposure[0] = 1.0;
      }

      if(verbose >= 2) {
	 printf("region %d: %3d distinct kernel entries, time = %.3f\n", 
		region, nkernel, inttime);
      }
      for(j=0; j<nkernel; j++) texp[j] = exposure[j] * TWO15;

      mx = regbox[region].nx;
      my = regbox[region].ny;
      sx = regbox[region].sx;
      sy = regbox[region].sy;

      for(i=0; i<mx*my; i++) sum[i] = 0.0;

      for(t=0.0, j=0; j<nkernel; j++) {
	 adder(mx,my,nx, &data[sx+sy*nx], 
	       mx,my,mx, sum, 
	       dx[j], (((region/2)&1)&&yswap?-1:1)*dy[j], 
	       texp[j]);
	 if(verbose >= 2) {
	    t += texp[j];
	    printf("\r%7d %8.1f %8.3f %10.1f", j, data[sx+sy*nx+mx/2+(my/2)*nx], 
		   t, sum[mx/2+(my/2)*mx]/TWO15);
	    fflush(stdout);
	 }
      }
      if(verbose >= 2) printf("\n");

      scale = 1.0 / TWO15;
      for(j=0; j<my; j++) {
	 for(i=0; i<mx; i++) {
	    data[sx+sy*nx+i+j*nx] = scale*sum[i+j*mx] - bias[region/2];
	 }
      }
   }
/* Write floating, convolved FITS file */
   chfitshead(&i, head, "BITPIX  ", NULL, -32, 0.0);
   wfitsreal(head, data, outname);
}


/* Adder adds src of size/stride (nx,ny,np) into dest of size/stride
 *  (mx,my,mp), offset by (dx,dy) and multiplied by weight wgt
 */
adder(int nx, int ny, int np, float *src, 
      int mx, int my, int mp, float *dest, int dx, int dy, double wgt)
{
   int j;
   register int count;
   register float *sp;		/* Source pointer */
   register float *dp;		/* Destination pointer */
   if(wgt <= 0) return(0);

   for(j=0; j<ny; j++) {
      if(j+dy < 0 || j+dy > my-1) continue;
      sp = src + j*np;
      dp = dest + (j+dy)*mp;
      if(dx < 0) sp -= dx;
      else       dp += dx;
      count = MIN(nx+dx, mx) - MAX(dx, 0);
      while(count--) *dp++ += *sp++ * wgt;
   }
   return(0);
}


/* Return a median derived from counting the specified subarray */
medcount(int sx, int sy, int nx, int ny, int nsx, float *data, int *sig)
{
   int i, j, k, ntot, median;
   int buf[65536];

   for(i=0; i<65536; i++) buf[i] = 0;
   for(j=sy, ntot=0; j<sy+ny; j++) {
      for(i=sx; i<sx+nx; i++) {
	 if(data[i+nsx*j] >= -0.5 && data[i+nsx*j] <= 65535.9) {
	    k = (int)data[i+nsx*j];
	    buf[k]++;
	    ntot++;
	 }
      }
   }
   for(median=1, i=j=0; median<65536; median++) {
      i += buf[median];
      if(i < ntot/6) j = median;
      if(i > ntot/2) break;
   }
   *sig = MAX(median-j, 1);

   return(median);
}