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



/*  computes statistics on an object in the specified BOX of an image
    usage:

			AXES source  [box=d] [sky=]

	prints out values on screen, and sets AXC and AXR symbol values.
	also calculates a (somewhat simple) FWHM estimate. 

	1/30/91 - added 'floor' to "sky" expression   - MWR

	11/18/91 - added a 'sec mom' to the last line of printout, which is
	           an approximation to the second moment of the light distribution
	           around the centroid.  It's the average of x-2nd-moment and 
	           y-2nd-moment.  - MWR

	11/18/91 - removed the restriction that a BOX must be specified - if not,
	           the whole image is used (surprise).  - MWR
	2/21/92  - added exit(0) -rrt
	2/25/92  - a. fixed bug in gotit=find(sky,argv)
				b. put 'continues in the argc loop
				c. added help option
				d. converted sky to a double and added skyflag
				e. added SMALL_SUM test before computing centroid
				f. a test that centroid lies within the box
				g. added matherr function
	9/23/92  - usage -rrt

	10/4/92  - changed centering function to 1-d gaussian fit along each
	             direction, instead of first moment along each direction.
	             also added "gaussian FWHM" estimates in output.  - MWR
	10/29/92 - added more printout to matherr -rrt

	12/11/92 - replaced NR gaussian-fitting routine with my own. MWR
	4/22/94  - replace nr.h with nr.h

	1/13/96  - subtract sky from the data before calling the routines
	              that calculate object properties.  Will simplify
	              code quite a bit, at the cost of rounding off
	              the sky value to an integer.  I think it's an
	              overall win.  MWR

	8/2/96   - add a "reflection" for every point, at an equal 
	              distance from the origin, but in negative direction. 
	              MWR
	
	8/13/96  - now we _don't_ subtract sky before calling other 
	              routines, but instead pass sky value to them.
	              avoids round-off errors.  MWR
*/

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

#ifdef PROTO
int main(int,char **);
static int16 guess_sky(int16 **p, int nr, int nc);
static void do_it(int16 **, int, int, int, int, double);
static void do_fwhm(int16 **, int, int, double, double, double);
#else
int main();
static int16 guess_sky();
static void do_it();
static void do_fwhm();
#endif

/*	this is the smallest total number of pixels above sky to accept
	I really don't know what a reasonable value should be 
	previously no testing at all was done
*/

#define SMALL_SUM 10.0	

#define USAGE "usage: axes file [box=] [sky=] [help]"


int
main(argc,argv)
int argc;
char *argv[];
{
	FITS_HANDLE fh;
	static int i, boxflag, skyflag;
	static int sc, sr, nr, nc, boxno, nrow, ncol;
	static int16 **p;
	char *cp, *gotit, fname[NBUF];
	double sky;

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

	for (i = 1; i < argc; i++) {
		if (keyword( "help", argv[i]) ){
			printf("%s\n", USAGE);
			exit(0);
		}
		if ((gotit = find("box", argv[i])) != NULL) {
			boxno = (int)(evaluate(gotit) + 0.5);
	 		boxflag = 1;
			continue;
		}
		if ((gotit = find("sky", argv[i])) != NULL) {
			sky = evaluate(gotit) ;
			skyflag=1;
			continue;
		} 
		if(i ==1){
			strcpy(fname, argv[1]);
		}else{
			fprintf(stderr,"unknown argument '%s'\n",argv[i]);
			exit(1);
		}
	}

	fh = fits_open(fname, "r", &nrow, &ncol);
	nr = nrow;
	nc = ncol;

	if (!boxflag) {
		boxno = -1;
	}else{
		if (getbox(boxno, &sr, &sc, &nr, &nc))
			error(-1,"box not found");
	 	check_box(boxno, (int) nrow, (int) ncol);
	}

	/* use the supplied value of SKY if one is given, or symbol-table */
	if (skyflag == 0) {	
		if ((cp = getsym("sky")) != NULL) {
			if (sscanf(cp, "%lf", &sky) == 1)
				skyflag = 1;
		}
	}

	/* read the data from the FITS image file */
	if ((p = boxdata(fh, boxno, (int) nrow, (int) ncol, &nr, &nc)) == NULL)
		error(-1,"couldn't read data in box from file");

	/* if still no sky value, guess one from the edges of the data */
	if (skyflag == 0) {
		sky = guess_sky(p, nr, nc);
		skyflag = 1;
	}

	printf("using a sky value of %.1f\n", sky);

	do_it(p, sr, sc, nr, nc, sky);
	exit(0);
}


	/*
	 * find a guess at the "sky" value by picking the largest
	 * pixel value around the edges of the box of data.
	 */

