/*----------------------------------------------------------------------------
   
   File name    :   sp_combine.c
   Author       :	Y. Jung
   Created on   :	Feb. 2001
   Description  :	ISAAC frame combination 

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

/*
	$Id: sp_combine.c,v 1.19 2001/10/22 12:23:34 ndevilla Exp $
	$Author: ndevilla $
	$Date: 2001/10/22 12:23:34 $
	$Revision: 1.19 $
 */

/*----------------------------------------------------------------------------
                                Includes
 ---------------------------------------------------------------------------*/
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <unistd.h>

#include "eclipse.h"
#include "isaacp_lib.h"

/*---------------------------------------------------------------------------
                            Function prototypes
 ---------------------------------------------------------------------------*/

static int sp_combine_engine(char *, char *, char *, int, int, char);
static framelist ** sp_combine_classify(char *, float *);

/* Static only because of the generic sort routine */
static char * key_offsety = NULL ;

/*----------------------------------------------------------------------------
                            Main code
 ---------------------------------------------------------------------------*/
int isaac_sp_combine_main(void * dict)
{
	dictionary	*	d ;
		
	char		*	flatname ;
	int				aver_flag ;
	int				diff_flag ;
	char			direction ;

	char			argname[10] ;
	char    	*	name_i ;
    char    	*	name_o ;
    int     		nfiles ;
	
	int				errors ;
	int				i ;
	 
	d = (dictionary*)dict ;
   	/* Get options */
	flatname  = dictionary_get(d, "arg.flatfield") ;
	key_offsety = dictionary_get(d, "arg.keyword") ;
	aver_flag = dictionary_getint(d, "arg.average", 0) ;
	diff_flag = dictionary_getint(d, "arg.difference", 0) ;
	direction = dictionary_getchar(d, "arg.direction", 'v') ;
	
	/* 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 += sp_combine_engine(name_i,
                            		name_o,
                           			flatname,
                            		aver_flag,
                            		diff_flag,
                            		direction) ;
		free(name_o) ;
	}
	return errors ;
}


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


/*-------------------------------------------------------------------------*/
/**
  @name		sp_combine_engine
  @memo		combines a series of fits images
  @param	inname			input ascii file name		
  @param	outname			output file base name
  @param	flatfield		flatfield name
  @param	aver_flag		average flag
  @param	diff_flag		difference flag
  @param	direction		direction of shifting (vertical by default)
  @return	int
 */
/*--------------------------------------------------------------------------*/
static int sp_combine_engine(
        char *      inname,
        char *      outname,
        char *      flatfield,
        int			aver_flag,
        int			diff_flag,
        char        direction)
{
    int             nfiles ;
    float           offset ;
    int             offset_int ;
    char            avgname[FILENAMESZ] ;
    char            out_diff[FILENAMESZ] ;
    char            out_comb[FILENAMESZ] ;
    image_t 	*   aver[2] ;
    image_t 	*   diff[2] ;
    image_t 	*   flat_image ;
    image_t 	*   shifted ;
    qfits_header*   hdr ;

	framelist	**	batch ;
	cube_t		*	batch_c[2] ;

	char		*	name1 ;
	char		*	name2 ;

    /* Classify the frames according to a criterion */
	e_comment(0, "classifying input frames...");
	batch = sp_combine_classify(inname, &offset);
	if (batch==NULL) {
        e_error("unable to classify the input frames") ;
        return -1 ;
    }
	
	/* Save total number of files */
	nfiles = batch[0]->n + batch[1]->n ;
	/* Save first name of first image in each batch */
	name1 = strdup(batch[0]->name[0]);
	name2 = strdup(batch[1]->name[0]);

	/* Load corresponding cube and average it immediately */
	e_comment(0, "loading frames and averaging them...");
	batch_c[0] = cube_load_strings(batch[0]->name, batch[0]->n);
	aver[0] = cube_avg_linear(batch_c[0]);
	cube_del(batch_c[0]);

	batch_c[1] = cube_load_strings(batch[1]->name, batch[1]->n);
	aver[1] = cube_avg_linear(batch_c[1]);
	cube_del(batch_c[1]);

	framelist_del(batch[0]);
	framelist_del(batch[1]);
	free(batch);

	if (aver_flag) {
		/* Save averaged frames */
		sprintf(avgname, "%s_aver1.fits", outname);
		e_comment(1, "saving check image %s in local directory", avgname);
		image_save_fits(aver[0], avgname, BPP_DEFAULT);

		sprintf(avgname, "%s_aver2.fits", outname);
		e_comment(1, "saving check image %s in local directory", avgname);
		image_save_fits(aver[1], avgname, BPP_DEFAULT);
	}

	e_comment(0, "computing cross-differences...");
    diff[0] = image_sub(aver[0], aver[1]) ;
    diff[1] = image_sub(aver[1], aver[0]) ;

    image_del(aver[0]) ;
    image_del(aver[1]) ;

    if (diff[0]==NULL || diff[1]==NULL) {
        e_error("subtraction failed - aborting") ;
        if (diff[0]) image_del(diff[0]) ;
        if (diff[1]) image_del(diff[1]) ;
		free(name1);
		free(name2);
        return -1 ;
    }

    /* Divide by the flatfield if available */
    if (flatfield != NULL) {
        if ((flat_image = image_load(flatfield)) == NULL) {
            e_error("in loading image [%s]: aborting", flatfield) ;
            image_del(diff[0]) ;
            image_del(diff[1]) ;
			free(name1);
			free(name2);
            return -1 ;
        }
		e_comment(0, "dividing by flat-field...");
        image_div_local(diff[0], flat_image) ;
        image_div_local(diff[1], flat_image) ;
        image_del(flat_image) ;
    }
    
	/* Output or not the difference images */
    if (diff_flag) {
        sprintf(out_diff, "%s_diff1.fits", outname) ;
		e_comment(1, "saving difference image %s", out_diff);
        image_save_fits_hdrcopy(diff[0], out_diff, name1, BPP_DEFAULT) ;
        sprintf(out_diff, "%s_diff2.fits", outname) ;
		e_comment(1, "saving difference image %s", out_diff);
        image_save_fits_hdrcopy(diff[1], out_diff, name2, BPP_DEFAULT) ;
    }

    /* Shift a frame */
    if (offset<0) offset_int = (int) (offset - 0.5) ;
    else offset_int = (int) (offset + 0.5) ;

	e_comment(0, "applying circular shift...");
    if (direction == 'v') {
        if ((shifted=image_shiftint_circular(diff[1], 0, offset_int)) == NULL) {
            e_error("cannot shift the image - aborting") ;
            image_del(diff[0]) ;
            image_del(diff[1]) ;
            return -1 ;
        }
    } else if (direction == 'h') {
        if ((shifted=image_shiftint_circular(diff[1], offset_int, 0)) == NULL) {
            e_error("cannot shift the image - aborting") ;
            image_del(diff[0]) ;
            image_del(diff[1]) ;
			free(name1);
			free(name2);
            return -1 ;
        }
    } else {
        e_error("unsupported shift direction - aborting") ;
        image_del(diff[0]) ;
        image_del(diff[1]) ;
		free(name1);
		free(name2);
        return -1 ;
    }
    image_del(diff[1]) ;

    /* Combine the images */
	e_comment(0, "recombining images...");
    if (image_add_local(diff[0], shifted)) {
        e_error("image addition failed - aborting") ;
        image_del(diff[0]) ;
        image_del(shifted) ;
		free(name1);
		free(name2);
        return -1 ;
    }
    image_del(shifted) ;

    /* Set output file name */
    sprintf(out_comb, "%s_comb.fits", outname) ;
	e_comment(0, "saving output file %s", out_comb);

    /* Read the FITS header of first image in first batch */
    if ((hdr = qfits_header_read(name1)) == NULL) {
        e_error("reading header from %s: aborting", name1);
		free(name1);
		free(name2);
        return -1 ;
    }
	free(name1);
	free(name2);

	/* Write the PRO keywords in the fits header */
    if (isaac_pro_fits(	hdr,
						out_comb,
						"REDUCED",
						NULL,
						isaac_spec_combine,
						"OK",
						"isaacp_sp_combine",
						nfiles,
						NULL,
						NULL)) {
		e_error("unable to write PRO keywords in the FITS header - aborting");
        image_del(diff[0]) ;
        qfits_header_destroy(hdr) ;
        return -1 ;
    }

    /* Output the combined image */
    image_save_fits_hdrdump(diff[0], out_comb, hdr, BPP_DEFAULT) ;
    image_del(diff[0]) ;
    qfits_header_destroy(hdr) ;

    e_comment(0, "in:  %s\n", inname) ;
    e_comment(0, "out: %s\n", outname) ;

    return 0 ;
}


/* Static only meant to be used with framelist_labelize */
static int sp_combine_offset_compare(char * file1, char * file2)
{
	char	*	cumoff_y1,
			*	cumoff_y2 ;
	int			ret ;

	/* Get offset value from both files */
	if (key_offsety==NULL) {
		cumoff_y1 = strdup(isaac_get_cumoffsety(file1));
		cumoff_y2 = strdup(isaac_get_cumoffsety(file2));
	} else {
		cumoff_y1 = strdup(qfits_query_hdr(file1, key_offsety));
		cumoff_y2 = strdup(qfits_query_hdr(file2, key_offsety));
	}
	if (!strcmp(cumoff_y1, cumoff_y2)) {
		ret = 1 ;
	} else {
		ret = 0 ;
	}
	free(cumoff_y1);
	free(cumoff_y2);
	return ret ;
}

/*-------------------------------------------------------------------------*/
/**
  @name		sp_combine_classify
  @memo		creates two ascii files with the separated files
  @param	inname			ascii file name
  @param	keyword			keyword name
  @param	offset_diff		differences between offsets array
  @param	outfile1		first outfile name
  @param	outfile2		second outfile name
  @return	total number of encountered valid frames
 */
/*--------------------------------------------------------------------------*/
static framelist ** sp_combine_classify(
        char 	*	inname,
        float 	*   offset_diff)
{
	framelist 	*	flist ;
	framelist 	**	batch ;
	int				nlabels ;
	char		*	first_offset,
				*	second_offset,
				*	keyval ;
	int				i ;

    /* Load all the file names contained in the list */
    if ((flist = framelist_load(inname)) == NULL) {
        e_error("reading ASCII list [%s]: aborting load", inname);
        return NULL ;
    }
	if (flist->n<2) {
		e_error("found %d files in framelist instead of minimum 2", flist->n);
		framelist_del(flist);
		return NULL ;
	}

	nlabels = framelist_labelize(flist, sp_combine_offset_compare);
    if (nlabels != 2) {
        e_error("found %d different values for Y offsets, should be 2",
                nlabels) ;
		framelist_del(flist);
        return NULL ;
    }
    /* Find first_offset, second_offset and offset_diff */
	if (key_offsety==NULL) {
		keyval = isaac_get_cumoffsety(flist->name[0]) ;
	} else {
		keyval = qfits_query_hdr(flist->name[0], key_offsety);
	}
    if (!qfits_is_float(keyval)) {
        e_error("reading Y offset for file %s", flist->name[0]);
		framelist_del(flist);
        return NULL ;
    }
    first_offset = strdup(keyval) ;

    for(i=1 ; i<flist->n ; i++) {
		if (key_offsety==NULL) {
			keyval = isaac_get_cumoffsety(flist->name[i]);
		} else {
			keyval = qfits_query_hdr(flist->name[i], key_offsety);
		}
        if (!qfits_is_float(keyval)) {
			e_error("reading Y offset for file %s", flist->name[i]);
			framelist_del(flist);
			free(first_offset) ;
            return NULL ;
        }
		if (strcmp(keyval, first_offset)) {
            second_offset = strdup(keyval) ;
            break ;
        }
    }
    *offset_diff = atof(second_offset) - atof(first_offset) ;
	free(first_offset);
	free(second_offset);

    /* Classify according to the keyword values  */

	batch = malloc(2 * sizeof(framelist*));
	batch[0] = framelist_select(flist, 0);
	batch[1] = framelist_select(flist, 1);
	framelist_del(flist);
    return batch ;
}