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


/* 
	more routines used by the STARS command
		
	do the actual computations on the data - find centroid, FWHM,
	and so forth...

	I added a sanity check in "star_axes()" to make sure star isn't
	outside of the tiny box.  Bad bleeding from a saturated star
	in one image caused that to happen, leading to core dump.
	   - MWR 10/4/92

	I changed the method by which stellar centers are calculated from
	first moment to 1-D gaussian fit to marginal light distribution.
	It'll probably slow things down a bit.  
	   - MWR 10/4/92

	I replaced the NR version of 'gauss_fit()' with my own version.
	   - MWR 12/11/92

	Try weighting the marginal sums to give points near the center of
	the box greater weight. I added two parameters, "row_c" and "col_c"
	to "star_axes()" for this purpose. MWR 12/12/92
	(Experiments show equal weights for all points work best, so I've
	 removed the parameters "row_c" and "col_c". MWR 12/13/92)

	use the "do_fwhm()" only if gaussian-fitting fails; otherwise, use the
	   average of gaussian-fit values for FWHM.
	Added checks in "do_fwhm()" to prevent out-of-bounds array references 
	Removed the "normalization" of data before fitting gaussian (but left
	code in place, with #define NORMALIZE)
	  MWR 12/13/92 

	Added sanity check just before calculating centroid: if the sum
	of all pixel values is <= 0, then immediately return with
	"fwhm", "xwidth" and "ywidth" all set to -1.
	  MWR 3/2/96

	Fixed the calculation of second moments in "star_axes()",
	around line 270, by adding sqrt().  Thanks to Shawn Gordon
	and Dick Treffers.
	  MWR 12/16/97

*/

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include "pcvista.h"
#include "stars.h"
 
#ifdef PROTO
static double do_fwhm(int16 **, int sr , int sc, 
	int nr, int nc,
	double cr,double cc,
	int sky); 

static double get_frac(int, int, int);
#else
void star_axes();
static double do_fwhm();
static double get_frac();
#endif


#undef DEBUG


#define NEW      /* means that we use "do_axes" instead of star_axes */


	/* note that star_axes() and do_fwhm() act on all the pixels in a box
	   which starts at (tsr, tsc) and continues for (nr, nc) pixels. */

