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

   File name    :	twflat.c
   Author       :	N. Devillard
   Created on   :	January 2002
   Description  :	CONICA imaging flat-field creation from twilight images

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

/*

	$Id: twflat.c,v 1.2 2002/01/23 13:59:42 ndevilla Exp $
	$Author: ndevilla $
	$Date: 2002/01/23 13:59:42 $
	$Revision: 1.2 $

*/

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

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

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

#define LO_THRESH_BADPIX	0.5
#define HI_THRESH_BADPIX	2.0

/*---------------------------------------------------------------------------
                            Function prototypes
 ---------------------------------------------------------------------------*/
static int conica_twilight_save(char *, char *, int, image_t **, int, int, 
		int, int, double, double) ;
static int conica_twflat_engine(char *, char	*, char	*, int, int, int, int, 
		double, double) ;

/*----------------------------------------------------------------------------
							Main code
 ----------------------------------------------------------------------------*/
int conica_twflat_main(void * dict)
{
	dictionary	*	d ;
	
	int				error_map_flag ;
	int				pixmap_flag ;
	int				intercepts_flag ;
	int				proportional_flag ;
	char		*	dark_name;
	double			lo_thresh ;
	double			hi_thresh ;

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

	char        *   tmp_string ;
	int             items ;
	int             errors ;
	int             i ;

	d = (dictionary*)dict ;
	/* Get options */
	tmp_string = dictionary_get(d, "arg.threshold") ; 
	if (tmp_string == NULL) {
		lo_thresh = LO_THRESH_BADPIX ;
		hi_thresh = HI_THRESH_BADPIX ;
	} else {
		items = sscanf(tmp_string, "%lg %lg", &lo_thresh, &hi_thresh) ;
		if (items != 2) {
			lo_thresh = LO_THRESH_BADPIX ;
			hi_thresh = HI_THRESH_BADPIX ;
		}
	}
	intercepts_flag   = dictionary_getint(d, "arg.intercepts", 0);
	error_map_flag    = dictionary_getint(d, "arg.errmap", 0);
	pixmap_flag       = dictionary_getint(d, "arg.pixmap", 0);
	proportional_flag = dictionary_getint(d, "arg.prop", 0);
	dark_name         = dictionary_get(d, "arg.dark") ;

	/* 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_twflat_engine(name_i,
									name_o,
									dark_name,
									intercepts_flag,
									error_map_flag,
									pixmap_flag,
									proportional_flag,
									lo_thresh,
									hi_thresh) ;

		free(name_o) ;
	}
	return errors ;
}


/*---------------------------------------------------------------------------
  							Function codes
 ---------------------------------------------------------------------------*/
static int conica_twflat_engine(
	char	*	name_in,
	char	*	name_out,
	char	*	dark_name,
	int			output_intercepts,
	int			output_errmap,
	int			output_pixmap,
	int			prop_fit,
	double		lo_thresh,
	double		hi_thresh)
{
	cube_t		*	in ;
	image_t	**	results ;
	image_t	*	dark_frame ;
	int				nplanes ;

	/* Load input cube */
	e_comment(0, "---> Part 1 of 4: loading input cube");
	in = cube_load(name_in) ;
	if (in == NULL)  {
		e_error("in loading cube [%s]: aborting", name_in) ;
		return -1 ;
	}
	nplanes = in->np ;

	/* Apply dark correction to all planes if requested */
	e_comment(0, "---> Part 2 of 4: subtracting dark");
	if (dark_name!=NULL) {
		e_comment(1, "loading dark frame [%s]", dark_name);
		dark_frame = image_load(dark_name);
		if (dark_frame==NULL) {
			e_error("cannot load dark frame [%s]", dark_name);
			e_error("skipped dark subtraction");
		} else {
			e_comment(1, "subtracting dark");
			cube_sub_im(in, dark_frame);
			image_del(dark_frame);
			e_comment(1, "switching to proportional fit");
			prop_fit = 1 ;
		}
	} else {
		e_comment(1, "no dark frame specified: skipped subtraction");
	}
	/* Fit slopes, get results */

	e_comment(0, "---> Part 3 of 4: fitting slopes");
	if (prop_fit) results = cube_create_gainmap_proportional(in);
	else results = cube_create_gainmap_robust(in) ;
	cube_del(in) ;

	if ((results == NULL) || (results[0] == NULL) || (results[1] == NULL)) {
		e_error("creating twilight flat-field: aborting") ;
		return -1 ;
	}

	e_comment(0, "---> Part 4 of 4: saving output");
	return conica_twilight_save(name_in,
                                name_out,
                                nplanes,
                                results,
                                output_intercepts,
                                output_errmap,
                                output_pixmap,
                                prop_fit,
                                lo_thresh,
                                hi_thresh) ;
}

static int conica_twilight_save(
    char		*	name_in,
    char		*	name_out,
    int				nplanes,
    image_t		**	results,
    int				output_intercepts,
    int				output_errmap,
    int				output_pixmap,
    int				prop_fit,
    double			lo_thresh,
    double			hi_thresh)
{
	qfits_header*	fh ;
	char		*	refname;
	pixelmap	*	badpixmap ;
	char			full_name[FILENAMESZ];
	image_t		*	promoted ;
	framelist	*	flist ;

	/* Save flat-field */
	e_comment(1, "saving flat-field");
	sprintf(full_name, "%s_flat.fits", name_out);
	
	/* Get name of first frame if ASCII list, or name of the cube if 3d data */
	if (is_ascii_list(name_in) == 1) {
        refname = framelist_firstname(name_in) ;
    } else {
        refname = name_in ;
    }
	e_comment(1, "using header from frame [%s]", refname);
	 
	/* Get FITS header from reference file */
	if ((fh = qfits_header_read(refname)) == NULL) {
		e_error("getting header from reference frame");
		return -1 ;
	}
	/* Prepare the header */
	conica_header_for_image(fh) ;
	
	/* Add DataFlow keywords */
	flist = framelist_load(name_in);
	conica_pro_fits(fh,
					full_name,
					"REDUCED",
					NULL,
					/* conica_imag_flat_result, */
					"Ok",
					"cal_twflats",
					nplanes,
					flist,
					NULL);
	
	/* Save list of input files as HISTORY in the header */
	qfits_header_add(fh, "COMMENT", "list of input files", NULL, NULL);
	if (is_ascii_list(name_in)==1) {
		conica_add_files_history(fh, flist) ;
	} else {
		qfits_header_add(fh, "COMMENT", get_basename(name_in), NULL, NULL);
	}
	e_comment(0, "saving twilight flat [%s]", full_name);
	image_save_fits_hdrdump(results[0], full_name, fh, BPP_DEFAULT);
	qfits_header_destroy(fh);

	/* Create and save badpixel map if requested */
	if (output_pixmap) {
		e_comment(1, "creating bad pixel map");
		badpixmap = image_threshold2pixelmap(results[0],
											 lo_thresh,
											 hi_thresh);
		if (badpixmap == NULL) {
			e_error("creating bad pixel map");
		} else {
			sprintf(full_name, "%s_badpix.fits", name_out);
			if ((promoted = pixelmap_2_image(badpixmap)) == NULL) {
				e_error("cannot promote pixelmap") ;
			} else {
				pixelmap_del(badpixmap);
				fh = qfits_header_read(refname);
				conica_header_for_image(fh) ;
				conica_pro_fits(fh,
						full_name,
						"REDUCED",
						NULL,
						/* conica_imag_flat_badpix, */
						"OK",
						"cal_twflats",
						nplanes,
						flist,
						NULL);
				image_save_fits_hdrdump(promoted,full_name,fh,BPP_8_UNSIGNED);
				image_del(promoted) ;
				qfits_header_destroy(fh) ;
				e_comment(1, "saved as [%s]", full_name);
			}
		}
	}
	
	/*
	 * Output results: different cases if the fit was linear or purely
	 * proportional
	 */
		
	/*
	 * Linear fit: results[1] has intercept map
	 *             results[2] has error map
	 */

	if (!prop_fit) {
		/* Save intercepts map if requested */
		if (output_intercepts) {
			e_comment(1, "saving intercept map");
			sprintf(full_name, "%s_intercept.fits", name_out);
			if (results[1]!=NULL) {
				fh = qfits_header_read(refname);
				conica_header_for_image(fh) ;
				conica_pro_fits(fh,
						full_name,
						"REDUCED",
						NULL,
						/* conica_imag_flat_interce, */
						"OK",
						"cal_twflats",
						nplanes,
						flist,
						NULL);
				image_save_fits_hdrdump(results[1], full_name, fh, BPP_DEFAULT) ;
				qfits_header_destroy(fh) ;
				e_comment(1, "saved as [%s]", full_name);
			} else {
				e_error("null intercept map: cannot save");
			}
		}
		/* Save error map if requested */
		if (output_errmap) {
			e_comment(1, "saving error map");
			sprintf(full_name, "%s_errmap.fits", name_out);
			if (results[2]!=NULL) {
				fh = qfits_header_read(refname);
				conica_header_for_image(fh) ;
				conica_pro_fits(fh,
						full_name,
						"REDUCED",
						NULL,
						/* conica_imag_flat_errmap, */
						"OK",
						"cal_twflats",
						nplanes,
						flist,
						NULL);
				image_save_fits_hdrdump(results[2], full_name, fh, BPP_DEFAULT) ;
				qfits_header_destroy(fh) ;
				e_comment(1, "saved as [%s]", full_name);
			} else {
				e_error("null error map: cannot save");
			}
		}
	} else {

	/*
	 * Proportional fit: results[1] has error map
	 *					 no intercept map
	 */
		if (output_intercepts) {
			e_warning("no intercept map for proportional fit");
		}

		if (output_errmap) {
			e_comment(1, "saving error map");
			sprintf(full_name, "%s_errmap.fits", name_out);
			if (results[1]!=NULL) {
				fh = qfits_header_read(refname);
				conica_header_for_image(fh) ;
                conica_pro_fits(fh,
                        full_name,
                        "REDUCED",
                        NULL,
                        /* conica_imag_flat_errmap, */
                        "OK",
                        "cal_twflats",
                        nplanes,
						flist,
						NULL);
                image_save_fits_hdrdump(results[1], full_name, fh, BPP_DEFAULT) ;
                qfits_header_destroy(fh) ;
                e_comment(1, "saved as [%s]", full_name);
			} else {
				e_error("null error map: cannot save");
			}
		}
	}

	framelist_del(flist);
	if (results[0]!=NULL) image_del(results[0]);
	if (results[1]!=NULL) image_del(results[1]);
	if (!prop_fit)
		if (results[2]!=NULL)
			image_del(results[2]);
	free(results) ;
	return 0 ;
}