/*---------------------------------------------------------------------------
   
   File name 	:	jitter.c
   Author 		:	N. Devillard
   Created on	:	July 2001
   Description	:	CONICA/NAOS jitter data reduction.

 *--------------------------------------------------------------------------*/

/*
	$Id: jitter.c,v 1.11 2002/02/19 14:20:03 yjung Exp $
	$Author: yjung $
	$Date: 2002/02/19 14:20:03 $
	$Revision: 1.11 $
*/

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

#include <stdio.h>
#include <stdlib.h>
#include <string.h>

#include "eclipse.h"
#include "conicap_lib.h"
#include "detect_ks.h"

#include "gnuplot_i.h"

/*---------------------------------------------------------------------------
   								Define
 ---------------------------------------------------------------------------*/
	
#define	NB_POINTS_XCORR	1
#define MEDIAN_SKY_METHOD 1
#define AVG_SKY_METHOD 2
#define LOW_RATE 0.1
#define HIGH_RATE 0.1

#define FIFTY_HZ_HW_SMOOTH    20
#define FIFTY_HZ_REJECTED_LOW_PIXELS    0
#define FIFTY_HZ_REJECTED_HIGH_PIXELS 900
#define FIFTY_HZ_THRESHOLD 5.0

/*---------------------------------------------------------------------------
							Private functions
 ---------------------------------------------------------------------------*/

static int conica_jitter_engine(char *, char *, int, int, int) ;
static int conica_remove_sky(cube_t **, int, double, double, int) ;
static framelist * conica_loadlist(char *) ;
static double3 * conica_getoffsets(framelist*) ;
static int conica_remove_50Hz(cube_t **) ;

/*---------------------------------------------------------------------------
  							Function codes
 ---------------------------------------------------------------------------*/
int conica_jitter_main(void * dict)
{
    dictionary  *   d ;

    char            argname[10] ;
    char        *   name_i ;
    char        *   name_o ;

    int             skysub_activate ;
    int             mfilt_activate ;
	int				remove_50hz ;

    int             nfiles ;
    int             errors ;
    int             i ;

    /* Initialize */
    skysub_activate = 1 ;
    mfilt_activate  = 1 ;
	remove_50hz     = 1 ;

    d = (dictionary*)dict ;
    /* Get options */
    skysub_activate = dictionary_getint(d, "arg.sky", 1) ;
    mfilt_activate  = dictionary_getint(d, "arg.filter", 1) ;
	remove_50hz     = dictionary_getint(d, "arg.fifty_hz", 1) ;
	
    /* Get input/output file names */
    nfiles = dictionary_getint(d, "arg.n", -1) ;
    if (nfiles<0) {
        e_error("missing input file name(s): aborting");
        return -1 ;
    }

    /* Loop on input file names */
    errors = 0 ;
    for (i=1 ; i<nfiles ; i++) {
        sprintf(argname, "arg.%d", i);
        name_i = dictionary_get(d, argname) ;
        name_o = dictionary_get(d, "arg.output") ;
        if (name_o == NULL) name_o = strdup(get_rootname(get_basename(name_i)));
        else name_o = strdup(get_rootname(name_o)) ;

        /*
         * Once command-line options have been cleared out, call the main
         * computing function.
         */
        errors += conica_jitter_engine(name_i, name_o, skysub_activate, 
				mfilt_activate, remove_50hz) ;
        free(name_o) ;
    }
    return errors ;
}

