# include <math.h>
# include <stdio.h>
# include <stdlib.h>
# include <string.h>
# include <float.h>
# include <xtables.h>	/* For the IRAF_INDEF macros */
# include "pedsub.h"

/* Global variables for use by fitting routines */
extern float A[11];
extern int MA;
float A[11];
int   MA;

/*   FINDBESTPED  --  Find best pedestal values for a NICMOS image.
**
**	Based on GETBESTPED routine from Roeland van der Marel's
**	"unpedestal" program.
**
**	This routine seeks the best pedestal (DC signal offset) value
**	for each quadrant of a NICMOS image by iteratively subtracting
**	trial pedestal values and finding the value that results in the
**	minimum spread in residual pixel values. A rough guess of the best
**	pedestal value in each quadrant is determined using the Numerical
**	Recipes BRENT routine. If asked for by the user, a further
**	refinement to the pedestal value may be obtained. The refined value
**	is obtained by computing the residual pixel spread for a range
**	of trial pedestal values, centered about the rough guess value,
**	fitting a polynomial to the resulting range of pedestal vs. spread
**	values, and again using the BRENT routine to find the minimum of
**	the fitted polynomial.
**
**	Each time a trial pedestal value is subtracted from an image
**	quadrant, there are options allowing for the resulting image to
**	be filtered in some way before computing the spread in pixel values.
**	These options include no filtering, median filtering, and unsharp
**	masking (unfiltered minus-median filtered image). The filtering
**	options are useful for removing signal from sources in the image,
**	which can bias the measured pixel spread.
**
**	When the pedestal values are subtracted from an image quadrant, it
**	is actually the pedestal times the flatfield value for each pixel
**	that is subtracted, because the pedestal signal is modulated by
**	the flatfield in the calibrated images that are used as input to
**	this routine.
**
**
**   Revision history:
**   ----------------
**   H. Bushouse	28-June-1999	Implementation
**
*/

/* This structure stores data used by the Brent routine. */
typedef struct {
	PedInfo *info;
	SingleNicmosGroup *input;
	SingleNicmosGroup *Flatim;
} PedData;
extern PedData peddata;
PedData peddata;

/*  ASPREAD  -  This function exists only to enable the Num. Rec. Brent
**		routine	to access the appropriate data needed for the pixel
**		spread computation, such as the input image, the flatfield
**		image, and some general task parameters. It must be scanned
**		by the compiler *before* the function that actually calls
**		the Brent routine.
*/

float aspread (float pedval) {

	float spread (PedInfo *, SingleNicmosGroup *, SingleNicmosGroup *,
		      float);

	return (spread (peddata.info, peddata.input, peddata.Flatim, pedval));
}

/*  FINDBESTPED  -  Main algorithm.
*/

