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


/****************************************************************************
 * A routine to find the sky value in an image, or subsection thereof.	
 *
 * This used to be part of the "sky.c" code, but I removed it and made
 * it a library routine so that other programs can call it as well.
 *
 * MWR 1/19/2003
 */


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

#ifndef DEBUG
#undef DEBUG					/* define this to get some diag output */
#endif

int find_sky(FITS_HANDLE fh, int sr, int sc, int nr, int nc, int nskip,
		    int binsize, int16 *sky, float *skysig);
static void make_hist(FITS_HANDLE fh, int sr, int sc, int nr, int nc, 
                      int nskip, int binsize, int min, int max);
static void find_range(int peak, double fraction, int *botbin, int *topbin);
static int find_fitted_values(int botbin, int topbin,
                              float *peak, float *sigma);
static float * vector(int nel);


#define FRACTION  0.40     /* use only bins w/ at least this fraction */
#define MINBIN    5	   /* use at least this many bins around the peak */
                   	   /*   on either side of the peak, when fitting */
                   	   /*   a gaussian to find sky sigma             */

static long hist[NHIST + 1];


	/********************************************************************
	 * given an image and the coords of a box within that image,
	 * find the a sky value for the image by looking at
	 * the peak of the histogram of data values.  Use only those pixels 
	 * in the general vicinity of the peak to find the variance of values
	 * around the peak.
	 *
	 * return 0 if all is well, or -1 if there's a problem.  
	 * If all went well, place the sky value into the 'sky' parameter and
	 * the skysigma into the 'skysig' parameter.  Otherwise, the values
	 * of the two arguments are undefined (but probably 0.0).
	 */

int
find_sky(fh, sr, sc, nr, nc, nskip, binsize, sky, skysig)
FITS_HANDLE fh;
int sr, sc, nr, nc, nskip, binsize;
int16 *sky;
float *skysig;
{
	int i, peak, botbin, topbin, min, max, ret;
	long binmax;
	float skypeak;

	ret = 0;

	/* now construct the first histogram, using the min and max values 
	   as the first and last bins */
	min = MIN_DATA_VAL;
	max = MAX_DATA_VAL;

	make_hist(fh, sr, sc, nr, nc, nskip, binsize, min, max);

#ifdef DEBUG
	printf("  initial hist %-5d bins, min=%-7d  max=%-7d\n",
				NHIST, min, max);
#endif

	/* find the peak of the histogram, and set the sky equal to the value
	   at the peak. However, if the peak bin is the first, or last, issue
	   a warning error message, because the job was probably not done
	   correctly. */
	peak = 0;
	binmax = 0;
	for (i = 0; i < NHIST; i++) {
		if (hist[i] > binmax) {
			peak = i;
			binmax = hist[i];
		}
	}

        /*
         * before 2/16/96, we used to set the sky value equal to this
         * bin with the maximum number of pixels, like so:
                *sky = (int16) (min + peak);
         * but now we use the center of a fitted gaussian ... see below
         */

	if ((peak == 0) || (peak == (NHIST - 1))) {
		ret = -1;
	}	


	/* find the bins which contain at least FRACTION of the number of pixels
	   as the peak bin, centered on the peak of the histogram */
	find_range(peak, FRACTION, &botbin, &topbin);
#ifdef DEBUG
	printf("  botbin %-5d  peak %-5d  topbin %-5d\n", botbin, peak, topbin);
#endif

        /*
         * now we fit a gaussian to the pixel values in this range,
         * and use the fitted peak and sigma of the gaussian
         * for output we'll return to the user.
	 * 
	 * Don't forget to convert the fitted "skypeak" and "skysig"
	 * from _binned_ units back into real units!
         */
        if (find_fitted_values(botbin, topbin, &skypeak, skysig) != 0) {
                ret = -1;
        }
        *sky = (int16) (min + skypeak*binsize);
	*skysig *= binsize;

	return(ret);
}
	
	
	/* make a histogram for the given image, ignoring any points
	   smaller than 'min' or larger than 'max' */

static void
make_hist(fh, sr, sc, nr, nc, nskip, binsize, min, max)
FITS_HANDLE fh;
int sr, sc, nr, nc, nskip, binsize, min, max;
{
	int i, j, bin;
	int16 val, data[NMAX + 1];

	for (i = 0; i < NHIST; i++)
		hist[i] = 0;

	for (i = 0; i < nr; i++) {
		fits_get_data(fh, sr + i, sc, data, nc);
		for (j = 0; j < nc; j += nskip) {
			val = data[j];
			if ((val >= min) && (val <= max)) {
				bin = (int) ((val - min)/binsize);
				hist[bin]++;
			}
		}
	}
}


	/* by going through the histogram array (assumed to have been
	   calculated already!), find the range of bin values, centered 
	   on the 'peak'th bin, between which the bins have at least FRACTION
	   times the number of counts of the peak bin. return the min and max 
	   bin in the 'botbin' and 'topbin' parameters. */


static void
find_range(peak, fraction, botbin, topbin)
int peak, *botbin, *topbin;
double fraction;
{
	int top, bot, cut;

	cut = FRACTION*hist[peak];

	for (bot = peak - MINBIN; bot > 0; bot--)
		if (hist[bot] < cut)
			break;
	for (top = peak + MINBIN; top < (NHIST - 1); top++)
		if (hist[top] < cut)
			break;

	*botbin = bot;
	*topbin = top;
}



        /* find the stdev of all pixels which fall into the bins between
         * 'botbin' and 'topbin', inclusive, and return the number.
         * to find the stdev, fit a gaussian to the histogram values
         * in the given range and return the gaussian's "peak" and "sigma".
         *
         * Return
         *     0  if all goes well, or
         *    -1  if the standard deviation is undefined, or if we
         *             can't allocate for memory.
         */

static int
find_fitted_values(botbin, topbin, peak, sigma)
int topbin, botbin;
float *peak, *sigma;
{
	int i, ndata;
	float *x, *y, *sig, amp, center, width, chisq;

	if ((ndata = (topbin - botbin) + 1) < 3) {
		error(0, "find_fitted_values: to few points to fit gaussian to sky");
		return(0.0);
	}	
	if ((x = vector(ndata)) == NULL) {
		error(0, "find_fitted_values: can't alloc for x");
		return(0.0);
	}
	if ((y = vector(ndata)) == NULL) {
		error(0, "find_fitted_values: can't alloc for y");
		free(x);
		return(0.0);
	}
	if ((sig = vector(ndata)) == NULL) {
		error(0, "find_fitted_values: can't alloc for sig");
		free(x);
		free(y);
		return(0.0);
	}
	for (i = 0; i < ndata; i++) {
		x[i] = i;
		y[i] = hist[i + botbin];
		sig[i] = 1.0;
	}
	if (gauss_fit(x, y, sig, ndata, &amp, &center, &width, &chisq) != 0) {
		error(0, "find_fitted_values: gaussFit returns with error");
		error(0, "find_fitted_values: setting sky sigma value to zero");
		width = 0.0;
		return(width);
	}

	free(x);
	free(y);
	free(sig);

	*peak = center + botbin;
	*sigma = width;

	return(0);
}
	

	/* return a pointer to an allocated array of floats, or NULL if 
	   we couldn't allocate memory */

static float *
vector(nel)
int nel;
{
	float *p;

	if ((p = (float *) malloc(nel*sizeof(float))) == NULL) {
		fprintf(stderr, "sky: can't allocate for array of %d floats", nel);
		return(NULL);
	}
	return(p);
}