static int conica_jitter_engine(
		char 	* 	name_i,
		char 	*	name_o,
		int			skysub_activate,
		int			mfilt_activate,
		int			remove_50hz)
{
	framelist * input_list ;
	cube_t	*	jit ;
	image_t	*	combined ;
	char	*	dark_name ;
	char	*	flat_name ;
	char	*	bp_name ;
	int         sky_method ;
	double3	*	offsets ;
	double3	*	xcorr_obj ;
	double3	*	first_brightest ;
	double3	*	brightest ;
	double3	*	xcorr_est ;
    int         i ;
    char        outname[FILENAMESZ] ;
	
	/* All settings temporarily hardcoded here. */
	dark_name = NULL ;
	flat_name = NULL ;
	bp_name = NULL ;
	sky_method = MEDIAN_SKY_METHOD ;

	/* Load input cube */
	e_comment(0, "loading input cube %s", name_i);
	if ((input_list = conica_loadlist(name_i)) == NULL) {
		e_error("loading %s: aborting", name_i);
		return 1 ;
	}
	if ((jit = cube_load_strings(input_list->name, input_list->n)) == NULL) {
		e_error("loading cube: aborting");
		framelist_del(input_list);
		return 1 ;
	}
	
	/* Apply dark subtraction / flatfield division if needed */
	e_comment(0, "applying dark/flat/badpix correction");
	conica_ff_dark_badpix_handling(&jit, flat_name, dark_name, bp_name);

	/* Filter out sky */
	if (skysub_activate) {
		e_comment(0, "applying sky subtraction") ;
		if (conica_remove_sky(&jit, 
					sky_method,
					LOW_RATE,
					HIGH_RATE,
					1) == -1) {
			e_warning("cannot remove sky - skipping") ;
		} 
	}

	/* Remove 50 Hz */
	if (remove_50hz) {
		e_comment(0, "Remove 50 Hz") ;
		if (conica_remove_50Hz(&jit) == -1) {
			e_warning("cannot remove 50 Hz - skipping") ;
		}
	}
	
	/* Get the NB_POINTS_XCORR brightest points in the first image */
	xcorr_obj = detected_ks_brightest_stars(jit->plane[0],
							NB_POINTS_XCORR,
							DETECTED_KAPPA) ;
	
    /* Get offset estimates from headers if possible */
    e_comment(0, "retrieving offsets from header...");
    xcorr_est = conica_getoffsets(input_list);
	framelist_del(input_list);
    if (xcorr_est == NULL) {
        e_comment(0, "using object position as offset estimate");
        xcorr_est = double3_new(jit->np);
		xcorr_est->x[0] = 0.0 ;
		xcorr_est->y[0] = 0.0 ;
		first_brightest = detected_ks_brightest_stars(jit->plane[0],
				1, DETECTED_KAPPA) ;
		for (i=1 ; i<jit->np ; i++) {
			brightest = detected_ks_brightest_stars(jit->plane[i], 
					1, DETECTED_KAPPA) ;
            xcorr_est->x[i] = brightest->x[0] - first_brightest->x[0];
            xcorr_est->y[i] = brightest->y[0] - first_brightest->y[0];
			double3_del(brightest) ;
		}
		double3_del(first_brightest) ;
	}

	/* Apply median filter to all frames to remove hot pixels */
	if (mfilt_activate) {
		e_comment(0, "applying median filter to all frames...");
		cube_filter(jit, "median", NULL);
	}

	/* Subtract the first offsets to the others */
	for (i=xcorr_est->n - 1 ; i>=0 ; i--) {
		xcorr_est->x[i] -= xcorr_est->x[0] ; 
		xcorr_est->y[i] -= xcorr_est->y[0] ; 
	}

	/* Apply refining with cross-correlation */
    e_comment(0, "applying cross-correlation pass...") ;
	if ((offsets = xcorr_with_objs(jit, jit->plane[0], 
					xcorr_est, xcorr_obj, 35, 35, 55, 55)) == NULL) {
		e_error("cannot apply the x-correlation") ;
    	double3_del(xcorr_obj) ;
    	double3_del(xcorr_est) ;
		cube_del(jit) ;
		return -1 ;
	}
    double3_del(xcorr_obj) ;
    double3_del(xcorr_est) ;
	
    e_comment(0, "refined offsets") ;
	for (i=0 ; i<offsets->n ; i++) {
		e_comment(1, "% 4.2f % 4.2f %g\n", offsets->x[i], offsets->y[i],
				offsets->z[i]) ;  
	}

	/* Mupltiply offsets by -1 */
	for (i=0 ; i<offsets->n ; i++) {
		offsets->x[i] *= -1 ;
		offsets->y[i] *= -1 ;
	}

	/* Apply shift and add with these measurements */
    e_comment(0, "applying shift and add") ;
	combined = cube_shiftandadd(jit, offsets, "default", 0, 0) ;
	double3_del(offsets) ;

	/* Save results */
	e_comment(0,"saving results") ;
    sprintf(outname, "%s.fits", name_o) ; 
	image_save_fits(combined, outname, BPP_DEFAULT) ;
	image_del(combined) ;
	cube_del(jit) ;
	e_comment(0,"done.\n") ;

	return 0 ;
}