int findBestPed (PedInfo *info, SingleNicmosGroup *input,
		 SingleNicmosGroup *Flatim, FILE *fp, Bool verbose) {

/* Arguments:
**	info	io: task info structure
**	input	io: input image
**	Flatim	 i: flatfield image
**	fp	 i: log file pointer
**	verbose	 i: verbose output switch
*/

	/* Local variables */
	float mean, median, mode, stdv,	/* image statistics */
	      min, max;
	float spreadmin;		/* min image spread value */
	float bestped;			/* best pedestal value */
	float ldata[4][12];		/* log file data storage */

	/* Function declarations */
	int  n_stats (SingleNicmosGroup *, int, int, int, int, float, float,
		      short, float *, float *, float *, float *, float *,
		      float *);
	float brent (float ax, float bx, float cx, float (*f)(float),
		     float tol, float *xmin);
	int RefinePed (PedInfo *, float, float ldata[4][12], FILE *fp, Bool);
	void transform (PedInfo *, SingleNicmosGroup *, SingleNicmosGroup *,
		       float, int, SingleNicmosGroup *);

	/* Retrieve the image exposure time for later use. */
	if (info->DoRefine) {
	    if (getKeyF (input->globalhdr, "EXPTIME", &(info->exptime))) {
		sprintf (MsgText, "Can't read EXPTIME keyword from %s",
			 input->filename);
		n_error (MsgText);
		return (1);
	    }
	}

	/* Store data in the spread routine structure */
	peddata.info   = info;
	peddata.input  = input;
	peddata.Flatim = Flatim;

	/* Loop over the image quadrants */
	for (info->quad = 0; info->quad < 4; info->quad++) {

	     /* Compute the mode and median of the image quadrant,
	     ** excluding flagged pixels */
	     if (n_stats (input, info->qx1[info->quad], info->qx2[info->quad],
			  info->qy1[info->quad], info->qy2[info->quad],
			  -FLT_MAX, FLT_MAX, info->BitMask, &mean, &median,
			  &mode, &stdv, &min, &max))
		 return (1);

	     /* Save info for log file, if requested */
	     if (info->KeepLog) {
		 ldata[info->quad][0]  = mode;
		 ldata[info->quad][1]  = median;
	     }

	     /* Call the Num. Rec. BRENT routine to find the pedestal value
	     ** that produces the minimum image spread */
	     bestped = 0.0;
	     spreadmin = brent (-10.0*mode, mode, 10.0*mode, aspread, 4*FTOL,
				&bestped);

	     /* Save the best pedestal value */
	     info->PedValue[info->quad] = bestped;

	     if (verbose) {
	         sprintf (MsgText,
			  "  Quad %d: Quick pedestal value = %8.6f",
			  info->quad+1, info->PedValue[info->quad]);
		 n_message (MsgText);
	     }

	     /* If requested, refine the best pedestal value that was
	     ** just determined */
	     if (info->DoRefine) {

		 if (RefinePed (info, bestped, ldata, fp, verbose))
		     return (1);

	     } else {
		 ldata[info->quad][2] = 0.0;
		 ldata[info->quad][7] = 0.0;
		 ldata[info->quad][8] = 0.0;
		 ldata[info->quad][9] = 0.0;
		 ldata[info->quad][10] = 0.0;
		 ldata[info->quad][11] = 0.0;
		 if (verbose) {
		     sprintf (MsgText, "\n");
		     n_message (MsgText);
		 }
	     }

	     /* Subtract the best pedestal value from the quadrant */
	     transform (info, input, Flatim, info->PedValue[info->quad], NONE,
			input);

	     /* Compute the mode and median of the corrected image quadrant,
	     ** excluding flagged pixels */
	     if (n_stats (input, info->qx1[info->quad], info->qx2[info->quad],
			  info->qy1[info->quad], info->qy2[info->quad],
			  -FLT_MAX, FLT_MAX, info->BitMask, &mean, &median,
			  &mode, &stdv, &min, &max))
		 return (1);

	     /* Save diagnostic infomation for log file, if requested */
	     if (info->KeepLog) {
		 ldata[info->quad][3]  = mode;
		 ldata[info->quad][4]  = median;
		 ldata[info->quad][5]  = bestped;
		 ldata[info->quad][6]  = spreadmin;
	     }

	} /* End of loop over quadrants */
	
	/* Write diagnostic information to log file, if requested */
	if (info->KeepLog) {
	    fprintf (fp, "\n#quad      mode    median  pedestal  new mode");
	    fprintf (fp, "   new med ped_bestD  spr_minD  ped_lowP");
	    fprintf (fp, " ped_highP  spr_minP  Pfit_rms  Pfit_max\n");
	    for (info->quad=0; info->quad<4; info->quad++) {
		 fprintf (fp,
		 "#%4d %9.6f %9.6f %9.6f %9.6f %9.6f %9.6f %9.6f %9.6f %9.6f %9.6f %9.6f %9.6f\n",
		 info->quad+1, ldata[info->quad][0], ldata[info->quad][1],
		 ldata[info->quad][2], ldata[info->quad][3], ldata[info->quad][4],
		 ldata[info->quad][5], ldata[info->quad][6], ldata[info->quad][7],
		 ldata[info->quad][8], ldata[info->quad][9], ldata[info->quad][10],
		 ldata[info->quad][11]);
	    }
	}

	/* Successful return */
	return (0);
}

/*  REFINEPED  -  Refine the pedestal value.
**
**    This routine will refine the initial pedestal value determined
**    for a NICMOS image quadrant by computing the pixel spread over a
**    range of pedestal values, centered on the initial guess, fitting
**    a polynomial to the resulting pedestal vs. spread relation, and
**    finding the minimum of that function.
*/

