/*-------------------------------------------------------------------------*/
/**
   @file	detector.c
   @author	Yves Jung
   @date	June 2001
   @version	$Revision: 1.8 $
   @brief	All detector check functions
*/
/*--------------------------------------------------------------------------*/

/*
	$Id: detector.c,v 1.8 2002/01/11 09:34:54 ndevilla Exp $
	$Author: ndevilla $
	$Date: 2002/01/11 09:34:54 $
	$Revision: 1.8 $
*/

/*---------------------------------------------------------------------------
   								Includes
 ---------------------------------------------------------------------------*/

#include "detector.h"
#include "doubles.h"
#include "dstats.h"
#include "random.h"

/*---------------------------------------------------------------------------
   								Defines
 ---------------------------------------------------------------------------*/

#define MIN_GAIN 	0.1	
#define MAX_GAIN	10.0
#define BIG_GAIN	(double)100.0 

/*---------------------------------------------------------------------------
   							Function codes	
 ---------------------------------------------------------------------------*/


/*-------------------------------------------------------------------------*/
/**
  @brief	Create a regression map from a set of exposures.
  @param	cube_in		Input cube.
  @return	1 newly allocated array of 3 newly allocated images.
  @see		cube_create_gainmap_robust() cube_create_gainmap_proportional()

  The input is assumed to be a cube containing planes of different
  intensities (usually increasing or decreasing). Typical inputs are:
  twilight data sets, halogen lamp, or skies of different airmasses in the
  thermal regime.

  The output is a set of 3 images. The first image contains a regression
  map, i.e. for each pixel position on the detector, a curve is plotted of
  the pixel intensity in each plane against the median intensity of the
  plane. A slope is fit, and the gain factor is stored into this first
  image.

  The second image contains the y-intercepts of the slope fit. It is
  usually good to check it out in case of failures.

  The third image is expected to contain the sum of squared errors for each
  fit. In the current implementation, a NULL pointer is returned (the error
  image is not computed in this version).

  The fit is using a basic least-squares criterion, without any refinement.
  If the input points are noisy, the resulting slope is affected
  immediately. For more precise results rejecting spurious signals, see
  cube_create_gainmap_robust().

  The returned result is an array of 3 image pointers, that must be
  deallocated using free(). Each of the returned image pointers must have
  been previously deallocated (if non-NULL) using image_del().

  Example:
  \begin{verbatim}
  image_t ** slopefit ;

  slopefit = cube_create_gainmap(cube);

  ...

  if (slopefit[0]!=NULL)
	  image_del(slopefit[0]);

  if (slopefit[1]!=NULL)
	  image_del(slopefit[1]);

  if (slopefit[2]!=NULL)
	  image_del(slopefit[2]);

  free(slopefit);
  \end{verbatim}
 */
