/**************************************************************************
 *                                                                        *
 * Copyright (c) 1996 Michael Richmond and Richard Treffers               *
 *                                                                        *
 *                    This software may be copied and distributed for     *
 *                    educational, research and not for profit services   *
 *                    provided that this copyright and statement are      *
 *                    included in all such copies.                        *
 *                                                                        *
 **************************************************************************/


/*
	make an image in which each pixel is the median of several images

      MEDIAN file1 file2 ... fileN [OUTFILE=] [BIN=] [NOMEAN] [VERBOSE] 
                                   [IQM] [help]

	creates a median frame from a set of images by 
         a. calculating the mean of each image
         b. scaling all to have the same mean
         c. for each pixel, choosing the median value
     the median frame can be used for flat-fielding, for example.

	7/13/92 -added copyheader.c so that the header of
	the first fits file is copied to the outfile -rrt

	9/12/92 - usage changed, commented code deleted -rrt
	9/13/92 - exten removed -rrt
	
	5/20/2001 - added 'iqm' = 'interquartile mean' option.
	            Instead of calculating the median pixel value
		    at each position, determine the mean of all
		    pixel values within second and third quartiles.
		    MWR

*/

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <string.h>
#include <unistd.h>
#include "pcvista.h"
#include "fits.h"

#ifdef PROTO
int main(int, char **);
void check_size(void);
void find_means(void);
void make_median(int iqmflag);
#else
void check_size();
void find_means();
void make_median();
int main();
#endif

#define MAXFILES   (NHANDLE-1)  /* max number of frames to compare -- see */
				/*  limit of # of FITS files open at once */

#define ENDI      (MAXFILES+1)  /* a number >> MAXFILES */
#define MEANFLAG       1 		/* if > 0, calc. mean for each file */
#define MEANBIN        8		/* use every nth pix to calc. mean */
#define IQMFLAG        0                /* if > 0, use interquartile mean */
#define OUTFILE  "median"		/* default output file name */

#define USAGE "usage: median file1 file2 ... fileN [outfile=] [pix=] [nomean] [iqm] [verbose] [help]"

/* some global variables */
static int nfiles;					/* number of files being compared */
static FITS_HANDLE fh[MAXFILES];	/*  their handles */
static char fname[MAXFILES][NBUF];	/*  their names */
static int nr[MAXFILES];			/*  number of rows in each */
static int nc[MAXFILES];			/*  number of cols in each */
static double mfactor[MAXFILES];	/*  scale factor for each rel. to first */
static int nrow, ncol;				/* number of rows, cols in output frame */
static int meanflag;				/* > 0 --> calculate mean for each file */
static int meanbin;					/* use every meanbin'th pixel to calc. mean */
static int iqmflag;                             /* use interquartile mean? */
static char buf[NBUF];
static char outfile[NBUF];			/* output file name */
static int verbose;					/* print out messages if desired */

int 
main(argc, argv)
int argc;
char *argv[];
{
	static int i, endi;
	char *gotit;

	meanflag = MEANFLAG;
	meanbin = MEANBIN;
	iqmflag = IQMFLAG;
	endi = ENDI;
	verbose = 0;
	
	if (argc < 3)
		error(-1, USAGE);

	for (i = 1; i < argc; i++) {
		if ((gotit = find("outfile", argv[i])) != NULL) {
		    if (i < endi)
		    	endi = i;
		    strcpy(outfile, gotit);
		    continue;
		}
		if ((gotit = find("pix", argv[i])) != NULL) {
			if (i < endi)
				endi = i;
			meanbin = (int) (evaluate(gotit) + 0.5);
			continue;
		}
		if (keyword("nomean", argv[i])) {
			if (i < endi)
				endi = i;
			meanflag = 0;
			continue;
		}
		if (keyword("iqm", argv[i])) {
			if (i < endi)
				endi = i;
			iqmflag = 1;
			continue;
		}
		if (keyword("verbose", argv[i])) {
			if (i < endi)
				endi = i;
			verbose = 1;
		}
		if (keyword("help", argv[i])) {
			printf("%s\n", USAGE);
			exit(0);
		}
	}

	if (endi < 3)
		error(-1, USAGE);
	else if (endi == ENDI)
		endi = argc;
	if ((nfiles = endi - 1) > MAXFILES) {
		sprintf(buf, "too many files - MEDIAN can only process %d\n", MAXFILES);
		error(-1, buf);
	}

	/* now we know all the file names start at argument 1 and go
	   up to (endi - 1). so let's process each in turn */
	for (i = 0; i < nfiles; i++) {
		strcpy(fname[i], argv[i + 1]);
		fh[i] = fits_open(fname[i], "r", &nrow, &ncol);
		nr[i] = nrow;
		nc[i] = ncol;
	}
	if (outfile[0] == '\0') {
		strcpy(outfile, OUTFILE);
	}
	
	/* check to make sure that all files are the same size */
	check_size();

	/* next, if NOMEAN has not been specified on the command line,
	   calculate the mean of each file */
	if (meanflag == 1)
		find_means();

	/* and now, create the median file, pixel by pixel */
	make_median(iqmflag);

	for (i = 0; i < nfiles; i++)
		fits_close(fh[i]);

	exit(0);
}

	/* check to make sure that all files are the same size -- if any
	   are not, print the size of all the files and then exit with 
	   an error message */

