/**************************************************************************
 *                                                                        *
 * 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.                        *
 *                                                                        *
 **************************************************************************/


/*************************************************************************
 * 
 * apply a median-ring filter to an image, creating an "unsharp mask"
 * of the original.  One will usually subtract the result from the 
 * original, yielding a "sharpened" version.  
 * 
 * The process is as follows:
 *    - form a square, NxN kernel specified by the user
 *    - for each pixel[i][j] in the image, create a set of NxN
 *          numbers by placing the kernel down at position [i][j]
 *          in the image and multiplying each element of the kernel
 *          by the pixel beneath it
 *    - place the median of those NxN values into the output[i][j]
 *
 * Thus, this program tends to _destroy_ small features, leaving
 * only large ones.
 *
 *    usage: RING file1 filter= [outfile=] [verbose]
 *
 * The "filter" argument must be the name of a file with an ASCII 
 * representation of the filter, something like this:
 *
 *      0 1 0 0 0 
 *      0 0 0 0 1    In this case, the filter is a 5x5 array, in which
 *      0 0 0 0 0    all non-zero elements have the same weight.  One
 *      1 0 0 0 0    could use fractional weights, if desired.
 *      0 0 0 1 0
 *
 *
 * 5/21/95 MWR
 **************************************************************************/

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

#undef DEBUG

#ifdef PROTO
int main(int, char **);
static void create_ring(char *file);
static void fill_element(int elem, int row, int col, double weight);
static void make_median(void);
#else
static void create_ring();
static void fill_element();
static void make_median();
#endif

	/* 
	 * we use this structure to hold offsets of each non-zero element
	 * of the ring filter in an array.
	 */
struct s_ring {
	int row;                /* offset from center in rows */
	int col;                /* offset from center in cols */
	double  weight;         /* weight assigned to this pixel */
};


#define USAGE "usage: ring file filter= [outfile=] [verbose]"

#define ENDI     100        /* max possible number of arguments */
#define OUTFILE  "ring"		/* default output file name */
#define MAXPTS   300        /* max # of non-zero elements in ring filter */


/* some global variables */
static FITS_HANDLE fh;				/*  input file handle */
static char fname[NBUF];			/*  input file name */
static int nr;						/*  number of rows in input file */
static int nc;						/*  number of cols in input file */
static int nrow, ncol;				/* number of rows, cols in output frame */
static char filterfile[NBUF];		/* name of ASCII filter file */
static char outfile[NBUF];			/* output file name */
static int verbose;					/* print out messages if desired */
static int nring;					/* num of non-zero elems in ring filter */
static struct s_ring ring_array[MAXPTS];   /* holds offset of non-zero elems */

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

	verbose = 0;
	
	if (argc < 2) {
		error(-1, USAGE);
        }
	if (strcmp(argv[1], "help") == 0) {
		printf("%s\n", USAGE);
		exit(0);
	}
        

	/* first argument must be input file name */
	strcpy(fname, argv[1]);

	for (i = 2; i < argc; i++) {
                if (keyword("help", argv[i])) {
                    printf("%s\n", USAGE);
                    exit(0);
                }
		if ((gotit = find("filter", argv[i])) != NULL) {
		    strcpy(filterfile, gotit);
		    continue;
		}
		if ((gotit = find("outfile", argv[i])) != NULL) {
		    strcpy(outfile, gotit);
		    continue;
		}
		if (keyword("verbose", argv[i])) {
			verbose = 1;
			continue;
		}
		/* if we got here, improper option was given */
		fprintf(stderr, "ring: unknown argument %s\n", argv[i]);
		exit(1);
	}

	if (filterfile[0] == '\0') {
		error(0, "no filter file name supplied");
		error(-1, USAGE);
	}

	fh = fits_open(fname, "r", &nrow, &ncol);
	nr = nrow;
	nc = ncol;
	if (outfile[0] == '\0') {
		strcpy(outfile, OUTFILE);
	}
	
	/* create the ring array filter */
	create_ring(filterfile);

#ifdef DEBUG
	for (i = 0; i < nring; i++) {
		printf("elem %2d  row %3d col %3d = %6.1f\n", 
				i, ring_array[i].row, ring_array[i].col,
				ring_array[i].weight);
	}
#endif


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

	fits_close(fh);

	exit(0);
}


	/*
	 * This function creates the ring filter; it sets "nring" to the
	 * number of non-zero elements in the filter, and then fills the
	 * first "nring" elements of ring_array with the offsets and weights
	 * of those non-zero elements.
	 *
	 * The given file must have an ASCII picture of the desired filter,
	 * which must be in the form of a square, with every element having
	 * a value.  Values can be integers or fractional, but the usual
	 * case is all ones and zeros, like so:
	 *
	 *     0 1 0 1 0
	 *     1 0 0 0 1
	 *     0 0 0 0 0
	 *     1 0 0 0 1
	 *     0 1 0 1 0
	 *
	 * The file must have no blank lines, and all lines must have the
	 * same number of elements.
	 *     
	 */

#define LINELEN 200