static int16
guess_sky(p, nr, nc)
int16 **p;
int nr, nc;
{
	static int i;
	static int16 sky;

	sky = p[0][0];
	for (i = 0; i < nc; i++) {
		if (p[0][i] > sky) {
			sky = p[0][i];
		}
	}
	for (i = 0; i < nc; i++) {
		if (p[nr-1][i] > sky) {
			sky = p[nr-1][i];
		}
	}
	for (i = 0; i < nr; i++) {
		if (p[i][0] > sky) {
			sky = p[i][0];
		}
	}
	for (i = 0; i < nr; i++) {
		if (p[nr-1][0] > sky) {
			sky = p[nr-1][0];
		}
	}

	return(sky);
}


	/* 
	 * Call the routines that do the actual computations on the data ---
	 * finding centroid, FWHM, and so forth... 
	 */

static void 
do_it(p, sr, sc, nr, nc, sky)
int16 **p;
int sr, sc, nr, nc;
double sky;
{
	double fitrad;
	double row_cen, col_cen, fwhm, peak, eccen, major_axis, minor_axis;
	double row_width, col_width;
	char buf[20];

	/* make the fit out as far as the smaller of the box dimensions */
	fitrad = (nr > nc ? nc : nr);

	/* compute stuff */
	if (do_axes(p, 0, 0, nr, nc, sky, fitrad, &row_cen, &col_cen, 
                    &row_width, &col_width, &fwhm,
	            &peak, &eccen, &major_axis, &minor_axis) != 0) {
		error(0, "axes: do_axes returns with error");
		return;
	}
	row_cen += sr;
	col_cen += sc;

	/* print out centroid of object to user, and put in symbol table */
	printf("centroid    row=%7.2f  col=%7.2f\n", row_cen, col_cen);
	sprintf(buf, "axr=%-7.2f", row_cen);
	putsym(buf);
	sprintf(buf, "axc=%-7.2f", col_cen);
	putsym(buf);

	/* print out the circular FWHM we found for user */
	printf("circular gaussian --> FWHM %.3f  peak %.1f\n", fwhm, peak);

	/* now use a simpler method to find FWHM, as a check */
	do_fwhm(p, nr, nc, row_cen - sr, col_cen - sc, sky);

	/* print out eccentricity, major anx minor axes */
	printf("eccen %5.3f maj axis %6.2f deg, min axis %6.2f deg\n",
		eccen, major_axis, minor_axis);
}

	/* 
	 * find a quick approximation to the FWHM of object 
	 * by walking down from the central pixel (whose coords are
	 * given as arguments) in each of the four cardinal directions,
	 * along rows and columns.  Figure out how many pixels one
	 * must go to reach a pixel value half that of the central
	 * value.  Use linear interpolation between the two pixels
	 * that bracket the half-maximum point.
	 *
	 * We subtract the supplied "sky" value from all pixels before
	 * we start the job -- this is just an approximate method, so
	 * rounding doesn't matter.
	 */