void
star_axes(image, sr, sc, tsr, tsc, nr, nc, sky, xc, yc, fwhm, xwidth, ywidth)
int16 **image;
int sr, sc, tsr, tsc, nr, nc, sky;
double *xc, *yc, *fwhm, *xwidth, *ywidth;
{
	static int i, row, col, er, ec, bigger, fw_flag;
	static int16 val;
	static float *x, *y, *sig, chisq, amp, center, width;
#ifdef NORMALIZE 
	static float minval;
#endif
	static long lsum, lsumx, lsumy, lsumxx, lsumyy;
	static double pix, sum, sumx, sumy, xc_ax, yc_ax, dxx;
	static double sumxx, sumyy;
#ifdef NEW
	static int fitrad;
	static double peak, eccen, major_axis, minor_axis;
#endif

	/* note that THERE IS SKY in the image supplied to this routine. */

	fw_flag = 0;
	er = tsr + nr;
	ec = tsc + nc;

#ifdef DEBUG
printf("into star_axes: sr=%d sc=%d tsr=%d tsc=%d nr=%d nc=%d er=%d ec=%d\n", sr, sc, tsr, tsc, nr, nc, er, ec);
#endif

#ifdef NEW
	fitrad = (nr > nc ? nc/2 : nr/2);
	if (do_axes(image, tsr, tsc, nr, nc, sky, fitrad, 
			xc, yc, xwidth, ywidth,
			fwhm, &peak, &eccen,
			&major_axis, &minor_axis) != 0) {
		*xc = tsr;
		*yc = tsc;
		*fwhm = -1.0;
		*xwidth = -1.0;
		*ywidth = -1.0;
		return;
	}
	return;
#endif
	



	lsum = lsumx = lsumy = lsumxx = lsumyy = 0;
	sum = 0.0;
	sumx = 0.0;
	sumy = 0.0;
	sumxx = 0.0;
	sumyy = 0.0;

	/* compute the centroid */
	for (row = tsr; row < er; row++) {
		for (col = tsc; col < ec; col++) {
			val = image[row][col] - sky;
			if (val > 0) {
				lsum += val;
				lsumx += val*row;
				lsumy += val*col;
			}
		}
	} 

	/* 
	 * make a sanity check: if the "lsum" is <= 0, then
	 * we have some seriously screwed-up stellar candidate.
	 * In such a case, immediately return
	 * with bogus values (-1.0) in fwhm, xwidth, ywidth.  That should
	 * ensure that this "star" is classified as some kind of noise.
	 */
	if (lsum <= 0) {
#ifdef DEBUG
		printf("lsum <= 0: xc = %lf yc = %lf \n", *xc, *yc);
#endif
		*xc = tsr;
		*yc = tsc;
		*fwhm = -1.0;
		*xwidth = -1.0;
		*ywidth = -1.0;
		return;
	}

	sum = (double) lsum;
	xc_ax = (double) lsumx/ (double) lsum;		/* centroid in x */
	yc_ax =  (double) lsumy/ (double) lsum;		/* centroid in y */


	/* upcoming section, which recalculates the center of light, added 10/4/92
	   by MWR. DO NOT remove the section above - we need it for calculating
	   the second moments later on */

	/* allocate the vectors we'll need */
	if (nr > nc)
		bigger = nr;
	else
		bigger = nc;	
	if ((x = (float *) malloc(bigger*sizeof(float))) == NULL) 
		error(1, "star_axes: can't malloc for X array");
	if ((y = (float *) malloc(bigger*sizeof(float))) == NULL) 
		error(1, "star_axes: can't malloc for Y array");
	if ((sig = (float *) malloc(bigger*sizeof(float))) == NULL) 
		error(1, "star_axes: can't malloc for sig array");
	for (i = 0; i < bigger; i++)
		sig[i] = 1.0;

	/* first do along the row (X) direction */
	for (row = tsr; row < er; row++) {
		x[(row - tsr)] = row;
		y[(row - tsr)] = 0;
		for (col = tsc; col < ec; col++) {
			pix = (double)(image[row][col] - sky);
			y[(row - tsr)] += pix;
		}
	}

#ifdef NORMALIZE
	/* here I "normalize" the marginal sums so that the lowest one is set
	   to zero.  the gaussian fitting routine assumes the baseline is zero */
	minval = y[0];
	for (i = 0; i < nr; i++) {
		if (y[i] < minval)
			minval = y[i];
	}
	for (i = 0; i < nr; i++)
		y[i] -= minval;
#endif

	if (gauss_fit(x, y, sig, nr, &amp, &center, &width, &chisq) == 0) {
		xc_ax = center;
		*xwidth = 2.35*width;	/* convert to FWHM */
	}
	else {
#ifdef DEBUG
		error(0, "star_axes: gauss_fit returns non-zero in X direction");
#endif
		*xwidth = -1;	/* a sign to use 2nd moment below */
	}
	
	/* now do along the col (Y) direction */
	for (col = tsc; col < ec; col++) {
		x[(col - tsc)] = col;
		y[(col - tsc)] = 0;
		for (row = tsr; row < er; row++) {
			pix = (double)(image[row][col] - sky);
			y[(col - tsc)] += pix;
		}
	}

#ifdef NORMALIZE
	/* here I "normalize" the marginal sums so that the lowest one is set
	   to zero.  the gaussian fitting routine assumes the baseline is zero */
	minval = y[1];
	for (i = 0; i < nc; i++) {
		if (y[i] < minval)
			minval = y[i];
	}
	for (i = 0; i < nr; i++)
		y[i] -= minval;
#endif

	if (gauss_fit(x, y, sig, nr, &amp, &center, &width, &chisq) == 0) {
		yc_ax = center;
		*ywidth = 2.35*width;	/* convert to FWHM */
	}	
	else {
#ifdef DEBUG
		error(0, "star_axes: fit_gauss returns non-zero in Y direction");
#endif
		*ywidth = -1;
	}

	free(x);
	free(y);
	free(sig);
	
	
	*xc = xc_ax + sr;
	*yc = yc_ax + sc;
#ifdef DEBUG
printf("xc_ax = %lf xc = %lf  yc_ax=%lf yc = %lf\n", xc_ax, yc_ax, *xc, *yc);
#endif

	/* now make a sanity check: if the centroid falls outside the
	   small box (tsr, er, tsc, ec), then immediately return
	   with bogus values (-1.0) in fwhm, xwidth, ywidth.  That should
	   assure this "star" of being classified as some kind of noise */
	if ((*xc < tsr) || (*xc > er) || (*yc < tsc) || (*yc > ec)) {
#ifdef DEBUG
printf("centroid out of smallbox bounds: xc = %lf yc = %lf \n", *xc, *yc);
#endif
		*xc = tsr;
		*yc = tsc;
		*fwhm = -1.0;
		*xwidth = -1.0;
		*ywidth = -1.0;
		return;
	}


	/* now calculate the second moments to get a measure for the
	   'width' of the image (in case the gaussian-fit failed) */
	for (row = tsr; row < er; row++) {
		dxx = (row - xc_ax)*(row - xc_ax);	
		for (col = tsc; col < ec; col++) {
			val = image[row][col] - sky;
			if (val > 0) {
				sumxx += dxx*val;
				sumyy += (col - yc_ax)*(col - yc_ax)*val;
			}
		}
	} 

	if (*xwidth == -1.0) {
		fw_flag = 1;
		*xwidth = sqrt(sumxx/sum);
	}
	if (*ywidth == -1.0) {
		fw_flag = 1;
		*ywidth = sqrt(sumyy/sum);
	}

	if (fw_flag == 1)
		*fwhm = do_fwhm(image, tsr, tsc, nr, nc, xc_ax, yc_ax, sky); 
	else
		*fwhm = 0.5*(*xwidth + *ywidth);

}

	/* find some sort of measurement of FWHM of object with
	   the given centroid, and return the value.  Remember that
	   the SKY is still there in this routine. */

	/* added checks to keep array references inside box boundaries
	     ... MWR 12/13/92 */

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