/*--------------------------------------------------------------------------*/
image_t ** cube_create_gainmap(cube_t * cube_in)
{
	image_t	**	result_array ;

	/*
	 * result_array[0] contains a pointer to the slopes image
	 * result_array[1] contains a pointer to the ordinates image
	 * result_array[2] contains a pointer to the errors image
	 */

	image_t	*	gain_map ;
	image_t	*	slopes_map ;
	double		*	sum_products ;
	double		*	sum_fluxes ;
	image_t	*	curr_img ;
	double			cur_mean,
					sum_mean,
					sq_sum_mean,
					curr_pix ;
	double			width ;
	double		*	slopes ;
	double			mean_slope ;
	int				n_plane ;
	int				curpos ;
	int				imsize ;
	double			inv_n_im ;
	int				nbpix ;

	/* Error handling : test entries	*/

	if (cube_in==NULL) return NULL ;

	imsize = cube_in->lx * cube_in->ly ;	
	inv_n_im = (double)1.0 / (double)cube_in->np ;

	/* Allocate tables	*/

	result_array = malloc(3 * sizeof(image_t*)) ;

	sum_products = (double*)malloc(imsize * sizeof(double)) ;
	sum_fluxes   = (double*)malloc(imsize * sizeof(double)) ;
	slopes       = (double*)malloc(imsize * sizeof(double)) ;
	
	sum_mean = sq_sum_mean = (pixelvalue)0.0 ;
	for (curpos=0 ; curpos<imsize ; curpos++)
		sum_products[curpos] =
		sum_fluxes[curpos] =
		slopes[curpos] = (double)0.0 ;

	/* Loop on all planes of input cube	*/
	for (n_plane=0 ; n_plane<cube_in->np ; n_plane++) {
		compute_status( 	"computing gain map", n_plane, cube_in->np, 2) ;
		curr_img = cube_in->plane[n_plane] ;

		/*
		 * Now we go: accumulate the mean value of the plane, along
		 * with the variation of each pixel.
		 * Actually, the median value is used: more robust than mean.
		 */

		cur_mean = (double)image_getmedian(curr_img) ;
		sum_mean += cur_mean ; 
		sq_sum_mean += (cur_mean * cur_mean) ; 
		
		/* Get each pixel variation in an accumulator 	*/

		for (curpos=0 ; curpos < (curr_img->lx * curr_img->ly) ; curpos++) {
			curr_pix = (double)curr_img->data[curpos] ;
			sum_products[curpos] += curr_pix * cur_mean ;
			sum_fluxes[curpos] += curr_pix ;
		}
	}

	/* Allocate for slopes map	*/
	slopes_map = image_new(cube_in->lx, cube_in->ly) ;
	if (slopes_map == NULL) {
		e_error("cannot allocate slopes map: aborting gain map creation") ;
		free(sum_products) ;
		free(sum_fluxes) ;
		free(slopes) ;
		return((image_t**)NULL) ;
	}

	gain_map = image_new(cube_in->lx, cube_in->ly) ;
	if (gain_map == NULL) {
		e_error("cannot allocate gain map: aborting gain map creation") ;
		image_del(slopes_map) ;
		free(sum_products) ;
		free(sum_fluxes) ;
		free(slopes) ;
		return((image_t**)NULL) ;
	}

	sq_sum_mean *= inv_n_im ;
	sum_mean *= inv_n_im ;
	
	width = (sq_sum_mean - (sum_mean * sum_mean)) ; 
	if (fabs((double)width) < (double)1e-4) {
		e_error("not enough variations (%g) to apply this method", width) ;
		image_del(slopes_map) ;
		image_del(gain_map) ;
		free(slopes) ;
		free(sum_products) ;
		free(sum_fluxes) ;
		return((image_t**)NULL) ;
	}
	width = (double)1.0 / width ;

	mean_slope = (double)0.0 ;
	nbpix = gain_map->lx * gain_map->ly ;
	for (curpos = 0 ; curpos < nbpix ; curpos++) {
		sum_products[curpos] *= inv_n_im ;
		sum_fluxes[curpos] *= inv_n_im ;

		slopes[curpos] = (sum_products[curpos] - sum_mean * sum_fluxes[curpos]) 
						* width ;

		mean_slope += slopes[curpos] ;

		/* Store result in slopes map	*/
		slopes_map->data[curpos] = (pixelvalue)
						(((sq_sum_mean * sum_fluxes[curpos]) -
						 (sum_mean * sum_products[curpos])) * width) ;
	}
	mean_slope /= (double)nbpix ;

	
	/* GAIN IMAGE NORMALIZATION */

	if (mean_slope != (double)0.0) {
		/* A multiplication is always faster than a division	*/
		mean_slope = (double)1.0 / mean_slope ;
		for (curpos=0 ; curpos<nbpix ; curpos++)
			gain_map->data[curpos] = (pixelvalue)(slopes[curpos] * mean_slope) ;
	}

	free(slopes) ;
	free(sum_products) ;
	free(sum_fluxes) ;

	result_array[0] = gain_map ;
	result_array[1] = slopes_map ;
	result_array[2] = NULL ;
	return result_array ;
}


/*-------------------------------------------------------------------------*/
/**
  @brief	Compute a flat-field out of a set of exposures.
  @param	twilight	Input cube.
  @return	1 newly allocated array of 3 newly allocated images.
  @see		cube_create_gainmap() cube_create_gainmap_proportional()

  The input is assumed to be a cube containing planes of different
  intensities (usually increasing or decreasing). Typical inputs are:
  twilight data sets, halogen lamp, or skies of different airmasses in the
  thermal regime.

  The output is a set of 3 images. The first image contains a regression
  map, i.e. for each pixel position on the detector, a curve is plotted of
  the pixel intensity in each plane against the median intensity of the
  plane. A slope is fit, and the gain factor is stored into this first
  image.

  The second image contains the y-intercepts of the slope fit. It is
  usually good to check it out in case of failures.

  The third image contains the sum of squared errors for each fit.

  The fit is using a robust least-squares criterion rejecting outliers.
  This is the algorithm to use with big telescopes like the VLT, which
  collect so much light that objects are actually seen in the twilight sky.
  It is also recommended to jitter the twilight acquisition in this case
  (this is what is done on ISAAC).

  The returned result is an array of 3 image pointers, that must be
  deallocated using free(). Each of the returned image pointers must have
  been previously deallocated using image_del().

  Example:
  \begin{verbatim}
  image_t ** slopefit ;

  slopefit = cube_create_gainmap_robust(cube);

  ...

  if (slopefit[0]!=NULL)
	  image_del(slopefit[0]);

  if (slopefit[1]!=NULL)
	  image_del(slopefit[1]);

  if (slopefit[2]!=NULL)
	  image_del(slopefit[2]);

  free(slopefit);
  \end{verbatim}
 */