static void
create_ring(filterfile)
char *filterfile;
{
	int length, index, pos, row, col;
	double val;
	char line[LINELEN + 1];
	FILE *fp;

	if ((fp = fopen(filterfile, "r")) == NULL) {
		error(1, "can't open filter file");
	}
	index = 0;
	row = 0;
	while (fgets(line, LINELEN, fp) != NULL) {
		length = strlen(line);
		pos = 0;
		while ((line[pos] == ' ') && (pos < length)) {
			pos++;
		}
		col = 0;
		while (pos < length) {
			if (sscanf(line + pos, "%lf", &val) != 1) {
				error(1, "bad format in filter file");
			}
			if (val != 0) {
				fill_element(index, row, col, val);
				if (++index >= MAXPTS) {
					error(1, "too many non-zero elements in filter");
				}
			}
			col++;
			/* now skip to the end of this number */
			while ((line[pos] != ' ') && (line[pos] != '\0')) {
				pos++;
			}
			/* and skip to start of next number */
			while (line[pos] == ' ') {
				pos++;
			}
		}
		row++;
	}

	/* 
	 * now, after having read all the elements, we finally know how many
	 * elements there are.  So, we can go back and modify the offsets
	 * to make them be relative to the center.
	 *
	 * We've read "row" lines from the file.  So, let's assume that there 
	 * are "row*row" elements.
	 */
	nring = index;
	for (index = 0; index < nring; index++) {
		ring_array[index].row -= row/2;
		ring_array[index].col -= row/2;
	}

	fclose(fp);
}

   /*
    * Fill the given element of "ring_array" with the given information.
    */

static void
fill_element(elem, row, col, weight)
int elem, row, col;
double weight;
{
   if (elem >= MAXPTS) {
		error(1, "fill_elements: too many ring elements");
   }
   ring_array[elem].row = row;
   ring_array[elem].col = col;
   ring_array[elem].weight = weight;
}


	/* create a new FITS file of the same size as all the
	   input and fill each pixel, one by one, with the
	   median of the values from the ring filter */

static void
make_median()
{
	int i, j, k, l, row, col, r, c;
	int npix, nbytes;
	FITS_HANDLE outfh;
	static int sortnum, sort[MAXPTS];
	static int16 mdata[MAXPTS];
	static int16 **in_data, **out_data;


	/* create the output file */
	outfh = fits_open(outfile, "w", &nrow, &ncol);
	fits_copyheader(fh, outfh, FITS_CHECK);	
	
	/* 
	 * we're going to try to read two copies of the input file into memory,
	 * so attempt to allocate the memory here.
	 */
	if ((in_data = (int16 **) malloc(nrow*sizeof(int16 *))) == NULL) {
		error(1, "can't allocate for input array of pointers ");
	}
	if ((out_data = (int16 **) malloc(nrow*sizeof(int16 *))) == NULL) {
		error(1, "can't allocate for output array of pointers ");
	}
	nbytes = sizeof(int16)*ncol;
	for (row = 0; row < nrow; row++) {
		if ((in_data[row] = (int16 *) malloc(nbytes)) == NULL) {
			error(1, "can't allocate for input data array");
		}
		if ((out_data[row] = (int16 *) malloc(nbytes)) == NULL) {
			error(1, "can't allocate for output data array");
		}
	}


	/* read the entire file into memory */
	for (row = 0; row < nrow; row++) {
		fits_get_data(fh, row, 0, in_data[row], ncol);
	}

	/*
	 * and now, with the entire image in memory, perform the filtering
	 * pixel-by-pixel.
	 */
#ifdef DEBUG
	printf("output about to appear\n");
#endif

	for (row = 0; row < nrow; row++) {

		if (verbose) {
			if (row % 20 == 0) {
					printf("\r   up to row %5d out of %5d", row, nrow);
					fflush(stdout);
			}
		}

		for (col = 0; col < ncol; col++) {

#ifdef DEBUG
			printf("pixel [%d][%d] has vals ", row, col);
#endif

			/* fill "mdata" with all valid values from the ring */
			npix = 0;
			for (i = 0; i < nring; i++) {
				r = row + ring_array[i].row;
				c = col + ring_array[i].col;
				if ((r >= 0) && (r < nrow) && (c >= 0) && (c < ncol)) {
					mdata[npix] = in_data[r][c]*ring_array[i].weight;
#ifdef DEBUG
					printf(" %6d", mdata[npix]);
#endif
					npix++;
				}
			}

			/* now find the median of values in array "mdata" */

			/* use insertion sort to sort the values of this pixel */
			sort[0] = mdata[0];
			for (j = 1; j < npix; j++) {
				sortnum = j;
				sort[j] = mdata[j];
				for (k = 0; k < sortnum; k++) {
					if (mdata[j] < sort[k]) {
						for (l = sortnum; l > k; l--)
							sort[l] = sort[l-1];
						sort[k] = mdata[j];
						break;
					}
				}
			}

			out_data[row][col] = sort[npix/2];
#ifdef DEBUG
			printf(" --> %6d \n", out_data[row][col]);
#endif
		}
#ifdef DEBUG
			printf("\n");
#endif
	}

	/*
	 * the filtering is done, so we can write out the output image
	 */
	for (row = 0; row < nrow; row++) {
		fits_put_data(outfh, row, 0, out_data[row], ncol);
	}

	if (verbose) {
		printf("\n");
		fflush(stdout);
	}

	fits_close(outfh);
}