#ifdef DEBUG
printf("into do_fwhm, sr=%d sc=%d cr=%lf cc=%lf  here comes the data\n", sr, sc, cr, cc);
#endif

	cent_row = (int) (cr + 0.5);		/* approximations to center */
	cent_col = (int) (cc + 0.5);		/* in box coords */
	/* check to make sure center pixel in bounds */
	if (cent_row >= sr + nr)
		cent_row = (sr + nr) - 1;
	if (cent_col >= sc + nc)
		cent_col = (sc + nc) - 1;
	halfmax = (image[cent_row][cent_col] - sky)/2 + sky;

	for (i = cent_col; i > sc; i--) 
		if (image[cent_row][i] < halfmax)
			break;
	if (i + 1 >= sc + nc) {
		d1 = 0.0;	
	}
	else {
		frac = get_frac(image[cent_row][i+1], image[cent_row][i], halfmax);
		d1 = (cent_col - (i+1)) + frac;
	}

	for (i = cent_col; i < (sc + nc) - 1; i++) 
		if (image[cent_row][i] < halfmax)
			break;
	if (i < sc) {
		d2 = 0.0;
	}
	else {
		frac = get_frac(image[cent_row][i-1], image[cent_row][i], halfmax);
		d2 = ((i - 1) - cent_col) + frac;
	}
	
	for (j = cent_row; j > sr; j--)
		if (image[j][cent_col] < halfmax)
			break;
	if (j + 1 >= sr + nr) {
		d3 = 0.0;
	}
	else {
		frac = get_frac(image[j+1][cent_col], image[j][cent_col], halfmax);
		d3 = (cent_row - (j+1)) + frac;
	}
	
	for (j = cent_row; j < (sr + nr) - 1; j++)
		if (image[j][cent_col] < halfmax)
			break;
	if (j - 1 < sr) {
		d4 = 0.0;
	}
	else {
		frac = get_frac(image[j-1][cent_col], image[j][cent_col], halfmax);
		d4 = ((j - 1) - cent_row) + frac;
	}

	return(((d1 + d2) + (d3 + d4)) / 2.0);
}

	/* given two numbers 'larger' and 'smaller', return the fraction of
	   the way from 'larger' to 'smaller' that you have to go to get
	   to the number 'middle'. 

	   i.e.  larger=10, smaller=4, middle=6.0 --> return 0.66 

	  ******/

static double 
get_frac(larger, smaller, middle)
int larger, smaller, middle;
{
	double a, b;

	a = (double) (larger - smaller);
	b = (double) (larger - middle);
	if ((a <= b) || (a <= 0.0))
		return(0.0);
	return(b/a);
}