/*--------------------------------------------------------------------------*/
image_t ** cube_create_gainmap_robust(cube_t * twilight)
{
	image_t		**	result ;
	image_t		*	gain ;
	image_t		*	intercept ;
	image_t		*	sq_err ;
	int				i, p ;
	double		*	slope ;
	double3		*	timeline ;
	double		*	plane_med ;

	if (twilight==NULL) return NULL ;

	/* Compute median for all planes */

	plane_med = malloc(twilight->np * sizeof(double));
	for (p=0 ; p<twilight->np ; p++) {
		compute_status("computing stats...", p, twilight->np, 1);
		plane_med[p] = (double)image_getmedian(twilight->plane[p]);
	}

	result = malloc(3 * sizeof(image_t*)) ;

	gain 		= image_new(twilight->lx, twilight->ly);
	intercept 	= image_new(twilight->lx, twilight->ly);
	sq_err 		= image_new(twilight->lx, twilight->ly);

	timeline = double3_new(twilight->np);

	/* Loop on all pixel positions */
	e_comment(1, "computing gains for all positions (long)...");
	for (i=0 ; i<(gain->lx*gain->ly) ; i++) {
		/* extract time line */
		for (p=0 ; p<twilight->np ; p++) {
			timeline->x[p] = plane_med[p] ;
			timeline->y[p] = (double)twilight->plane[p]->data[i] ;
		}

		/* fit slope to this time line */
		slope = fit_slope_robust(timeline);
		/* set results in output images */
		intercept->data[i] 	= slope[0] ;
		gain->data[i] 		= slope[1] ;
		sq_err->data[i]		= slope[2] ;
		free(slope);
	}
	free(plane_med);
	double3_del(timeline);

	/* Normalize output flat-field image to unity mean */
	result[0] = image_normalize(gain, NORM_MEAN);
	image_del(gain);

	result[1] = intercept ;
	result[2] = sq_err ;

	return result ;
}


/*-------------------------------------------------------------------------*/
/**
  @brief	Compute a flat-field out of a set of exposures.
  @param	twilight	Input cube.
  @return	1 newly allocated array containing 2 newly allocated images.
  @see 		cube_create_gainmap
  @see		cube_create_gainmap_robust

  The input is assumed to be a cube containing planes of different
  intensities (usually increasing or decreasing), from which any source of
  bias has been removed. Typical inputs are: twilight data sets, halogen
  lamp, or skies of different airmasses in the thermal regime. The input
  frame should have been dark-subtracted or de-biased before entering this
  function.

  The output is an array of 2 images. The first image contains a regression
  map, i.e. for each pixel position on the detector, a curve is plotted of
  the pixel intensity in each plane against the median intensity of the
  plane. A slope is fit assuming a zero y-intercept, and the gain factor is
  stored into this first image.

  The second image contains the sum of squared errors for each fit.

  The fit is using a robust slope fit criterion rejecting outliers. The
  slope of each pixel is computed in all the input planes, and only the
  median slope is stored in output.

  The returned result is an array of 2 image pointers, that must be
  deallocated using free(). Each of the returned image pointers must have
  been previously deallocated using image_del().

  Example:
  \begin{verbatim}
  image_t ** slopefit ;

  slopefit = cube_create_gainmap_proportional(cube);

  ...

  if (slopefit[0]!=NULL)
	  image_del(slopefit[0]);

  if (slopefit[1]!=NULL)
	  image_del(slopefit[1]);

  free(slopefit);
  \end{verbatim}

 */