static framelist * conica_loadlist(char * name_i)
{
	framelist * flist ;
	framelist * objs ;

	/* Load input framelist */
	flist = framelist_load(name_i);
	if (flist==NULL) {
		e_error("loading list %s", name_i);
		return NULL ;
	}
	/* Select only frames with DPR.TYPE=='OBJECT' */
	objs = framelist_select_tokenget(flist, "OBJECT", conica_get_dpr_type);
	framelist_del(flist);
	if (objs==NULL) {
		e_error("no object frame found in input");
	}
	return objs ;
}

static double3 * conica_getoffsets(framelist * flist)
{
    char        *   xval,
                *   yval ;
    double3     *   offs ;
    int             i ;
    int             header_ok ;

    offs = double3_new(flist->n);
    header_ok = 1 ;
    for (i=0 ; i<flist->n ; i++) {
        if ((xval = conica_get_cumoffsetx(flist->name[i])) == NULL) {
			e_warning("no offset in headers");
			header_ok=0 ;
			break ;
		} else offs->x[i]=(double)atof(xval);
		
        if ((yval = conica_get_cumoffsety(flist->name[i])) == NULL) {
			e_warning("no offset in headers");
			header_ok=0 ;
			break ;
		} else offs->y[i]=(double)atof(yval);

    }
    if (!header_ok) {
        double3_del(offs) ;
        offs=NULL ;
    }
    return offs ;
}

static int conica_remove_sky(
		cube_t	**	incube,
		int        	method,
		double 		low_rate,
		double 		high_rate,
		int 		output_sky)
{
	image_t	*	sky ;
	int 		low ;
	int 		high ;
  
	/* At least 3 planes for sky correction */
	if ((*incube)->np < 3) {
		e_warning("less than 3 planes - no sky subtraction") ;
		return -1 ;
	}

	/* Set low and high */
	low = (int)(low_rate * (*incube)->np) ;
	high = (int)(high_rate * (*incube)->np) ;

	/* Test which method is used */
	switch (method) {
		case AVG_SKY_METHOD:
			if ((sky=cube_avg_reject(*incube, low, high)) == NULL) {
				e_warning("cannot estimate the sky") ;
				return -1 ;
			}
			break ;
		case MEDIAN_SKY_METHOD:
			if ((sky=cube_avg_medreject(*incube, low, high)) == NULL) {
				e_warning("cannot estimate the sky") ;
				return -1 ;
			}
			break ;
		default:
			e_warning("unrecognized method - skip") ;
			return -1 ;
			break ;
	}

	/* At this point, the sky is computed */
	if (output_sky) {
		image_save_fits(sky, "sky.fits", BPP_DEFAULT) ;
	}
  
	/* Subtract the sky */
	cube_sub_im(*incube, sky) ;

	/* Free and return */
	image_del(sky) ;
	return 0 ;
}

static int conica_remove_50Hz(
		cube_t	**	incube)
{
	image_t 	*	collapsed ;
	image_t 	*	highfreq ;
	pixelvalue	* 	lowfreq ;
	int				i, j ;

	for (i=0 ; i<(*incube)->np ; i++) {
		/* Collapse the current image */
		collapsed = image_collapse_median((*incube)->plane[i], 1, 
				FIFTY_HZ_REJECTED_LOW_PIXELS, 
				FIFTY_HZ_REJECTED_HIGH_PIXELS);
	
		/* Extract the low frequency signal */
		lowfreq = function1d_median_smooth(collapsed->data, collapsed->ly, 
				FIFTY_HZ_HW_SMOOTH) ;

		/* Subtract the lowfreq to the collapsed image */
		/* Set to 0 everything above 5 or under -5 in high frequency sig */
		highfreq = image_copy(collapsed) ;
		image_del(collapsed) ;
		for (j=0 ; j<highfreq->ly ; j++) {
			highfreq->data[j] -= lowfreq[j] ;
			if (fabs(highfreq->data[j]) > FIFTY_HZ_THRESHOLD) 
				highfreq->data[j] = 0.0 ;
		}
		free(lowfreq) ;	
		
		/* Correct the input cube */
		image_sub_1d_local((*incube)->plane[i], highfreq) ;
		image_del(highfreq) ;
	}

	return 0 ;
}