void
check_size()
{
	int i, j;

	for (i = 0; i < nfiles; i++) {
		if ((nr[i] != nr[0]) || (nc[i] != nc[0])) {
			error(0, "ERROR: at least one of the files is not the same size\n");
			for (j = 0; j < nfiles; j++) {
				sprintf(buf, "file %20s: %5d rows  %5d cols\n", fname[j], 
								nr[j], nc[j]);
				error(0, buf);
			}
			error(-1, "execution terminated\n");
		}
	}
	nrow = nr[0];
	ncol = nc[0];
}


	/* calculate the mean of each file (or use the value given 
	   in the symbol table */

void
find_means()
{
	int i, j, row;
	static int16 data[NMAX];
	static char str[NBUF];
	char *cp;
	double mean, mean0;
	long sum, pixnum;

	/* placate compiler */
	mean0 = 1.0;

	for (i = 0; i < nfiles; i++) {

		/* first, look for value in the symbol table */
		sprintf(str, "mn_%s", fname[i]);
		if ((cp = getsym(str)) != NULL) {
			if (sscanf(cp, "%lf", &mean) == 1) {
				if (verbose) {
					printf("file %15s has mean %9.2f   (in symbol table)\n", 
							fname[i], mean);
					fflush(stdout);
				}
				if (i == 0) {
					mfactor[i] = 1.0;
					mean0 = mean;
				}
				else
					mfactor[i] = mean0/mean;
				continue;
			}
		}
			
		/* otherwise, we have to calculate it */		
		sum = 0;
		pixnum = 0;
		for (row = 0; row < nrow; row++) {
			if (verbose) {
				if (row % 20 == 0) {
					printf("\r  mean of file %12s: up to row %5d out of %5d", fname[i], row, nrow);
					fflush(stdout);
				}
			}
			fits_get_data(fh[i], row, 0, data, ncol);
			for (j = 0; j < ncol; j += meanbin) {
					sum += (long) data[j];
					pixnum++;
			}
		}
		mean = (double) (sum/pixnum);
		if (verbose) {
			printf("\rfile %15s has mean %9.2f   (just calculated)\n", 
					fname[i], mean);
			fflush(stdout);
		}

		if (i == 0) {
			mfactor[i] = 1.0;
			mean0 = mean;
		}
		else
			mfactor[i] = mean0/mean;

		/* place the value of the mean into the PCVISTA symbol table */
		sprintf(buf, "mn_%s=%9.2f", fname[i], mean);   /* update vista file */
		putsym(buf);
		sprintf(buf,"mean=%9.2f", mean);
		putsym(buf);
	}
}


	/*
	 * create a new FITS file of the same size as all the
	 * others and fill each pixel, one by one, with the
	 * median of the values of that pixel from all the
	 * files 
	 *
	 * If the 'iqmflag' is equal to 1, then use the mean
	 * of all pixel values within second and third quartiles,
	 * intead of the true median.
	 *
	 */

void
make_median(int iqmflag)
{
	int i, j, k, l, n_med, row;
	int second_quartile_start, third_quartile_end;
	int num_interquartile_values;
	double sum;
	FITS_HANDLE outfh;
	static int sortnum, sort[MAXFILES];
	static int16 mdata[MAXFILES][NMAX];
	static int16 outdata[NMAX];
	static double mf;

	n_med = ((nfiles + 1)/2) - 1;	  /* the median position of 'nfiles' sorted values */

	/* 
	 * calculate the position at which the second quartile starts,
	 * and the third quartile ends.
	 */
	second_quartile_start = floor(nfiles/4.0);
	third_quartile_end = floor((3.0*nfiles)/4.0) - 1;
	num_interquartile_values = (third_quartile_end - second_quartile_start) + 1;

	/* create the output file */
	outfh = fits_open(outfile, "w", &nrow, &ncol);
	fits_copyheader(fh[0],outfh, FITS_CHECK);	/* copy the first files header to output */
	
	/* now go into a big loop - process the files one row at 
	   a time and write one row at a time into the new file */
	for (row = 0; row < nrow; row++) {

		if (verbose) {
			if (row % 20 == 0) {
					printf("\r   median file %12s: up to row %5d out of %5d",
							 outfile, row, nrow);
					fflush(stdout);
			}
		}
		
		/* read in a row from each file and bring it to the
		   same mean as all the others */
		for (i = 0; i < nfiles; i++) {
			fits_get_data(fh[i], row, 0, mdata[i], ncol);
			if ((i != 0) && (meanflag == 1)) {
				mf = mfactor[i];
				for (j = 0; j < ncol; j++)
					mdata[i][j] *= mf;
			}
		}

		/* now put the median value of each pixel in this row
		   into the output row, and write it out to disk */
		for (j = 0; j < ncol; j++) {

			/* use insertion sort to sort the values of this pixel */
			sort[0] = mdata[0][j];
			for (i = 1; i < nfiles; i++) {
				sortnum = i;
				sort[i] = mdata[i][j];
				for (k = 0; k < sortnum; k++) {
					if (mdata[i][j] < sort[k]) {
						for (l = sortnum; l > k; l--)
							sort[l] = sort[l-1];
						sort[k] = mdata[i][j];
						break;
					}
				}
			}
			/* choose either the true median ... */
			if (iqmflag == 0) {
			 	outdata[j] = sort[n_med];
			}
			/* ... or the interquartile mean */
			else {
				sum = 0.0;
				for (i = second_quartile_start; 
						i <= third_quartile_end; i++) {
					sum += sort[i];
				}
				sum /= num_interquartile_values;
				outdata[j] = sum;
			}


		}
		fits_put_data(outfh, row, 0, outdata, ncol);
	}
	if (verbose) {
		printf("\n");
		fflush(stdout);
	}

	fits_close(outfh);
}