/*--------------------------------------------------------------------------*/
image_t ** cube_create_gainmap_proportional(cube_t * twilight)
{
	image_t		**	result ;
	image_t		*	gain ;
	image_t		*	sq_err ;
	int				i, p ;
	double		*	slope ;
	double3		*	timeline ;
	double		*	plane_med ;
	int				step ;
	int				ly_step ;

	if (twilight==NULL) return NULL ;

	/* Compute median for all planes */
	plane_med = malloc(twilight->np * sizeof(double));
	for (p=0 ; p<twilight->np ; p++) {
		compute_status("computing stats...", p, twilight->np, 1);
		plane_med[p] = (double)image_getmedian(twilight->plane[p]);
	}

	result = malloc(2 * sizeof(image_t*)) ;

	gain 		= image_new(twilight->lx, twilight->ly);
	sq_err 		= image_new(twilight->lx, twilight->ly);

	timeline = double3_new(twilight->np) ;

	/* Loop on all pixel positions */
	step = 0 ;
	ly_step = 0 ;
	for (i=0 ; i<(gain->lx*gain->ly) ; i++) {
		/* Compute status message is handled here with 'step' */
		if (!step) {
			compute_status("computing gain...", ly_step, gain->ly, 1);
		}
		step++ ;
		if (step==gain->lx) {
			step=0 ;
			ly_step++ ;
		}
		/* extract time line */
		for (p=0 ; p<twilight->np ; p++) {
			timeline->x[p] = plane_med[p] ;
			timeline->y[p] = (double)twilight->plane[p]->data[i] ;
		}

		/* fit slope to this time line */
		slope = fit_proportional(timeline);
		/* set results in output images */
		gain->data[i] 		= slope[0] ;
		sq_err->data[i]		= slope[1] ;
		free(slope);
	}
	free(plane_med);
	double3_del(timeline);

	/* Normalize output flat-field image to unity mean */
	result[0] = image_normalize(gain, NORM_MEAN);
	image_del(gain);
	result[1] = sq_err ;

	return result ;
}


/*-------------------------------------------------------------------------*/
/**
  @brief    Compute the readout noise in an image.
  @param    diff    Input image, usually a difference frame.
  @param    zone    Zone where the readout noise is to be computed.
  @param    noise   Output parameter: noise in the frame.
  @param    error   Output parameter: error on the noise.
  @return   int 0 if Ok, -1 otherwise.

  This function is meant to compute the readout noise in a frame by
  means of a MonteCarlo approach. The input is a frame, usually a
  difference between two frames taken with the same settings for the
  acquisition system, although no check is done on that, it is up to
  the caller to feed in the right kind of frame.

  The provided zone is an array of four integers specifying the zone
  to take into account for the computation. The integers specify
  ranges as xmin,xmax,ymin,ymax, where these coordinates are given
  in the FITS notation (x from 1 to lx, y from 1 to ly and bottom to
  top). Specify NULL instead of an array of four values to use the
  whole frame in the computation.

  The algorithm will create typically 100 9x9 windows on the frame,
  scattered optimally using a Poisson law. In each window, the standard
  deviation of all pixels in the window is computed and this value is
  stored. The readout noise is the median of all computed standard
  deviations, and the error is the standard deviation of the standard
  deviations.

  Both values (noise and error) are returned by modifying a passed
  double. If you do not care about the error, pass NULL.
 */
/*--------------------------------------------------------------------------*/
#define RON_REG_HLX     4
#define RON_REG_HLY     4
#define RON_SAMPLES     100
int image_readout_noise(
	    image_t *   diff,
	    int     *   zone_def,
	    double  *   noise,
	    double  *   error)
{
    double3     *   sample_reg ;
    int				rect[4] ;
    image_stats *   vig_stats ;
    double      *   rms_list ;
    double          rms_median ;
    double          rms_error ;
    int             zone[4];
    int             i ;

    if (diff==NULL || noise==NULL) return -1 ;

    /* Generate RON_SAMPLES window centers in the image */
    if (zone_def!=NULL) {
        rect[0] = zone_def[0] + RON_REG_HLX + 1 ; /* xmin */
        rect[1] = zone_def[1] - RON_REG_HLX - 1 ; /* xmax */
        rect[2] = zone_def[2] + RON_REG_HLY + 1 ; /* ymin */
        rect[3] = zone_def[3] - RON_REG_HLY - 1 ; /* ymax */
    } else {
        rect[0] = RON_REG_HLX + 1 ;	            /* xmin */
        rect[1] = diff->lx - RON_REG_HLX - 1 ;	/* xmax */
        rect[2] = RON_REG_HLY + 1 ;	            /* ymin */
        rect[3] = diff->ly - RON_REG_HLY - 1 ;	/* ymax */
    }

    if ((rect[0]>=rect[1]) || (rect[2]>=rect[3])) {
        e_error("in readout noise: invalid region def");
        return -1 ;
    }

    /* Generate n+1 regions, because the first region is always at (0,0) */
    /* and it would bias the measurement. */
    sample_reg = generate_poisson_points(rect, RON_SAMPLES+1, RON_SAMPLES+1) ;

    /* Now, for each window center, extract a vignette and compute the */
    /* signal RMS in it. Store this rms into a table. */
    rms_list = malloc(RON_SAMPLES * sizeof(double));
    for (i=0 ; i<RON_SAMPLES ; i++) {
        zone[0] = (int)sample_reg->x[i+1] - RON_REG_HLX ;
        zone[1] = (int)sample_reg->x[i+1] + RON_REG_HLX ;
        zone[2] = (int)sample_reg->y[i+1] - RON_REG_HLY ;
        zone[3] = (int)sample_reg->y[i+1] + RON_REG_HLY ;
        if ((vig_stats = image_getstats_opts(diff,
                                         NULL,
                                         NULL,
                                         zone)) == NULL) {
            e_error("computing stats on vignette");
            free(rms_list);
            double3_del(sample_reg);
            return -1 ;
        }
        rms_list[i] = (double)vig_stats->stdev;
        free(vig_stats);
    }
    double3_del(sample_reg);

    /* The final computed RMS is the median of all values.  */
    rms_median = double_median(rms_list, RON_SAMPLES);
    (*noise) = rms_median ;
    if (error!=NULL) {
        /* The error is the rms of the rms */
        rms_error = double_rms(rms_list, RON_SAMPLES);
        (*error) = rms_error ;
    }
    free(rms_list);
    return 0 ;
}
#undef RON_REG_HLX
#undef RON_REG_HLY
#undef RON_SAMPLES