int RefinePed (PedInfo *info, float bestped, float ldata[4][12], FILE *fp,
	       Bool verbose) {

/* Arguments:
**	info	io: task info structure
**	bestped	 i: initial guess best pedestal value
**	ldata	 o: log file data
**	fp	 i: log file pointer
**	verbose	 i: verbose output switch
*/

	/* Local variables */
	int   i;			/* loop index */
	int   Nsamp;			/* number of samples in fit */
	float chisq;			/* chi-square of fit */
	float pedlow, pedhigh;		/* pedestal error limits */
	float pedstep;			/* pedestal search step size */
	float spfit;			/* spread fit value */
	float res, resmax, resrms;	/* fit residuals */
	float spreadpolymin;		/* min spread value from fit */
	float bestpedpoly;		/* best pedestal value from poly fit */
	float *pedval, *spreadval,	/* pedestal & spread fit arrays */
	      *spreadsig;
	float **u, **v, *w;		/* poly fit work arrays */
	int *ia;
	float **covar;

	/* Function declarations */
	float aspread (float pedval);
	float spoly (float pedval);
	void  fpoly (float x, float p[], int np);
	void  svdfit (float x[], float y[], float sig[], int ndata, float a[],
		      int ma, float **u, float **v, float *w, float *chisq,
		      void (*funcs)(float, float[], int));
	float brent (float ax, float bx, float cx, float (*f)(float),
		     float tol, float *xmin);
	float rtbis(float (*func)(float), float soff, float x1, float x2,
		    float xacc);
	float *vector (int nl, int nh);
	void free_vector (float *v, int nl, int nh);
	int *ivector (int nl, int nh);
	void free_ivector (int *v, int nl, int nh);
	float **matrix (int nrl, int nrh, int ncl, int nch);
	void free_matrix (float **m, int nrl, int nrh, int ncl, int nch);

	/* Initialize the fitting params */
	Nsamp = info->Nrefine;
	MA    = info->Morder + 1;

	/* Allocate memory for the fitting arrays */
	pedval    = (float *)calloc(Nsamp+1, sizeof(float));
	spreadval = (float *)calloc(Nsamp+1, sizeof(float));
	spreadsig = (float *)calloc(Nsamp+1, sizeof(float));
	u = matrix (1, Nsamp, 1, MA);
	v = matrix (1, MA, 1, MA);
	w = vector (1, MA);
	covar = matrix (1, MA, 1, MA);
	ia = ivector (1, MA);
	for (i = 1; i <= MA; i++)
	     ia[i] = 1;

	/* Set the pedestal value step size */
	if (info->PedStep == (float)IRAF_INDEFR || info->PedStep == 0) {
	    pedstep = MAX (0.5/info->exptime, 0.03*fabs(bestped));
	} else
	    pedstep = info->PedStep;

	/* Compute image spread over range of pedestal values */
	for (i = 1; i <= Nsamp; i++) {
	     pedval[i]    = bestped + pedstep*(i-(info->Nrefine/2)-1);
	     spreadval[i] = aspread(pedval[i]);
	     spreadsig[i] = spreadval[i]/10.0;
	}

	/* Fit a polynomial to the computed spread values;
	** The fit is performed using the Num. Recipes SVDFIT routine. */
	svdfit (pedval, spreadval, spreadsig, Nsamp, A, MA, u, v, w, &chisq,
		fpoly);

	/* Write log file header line, if necessary */
	if (info->KeepLog) {
	    fprintf (fp,
		     "#quad   ped_trial      spread  spread_fit    residual\n");
	}

	/* Calculate the rms and max residual of the poly fit and
	** print this information to the log file (if requested) */
	resmax = -FLT_MAX;
	for (resrms=0, i = 1; i <= Nsamp; i++) {
	     spfit = spoly (pedval[i]);
	     res = spreadval[i] - spfit;
	     if (fabs(res) > resmax) resmax = fabs(res);
	     resrms += res*res;

	     if (info->KeepLog) {
		 fprintf (fp, "%5d %11.7g %11.7g %11.7g %11.7g\n",
			  info->quad+1, pedval[i], spreadval[i], spfit, res);
	     }
	}
	resrms = sqrt (resrms/(float)Nsamp);

	/* Find the minimum in the polynomial fit using the
	** Num. Rec. BRENT routine */
	bestpedpoly = 0.0;
	spreadpolymin = brent (pedval[1], bestped, pedval[Nsamp], spoly, FTOL,
			&bestpedpoly);

	/* Save this as the best pedestal value */
	info->PedValue[info->quad] = bestpedpoly;

	/* Report this value */
	if (verbose) {
	    sprintf (MsgText, "; Refined value = %8.6f\n",
		     info->PedValue[info->quad]);
	    n_message (MsgText);
	}

	/* Find the points where spreadpoly = spreadpolymin + resmax */
	if (spoly(pedval[1]) <= spreadpolymin + resmax)
	    pedlow = pedval[1];
	else
	    pedlow = rtbis (spoly, spreadpolymin+resrms, pedval[1],
			    bestpedpoly, FTOL);
	if (spoly(pedval[Nsamp]) <= spreadpolymin + resmax)
	    pedhigh = pedval[Nsamp];
	else
	    pedhigh = rtbis (spoly, spreadpolymin+resrms, bestpedpoly,
			     pedval[Nsamp], FTOL);

	/* Save diagnostic data for the log file, if requested */
	if (info->KeepLog) {
	    ldata[info->quad][2]  = bestpedpoly;
	    ldata[info->quad][7]  = pedlow;
	    ldata[info->quad][8]  = pedhigh;
	    ldata[info->quad][9]  = spreadpolymin;
	    ldata[info->quad][10] = resrms;
	    ldata[info->quad][11] = resmax;
	}

	/* Free local arrays */
	free (pedval);
	free (spreadval);
	free (spreadsig);
	free_matrix (u, 1, Nsamp, 1, MA);
	free_matrix (v, 1, MA, 1, MA);
	free_vector (w, 1, MA);
	free_ivector (ia, 1, MA);
	free_matrix  (covar, 1, MA, 1, MA);

	return (0);

}