static void
do_fwhm(p, nr, nc, cr, cc, sky)
int16 **p, nr, nc;
double cr, cc, sky;
{
	int i, j, cent_col, cent_row;
	double halfmax, frac, d1, d2, d3, d4, fwhm;

	/* subtract sky from each pixel -- will make life simpler below */
	for (i = 0; i < nr; i++) {
		for (j = 0; j < nr; j++) {
			d1 = p[i][j] - sky;
			if (d1 < MIN_DATA_VAL) {
				p[i][j] = MIN_DATA_VAL;
			}
			else if (d1 > MAX_DATA_VAL) {
				p[i][j] = MAX_DATA_VAL;
			}
			else {
				p[i][j] = d1;
			}
		}
	}
	
	cent_row = (int)(cr + 0.5);		/* approximations to center */
	cent_col = (int)(cc + 0.5);
	halfmax = p[cent_row][cent_col]/2.0;

	/* 
	 * first, we walk backwards along the central row.
	 * Note that if we don't reach the half-max point by the
	 * time we hit the edge of the data, we extrapolate
	 * linearly beyond the data's edge. 
	 */
	for (i = cent_col; i > 0; i--) {
		if (p[cent_row][i] < halfmax) {
			break;
		}
	}

	/*
	 * to do the linear interpolation, let's imagine that we are
	 * looking on an x-y graph between points (x1,y1) and (x2,y2),
	 * where x2 > x1 are adjacent integers.  What 
	 * we know is the value of some value y, with y1 < y < y2.
	 * We want to figure out the value of x which corresponds to it.
	 * In general, we can write
	 *
	 *                   x = x2 - (x2 - x1) * [ (y2 - y)/(y2 - y1) ]
	 *
	 * but since x2 - x1 = 1, we get
	 *
	 *                   x = x2 - (y2 - y)/(y2 - y1)
	 *
	 * Thus, the fraction of a pixel we have to move from pixel
	 * p[cent_row][i + 1] outwards from the center is 
	 *
	 *                      (y2 - y)/(y2 - y1)
	 *
	 * and the total distance from the central pixel is therefore
	 * 
	 *       (cent_col - (i + 1)) + (y2 - y)/(y2 - y1)
	 *
	 * We have to make a check to see if the two adjacent pixels
	 * have exactly equal values, so that we don't divide by zero.
	 * in such a case, just set "frac" = 0.5.
	 */
	if (p[cent_row][i + 1] == p[cent_row][i]) {
		frac = 0.5;
	} 
	else {
		frac = ((double) (p[cent_row][i + 1] - halfmax)) / 
				 (p[cent_row][i + 1] - p[cent_row][i]);
	}
	d1 = (cent_col - (i + 1)) + frac;

	/* now we walk from center forwards along the row. */
	for (i = cent_col; i < nc - 1; i++)  {
		if (p[cent_row][i] < halfmax)
			break;
	}
	if (p[cent_row][i - 1] == p[cent_row][i]) {
		frac = 0.5;
	}
	else {
		frac = ((double) (p[cent_row][i - 1] - halfmax)) / 
				 (p[cent_row][i - 1] - p[cent_row][i]);
	}
	d2 = ((i - 1) - cent_col) + frac;
	
	/* now we walk from center backwards along a column */
	for (j = cent_row; j > 0; j--) {
		if (p[j][cent_col] < halfmax)
			break;
	}
	if (p[j + 1][cent_col] == p[j][cent_col]) {
		frac = 0.5;
	}
	else {
		frac = ((double) (p[j + 1][cent_col] - halfmax)) / 
				 (p[j + 1][cent_col] - p[j][cent_col]);
	}
	d3 = (cent_row - (j + 1)) + frac;
	
	/* finally, walk forward along the central column */
	for (j = cent_row; j < nr - 1; j++) {
		if (p[j][cent_col] < halfmax)
			break;
	}
	if (p[j - 1][cent_col] == p[j][cent_col]) {
		frac = 0.5;
	}
	else {
		frac = ((double) (p[j - 1][cent_col] - halfmax)) / 
				 (p[j - 1][cent_col] - p[j][cent_col]);
	}
	d4 = ((j - 1) - cent_row) + frac;

	fwhm = ((d1 + d2) + (d3 + d4)) / 2.0;

	printf("get HM of %.3f %.3f %.3f %.3f ---> FWHM %.3f\n", 
			d1, d2, d3, d4, fwhm);

}

#ifndef LINUX

#if 0
int matherr(x)
/*register*/ struct exception *x;
{
	char errbuf[80];

	sprintf(errbuf, "math error in %s function argument=%3f",
			x->name, x->arg1);
	error(1,errbuf);
	return 0;
}
#endif 

#endif  /* #ifndef LINUX */