/*-------------------------------------------------------------------------*/
/**
  @brief	Compute the linearity of the detector
  @param	in	input cube
  @param	dit	list of DIT values
  @return	cube with 4 images (3 coeffs & rms)
 */
/*--------------------------------------------------------------------------*/
cube_t * detector_linearity_fit(
		cube_t	*	in, 
		double	*	dit)
{
    cube_t      *   fitres ;
    image_t     *   timeline ;
    int             i, j, k ;
    int             pos ;
    int             npix ;
    double          f ;
    double          y, err, sq_err ;
    matrix      *   mx, 
				*	ma, 
				*	mb ;

    if (in==NULL || dit==NULL) return NULL ;

    /* Allocate 4 images to store the results */
    fitres = cube_new(in->lx, in->ly, 4);
    for (k=0 ; k<4 ; k++) {
        fitres->plane[k] = image_new(in->lx, in->ly);
    }

    /* Loop over all pixels */
    npix = in->lx * in->ly ;

    /* Create matrices */
    ma = matrix_new(3, in->np);
    mb = matrix_new(1, in->np);

    /*
     * Sorry for the double loop, it is not truly needed.
     * The only interest is to put a compute_status that
     * does not print out too often, otherwise it slows down
     * the already heavy computation.
     */
    for (j=0 ; j<in->ly ; j++) {
        compute_status("fitting polynomial...", j, in->ly, 1);
        for (i=0 ; i<in->lx ; i++) {
            pos = i + j * in->lx ;

            /* Extract time line */
            timeline = cube_get_z(in, pos);
            if (timeline==NULL) {
                e_error("extracting time line at pos %d: aborting", pos);
                cube_del(fitres);
                return NULL ;
            }

            /* Fill in matrices */
            for (k=0 ; k<in->np ; k++) {
                f = (double)timeline->data[k] ;

                ma->m[k]            = f ;
                ma->m[k+in->np]     = f*f ;
                ma->m[k+2*in->np]   = f*f*f ;

                mb->m[k] = dit[k] ;
            }

            /* Solve least-squares */
            mx = matrix_leastsq(ma,mb) ;
            if (mx == NULL) {
                fitres->plane[0]->data[pos] = 0 ;
                fitres->plane[1]->data[pos] = 0 ;
                fitres->plane[2]->data[pos] = 0 ;
                fitres->plane[3]->data[pos] = 0 ;
            } else {
                /* Store results in a, b, c images */
                fitres->plane[0]->data[pos] = (pixelvalue)mx->m[0] ;
                fitres->plane[1]->data[pos] = (pixelvalue)mx->m[1] ;
                fitres->plane[2]->data[pos] = (pixelvalue)mx->m[2] ;

                /* Goodness of fit */
                sq_err=0 ;
                for (k=0 ; k<in->np ; k++) {
                    /* Computed by model */
                    y = mx->m[0] * ma->m[k] +
                        mx->m[1] * ma->m[k+in->np] +
                        mx->m[2] * ma->m[k+2*in->np] ;
                    /* Error between model and reality */
                    err = y - mb->m[k] ;
                    sq_err += err * err ;
                }
                sq_err /= (double)in->np ;
                fitres->plane[3]->data[pos] = (pixelvalue)sq_err ;
                matrix_del(mx);
            }
            image_del(timeline);
        }
    }

    /* Delete matrices */
    matrix_del(ma);
    matrix_del(mb);

    return fitres ;
}