/*----------------------------------------------------------------------------
   
   File name    :   utils.c
   Author       :	N. Devillard
   Created on   :	August 1998
   Description  :	ISAAC various utilities 

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

/*
 $Id: utils.c,v 1.5 2001/10/22 12:23:31 ndevilla Exp $
 $Author: ndevilla $
 $Date: 2001/10/22 12:23:31 $
 $Revision: 1.5 $
 */

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

#include "utils.h"
#include "classif.h"

/*----------------------------------------------------------------------------
						Function ANSI C code 
 ---------------------------------------------------------------------------*/

/*-------------------------------------------------------------------------*/
/**
  @name		isaac_check_instrument
  @memo		Check that a FITS file or ASCII list is an ISAAC product.
  @param	name		Name of the frame or list to check
  @return	1 if ISAAC product, 0 if not, -1 if error occurred.
  @doc

  This function checks that a FITS frame or an ASCII list are indeed ISAAC
  products (i.e. INSTRUME='ISAAC'). If it is not the case, a warning is
  produced by this function to signal this, and the function returns 0. If
  all frames check Ok, the function is silent and returns 1. If an error
  occurs (non-existing frame or list), an error message is produced and -1
  is returned.
 */
/*--------------------------------------------------------------------------*/

int isaac_check_instrument(char * name)
{
	char	  *	instrum ;
	framelist *	flist ;
	int			i ;
	int			status ;

	/* Set status to valid */
	status = 1 ;

	/* Is the input name the name of a FITS frame? */
	if (is_fits_file(name)==1) {
		/* Get INSTRUME keyword from header */
		instrum = qfits_query_hdr(name, "INSTRUME");
		if (instrum==NULL) {
			/* No INSTRUME keyword found: not ISAAC */
			e_warning("file [%s] is not an ISAAC product", name);
			status = 0 ;
		} else {
			/* Check value of INSTRUME keyword */
			instrum = qfits_pretty_string(instrum);
			if (strncmp(instrum, "ISAAC", 5)) {
				/* Not 'ISAAC' */
				e_warning("file [%s] is not an ISAAC product", name);
				status = 0 ;
			}
		}
	} else if (is_ascii_list(name)==1) {
		/* The input is an ASCII list: load the names */
		flist = framelist_load(name);
		if (flist==NULL) {
			/* Error loading list */
			e_error("loading ASCII list [%s]", name);
			status = -1 ;
		} else {
			/* Browse through FITS files */
			for (i=0 ; i<flist->n ; i++) {
				/* Get INSTRUME keyword */
				instrum = qfits_query_hdr(flist->name[i], "INSTRUME");
				if (instrum==NULL) {
					/* INSTRUME not found: not ISAAC */
					e_warning("file [%s] is not an ISAAC product",
								flist->name[i]);
					status = 0 ;
				} else {
					/* Check value of INSTRUME keyword */
					instrum = qfits_pretty_string(instrum);
					if (strncmp(instrum, "ISAAC", 5)) {
						/* Not 'ISAAC' */
						e_warning("file [%s] is not an ISAAC product",
									flist->name[i]);
						status = 0 ;
					}
				}
			}
			framelist_del(flist);
		}
	} else {
		/* Not a FITS frame nor an ASCII list */
		e_error("[%s] is not a FITS frame nor an ASCII list", name);
		status = -1 ;
	}
	return status ;
}


/*-------------------------------------------------------------------------*/
/**
  @name		isaac_ff_dark_badpix_handling
  @memo		apply dark subtraction and ff division, replace bad pixels
  @param	in	pointer to allocated cube
  @param	ff_name	flat field name
  @param	dark_name	dark name
  @return	badpix_name	bad pixel name
  @doc
	The input cube is always modified, to ensure that the
	returned value is read/write.
 */
/*--------------------------------------------------------------------------*/
void isaac_ff_dark_badpix_handling(
		cube_t	**	in,
		char	*	ff_name,
		char	*	dark_name,
		char	*	badpix_name)
{
	image_t		*	dark ;
	image_t		*	ff ;
	pixelmap	*	badpix ;
	long			i, p ;
	pixelvalue	*	pix ;
	int				has_dark,
					has_ff ;
	int				nbpix ;

	/* Initialize */
	dark = NULL ;
	ff   = NULL ;

	/*
	 * see if dark was given and can be found
	 * note: further tests should be done to ensure consistency, too
	 */

	has_dark=0 ;
	if (dark_name && *dark_name) {
		dark = image_load(dark_name) ;
		if (dark==NULL) {
			e_error("cannot load dark file [%s]", dark_name) ;
		} else if ((dark->lx != (*in)->lx) || (dark->ly != (*in)->ly)) {
			e_error("incompatible sizes for dark and jitter cube") ;
			e_error("dark image size is [%d x %d]", dark->lx, dark->ly);
			image_del(dark) ;
		} else {
			has_dark = 1 ;
		}
	}

	/*
	 * see if flat-field was given and can be found
	 * note: further tests shouls be done to ensure consistency, too
	 */

	has_ff=0 ;
	if (ff_name && *ff_name) {
		ff = image_load(ff_name) ;
		if (ff==NULL) {
			e_error("cannot load flat-field file [%s]", ff_name) ;
		} else if ((ff->lx != (*in)->lx) || (ff->ly != (*in)->ly)) {
			e_error("incompatible sizes for flat-field and jitter cube") ;
			e_error("flat-field image size is [%d x %d]", ff->lx, ff->ly);
			image_del(ff) ;
		} else {
			has_ff = 1 ;
		}
	}

	/*
	 * If no input was provided, return the input cube untouched
	 */
	if (!has_dark && !has_ff) {
		e_comment(1, "flat-field division and dark subtraction skipped");
		return ;
	}

	/*
	 * Different cases:
	 * dark provided
	 * flat-field provided
	 * dark & flat-field provided
	 */

	/* only dark was provided */
	if (has_dark && !has_ff) {
		e_comment(1, "applying dark subtraction") ;
		cube_sub_im((*in), dark) ;
		image_del(dark) ;
		e_comment(1, "no flat-field provided: skipping") ;
	}
	
	/* only flat-field was provided */
	if (!has_dark && has_ff) {
		e_comment(1, "no dark provided: skipped") ;
		e_comment(1, "applying flat-field division") ;
		cube_div_im((*in), ff) ;
		image_del(ff) ;
	}

	/* both dark and flat-field have been provided */
	if (has_dark && has_ff) {
		e_comment(1, "applying dark subtraction and flat-field division");
		for (p=0 ; p<(*in)->np ; p++) {
			pix = (*in)->plane[p]->data ;
			nbpix = (*in)->plane[p]->lx * (*in)->plane[p]->ly ;
			for (i=0 ; i<nbpix ; i++) {
				*pix = (*pix - dark->data[i]) / ff->data[i] ;
				pix++ ;
			}
		}
		image_del(dark) ;
		image_del(ff) ;
	}

	/*
	 * Apply now bad pixel correction if needed
	 */

	if (badpix_name && *badpix_name) {
		badpix = pixelmap_load(badpix_name) ;
		if (badpix == NULL) {
			e_error("cannot load bad pixel map [%s]: skipping",
					badpix_name) ;
		} else {
			e_comment(1, "applying dead pixel correction") ;
			cube_clean_deadpix(*in, badpix) ;
			pixelmap_del(badpix) ;
		}
    } else {
		e_comment(1, "bad pixel replacement: skipped") ;
	}
	return ;
}


/*-------------------------------------------------------------------------*/
/**
  @name     spectro_jitter_average
  @memo     shift(?) and average images
  @param    cube    cubes array 
  @param    nbcubes	nb of cubes
  @param	offsets_list	offsets list
  @param	t_switch		array that associates images to cubes
  @param	noffsets	nb entries in t_switch and offsets_list
  @param	average_status	status flag. MODIFIED.
  @return   averaged images (1 image per input cube)
 */
/*--------------------------------------------------------------------------*/
image_t ** spectro_jitter_average(
        cube_t         	**  cube,
		int					nbcubes,
		double			*	offsets_list,
		int				*	t_switch,
		int					noffsets,
		int				*	average_status)
{
    cube_t         **  shifted ;
    image_t        **  averaged ;
    double              shift_val ;
    double              off_val1 ;
    double              off_val2 ;
    int                 flag_shifted ;

    int                 i, j, k ;

    flag_shifted = 0 ;

    /* Allocate shifted */
    shifted = malloc(nbcubes*sizeof(cube_t*)) ;

    /* For each cube, shift the frames according to the first frame */
    for (i=0 ; i<nbcubes ; i++) {
        compute_status("Shift frames in each cube (if necessary)", i, nbcubes, 1) ;
        /* Create the shifted cubes */
        shifted[i] = cube_new(cube[i]->lx, cube[i]->ly, cube[i]->np) ;
        if (shifted[i] == NULL) {
            e_error("cannot create new cube") ;
            for (j=0 ; j<i ; j++) {
                cube_del(shifted[j]) ;
            }
            free(shifted) ;
            return NULL ;
        }
        for (j=0 ; j<cube[i]->np ; j++) {
            off_val1 = get_offset_in_list(offsets_list,
										noffsets,
										t_switch,
										i,
										0) ;
            off_val2 = get_offset_in_list(offsets_list,
										noffsets,
										t_switch,
										i,
										j) ;
            shift_val = off_val2 - off_val1 ;
            if (shift_val != 0.0) {
                flag_shifted = 1 ;
                /* Shift the current image   */
				shifted[i]->plane[j] = image_shift(cube[i]->plane[j],
												0,
												shift_val,
												(double*)NULL) ;
                if (shifted[i]->plane[j] == NULL) {
                    e_error("while shifting the frames") ;
                    for (k=0 ; k<i ; k++) {
                        cube_del(shifted[k]) ;
                    }
                    free(shifted) ;
                    return NULL ;
                }
            } else {
				shifted[i]->plane[j] = image_copy(cube[i]->plane[j]) ;
            }
        }
    }

    if (flag_shifted == 0) {
        e_comment(2, "No shift computed") ;
    }

    /* Allocate memory for averaged images */
    averaged = malloc(nbcubes*sizeof(image_t*)) ;

    /* Average the cubes   */
    for (i=0 ; i<nbcubes ; i++) {
        if ((averaged[i]=cube_avg_linear(shifted[i])) == NULL) {
            *average_status = ALGO_FAILED ;
            for (j=0 ; j<i ; j++) {
                image_del(averaged[j]) ;
            }
            free(averaged) ;
            return NULL ;
        }
	}
    *average_status = ALGO_OK ;

    for (i=0 ; i<nbcubes ; i++) {
        cube_del(shifted[i]) ;
    }
    free(shifted) ;
    return averaged ;
}


/*-------------------------------------------------------------------------*/
/**
  @name		spectro_jitter_extract
  @memo		Specrum detection and extraction for jitter-like recipes 
  @param	combined	combined image in which the spec is extracted
  @param	spec_width	Width of the spectrum
  @param	lo_dist	Dist. spec -> lower resid. sky
  @param	hi_dist	Dist. spec -> high resid. sky
  @param	lo_width	Width of the low residual sky
  @param	hi_width	Width of the high residual sky
  @param	apply_filter	Flag to apply a median filter before extraction
  @param	sky_frame	Frame for sky extraction (not done if NULL provided)
  @param	main_offset_diff	Offsets between black and white lines
  @param	nb_offsets	Nb of offsets
  @param	status_wavecal_done	Flag to specify that the disprel is available
  @param	wavecal	Dispersion relation
  @param	nb_coeffs	Nb of coeffs of the polynomial disp relation
  @param	spec_position	Spectrum position. Found if <0 provided. MODIFIED
  @param	spec_detected	Flag. MODIFIED
  @param	spec_extracted	Flag. MODIFIED
  @param	status	Status of algorithm. MODIFIED
  @return	A double3 array: x=xcoordinate, y=extracted spec, z=sky spec
  @doc
  The spectra are supposed to be horizontal
 */
/*--------------------------------------------------------------------------*/
double3 * spectro_jitter_extract(
        image_t        *   combined,
		int					spec_width,
		int					lo_dist,
		int					hi_dist,
		int					lo_width,
		int					hi_width,
		int					apply_filter,
		char			*	sky_frame,
		double			*	main_offset_diff,
		int					nb_offsets,
		int					status_wavecal_done,
		double			*	wavecal,
		int					nb_coeffs,
		int				*	spec_position,
		int				*	spec_detected,
		int				*	spec_extracted,
		int				*	status)	
{
	double3			*	out_array ;
    int                 slit_length,
						spec_length ;
    double3         *   position ;
    int                 offset_diff ;
    int                 low_side,
                        up_side ;
    int                 sky_pos[4] ;
    image_t        	*   filtered ;
    image_t        	*   skyframe ;
    pixelvalue          res_sky_estim ;
    pixelvalue          sky_signal ;
    double              median_1 ;
    double              median_2 ;
    image_t        *   extr_line ;

    int                 i ;

	/* Initialize variables */
	slit_length = combined->ly ;
	spec_length = combined->lx ;
	skyframe = NULL ;
	sky_signal = 0.0 ;

	/* Try to detect the brightest spectrum if its position is unknown */
    if (*spec_position < 0) {
        i = 0 ;
        /* Try to find 2 black shadows */
        position = NULL ;
        while ((position == NULL) && (i < nb_offsets)) {
            offset_diff = (int)fabs(main_offset_diff[i]) ;
            if ((position=find_brightest_spectrum_1d(combined,
											offset_diff,
											EQUALLY_SPACED_SHADOW_SPECTRA,
											0)) == NULL) {
                if (i == nb_offsets - 1) {
                    e_warning("Detection failed - try with lower criteria") ;
                } else {
                    e_warning("Detection failed - try with next offset") ;
                }
                i++ ;
            }
        }

        /* Try with 1 black shadow   */
        if (position == NULL) {
            i = 0 ;
            while ((position == NULL) && (i < nb_offsets)) {
                offset_diff = (int)fabs(main_offset_diff[i]) ;
                if ((position=find_brightest_spectrum_1d(combined,
                                                    offset_diff,
                                                    ONE_SHADOW_SPECTRUM,
                                                    0)) == NULL) {
                    if (i < nb_offsets - 1) {
                        e_warning("Detection failed - try with next offset") ;
                    }
                    i++ ;
                }
            }
        }

		/* Exit if detection fails */
        if (position == NULL) {
            *spec_detected  = 0 ;
            *spec_extracted = 0 ;
            *status         = ALGO_SKIPPED ;
            return NULL ;
        } else *spec_detected = 1 ;

		/* Set the spectrum position */
        *spec_position = (int)(position->y[0]) ;
        double3_del(position) ;
        
	/* Else use the position specified in the INI file */
    } else *spec_detected = 0 ;

    /* Set the parameters for the extraction */

    /* Spectrum position */
    low_side = (int)(*spec_position - (spec_width/2)) ;
    up_side = low_side + spec_width ;
    if ((low_side < 1) || (up_side > slit_length)) {
        e_error("spectrum position out of the image - aborting") ;
        *spec_extracted = 0 ;
        *status = ALGO_SKIPPED ;
        return NULL ;
    }

    /* Residual Sky position */
    if (lo_dist < 0) lo_dist = 2*spec_width ;
    if (hi_dist < 0) hi_dist = 2*spec_width ;
    if (lo_width < 0) lo_width = 10 ;
    if (hi_width < 0) hi_width = 10 ;
    sky_pos[1] = (int)(*spec_position - lo_dist) ;
    sky_pos[0] = (int)(sky_pos[1] - lo_width) ;
    sky_pos[2] = (int)(*spec_position + hi_dist) ;
    sky_pos[3] = (int)(sky_pos[2] + hi_width) ;

    /* Allocate extracted array */
	out_array = double3_new(spec_length) ;

    if (apply_filter == 1) {
        /* Filter the combined image */
        if ((filtered=image_filter_median(combined)) == NULL) {
            e_warning("cannot filter the combined image") ;
            filtered = image_copy(combined) ;
        }
    } else filtered = image_copy(combined) ;

    /* Load the sky frame (if provided) to extract the sky signal */
    if (sky_frame != NULL) skyframe = image_load(sky_frame) ;

    /* Extract the spectrum and get rid of the residual sky */
    for (i=0 ; i<spec_length ; i++) {
		/* Estimate the SKY */
		if (((sky_pos[0] < 1) || (lo_width == 0)) &&
				((sky_pos[3] <= slit_length) && (hi_width > 0))) {
			res_sky_estim = image_getmedian_vig(filtered,
					i+1,
					sky_pos[2],
					i+1,
					sky_pos[3]) ;
			if (sky_frame != NULL) {
				sky_signal = image_getmedian_vig(skyframe,
						i+1,
						sky_pos[2],
						i+1,
						sky_pos[3]) ;
			}
		} else if (((sky_pos[3] > slit_length) || (hi_width == 0)) 
				&& ((sky_pos[0] > 0) && (lo_width > 0))) {
			res_sky_estim = image_getmedian_vig(filtered,
					i+1,
					sky_pos[0],
					i+1,
					sky_pos[1]) ;
			if (sky_frame != NULL) {
				sky_signal = image_getmedian_vig(skyframe,
						i+1,
						sky_pos[0],
						i+1,
						sky_pos[1]) ;
			}
		} else if ((lo_width != 0) && (hi_width != 0)
				&& (sky_pos[0] > 0) && (sky_pos[3] <= slit_length)) {
			median_1 = image_getmedian_vig(filtered,
					i+1,
					sky_pos[0],
					i+1,
					sky_pos[1]) ;
			median_2 = image_getmedian_vig(filtered,
					i+1,
					sky_pos[2],
					i+1,
					sky_pos[3]) ;
			res_sky_estim=(median_1 + median_2)/2 ;
			if (sky_frame != NULL) {
				median_1 = image_getmedian_vig(skyframe,
						i+1,
						sky_pos[0],
						i+1,
						sky_pos[1]) ;
				median_2 = image_getmedian_vig(skyframe,
						i+1,
						sky_pos[2],
						i+1,
						sky_pos[3]) ;
				sky_signal=(median_1 + median_2)/2 ;
			}
		} else {
			res_sky_estim = 0 ;
			sky_signal = 0 ;
		}

		/* Estimate the spectrum */
		if ((extr_line = image_getvig(filtered,
						i+1,
						low_side,
						i+1,
						up_side)) == NULL) {
			e_error("error in line extraction - aborting") ;
			*spec_extracted = 0 ;
			*status = ALGO_SKIPPED ;
			if (sky_frame != NULL) image_del(skyframe) ;
			image_del(filtered) ;
			double3_del(out_array) ;
			return NULL ;
		}

		/* Write in output double3 */
        out_array->y[i] = (double)image_getsumpix(extr_line) ;
        image_del(extr_line);
		out_array->y[i] -= spec_width * res_sky_estim ;
        if (status_wavecal_done == ALGO_OK) {
			out_array->x[i]=(double)(wavecal[0]+wavecal[1]*(i+1)) ;
        } else out_array->x[i]=(double)(i+1) ;
        if (sky_frame != NULL) out_array->z[i] = (double)sky_signal ;
    }
    *spec_extracted = 1 ;
    *status = ALGO_OK ;
	
	/* Free and return */
    image_del(filtered) ;
    if (sky_frame != NULL) image_del(skyframe) ;
    return out_array ;
}


/*-------------------------------------------------------------------------*/
/**
  @name     divide_by_fit
  @memo     Fit the input image and divide it by the fit
  @param    in          input image
  @param    order       order of the fit
  @param    xsize       central X size to be used for the fit
  @param    offset      offset
  @param    setting_nb  setting id
  @param    pair_nb     pair id
  @param    output_poly_image   flag to output the polynomial images
  @return   corrected image
 */ 
/*--------------------------------------------------------------------------*/
image_t * divide_by_fit(
        image_t    *   in,
        int             order,
        int             xsize,
        int             offset,
        int             setting_nb,
        int             pair_nb,
        int             output_poly_image)
{
    image_t    *   out ;
    int             xstart,
                    xend,
                    ystart,
                    yend ;
    image_t    *   extracted ;
    image_t    *   collapsed ;
    image_t    *   extracted2 ;
    image_t    *   fit_image ;

    int             nb_samples ;
    double3     *   to_fit ;
    double      *   coeffs ;
    char            poly_string[FILENAMESZ] ;
    char            poly_string2[FILENAMESZ] ;
    char            poly_name[FILENAMESZ] ;

    int             i ;

    /* Determine the zone to extract */
    xstart = (int)((in->lx - xsize)/2) + 1 ;
    xend = xstart + xsize - 1 ;
    if ((xstart<1) || (xend>in->lx)) {
        e_error("bad X size specified") ;
        return NULL ;
    }

    /* Extract the central zone */
    if ((extracted=image_getvig(in,
                                    xstart,
                                    1,
                                    xend,
                                    in->ly)) == NULL) {
        e_error("cannot extract image") ;
        return NULL ;
    }

    /* Collapse the extracted zone */
    if ((collapsed=image_collapse(extracted, 1)) == NULL) {
        e_error("cannot collapse the image") ;
        image_del(extracted) ;
        return NULL ;
    }
    image_del(extracted) ;

    /* Find the 'valid' zone in the 1D image */
    ystart = 1 ;
    while ((collapsed->data[ystart-1] == (pixelvalue)0) &&
            (ystart < in->lx)) ystart++ ;
    ystart += offset ;

    yend = collapsed->ly ;
    while ((collapsed->data[yend-1] == (pixelvalue)0) &&
            (yend > 1)) yend-- ;
    yend -= offset ;

    if (ystart > yend) {
        e_error("invalid coordinates of the zone to extract") ;
        image_del(collapsed) ;
        return NULL ;
    }

    /* Extract the 1D signal to fit */
    if ((extracted2=image_getvig(collapsed,
                                            1,
                                            ystart,
                                            1,
                                            yend)) == NULL) {
        e_error("cannot extract 1D image") ;
        image_del(collapsed) ;
        return NULL ;
    }
    image_del(collapsed) ;

    nb_samples = extracted2->ly ;

    /* Fill to_fit */
    to_fit = double3_new(nb_samples) ;
    for (i=0 ; i<nb_samples ; i++) {
        to_fit->x[i] = (double)(ystart + i) ;
        to_fit->y[i] = (double)(extracted2->data[i] / xsize) ;
    }
    image_del(extracted2) ;

    /* Find polynomial coefficients */
    coeffs=fit_1d_poly(order-1, to_fit, NULL);
	double3_del(to_fit);
	if (coeffs==NULL) {
        e_error("cannot fit the 1D signal") ;
        return NULL ;
    }

    /* Build the fit image   */
	sprintf(poly_string, "(0,0)") ;
    for (i=1 ; i<order ; i++) {
        sprintf(poly_string2, "%s (0,%d)", poly_string, i) ;
        sprintf(poly_string, poly_string2) ;
    }
    if ((fit_image=image_gen_polynomial(in->lx,
                                in->ly,
                                coeffs,
                                order,
                                order-1,
                                poly_string)) == NULL) {
        e_error("cannot generate polynomial image") ;
        free(coeffs) ;
        return NULL ;
    }
    free(coeffs) ;

    if (output_poly_image) {
        if (pair_nb == 0) {
            sprintf(poly_name, "poly_set%d.fits", setting_nb) ;
        } else {
            sprintf(poly_name, "poly_set%d_pair%d.fits", setting_nb, pair_nb) ;
        }
        image_save_fits(fit_image, poly_name, BPP_DEFAULT) ;
    }


    /* Divide the input image by the polynomial image */
    if ((out=image_div(in, fit_image)) == NULL) {
        e_error("cannot divide the images") ;
        image_del(fit_image) ;
        return NULL ;
    }

    /* Free and return */
    image_del(fit_image) ;
    return out ;
}


/*
 * This function only converts an integer into a character string...
 */
char * algo_stat(int s)
{
    switch (s) {
        case ALGO_NOTREACHED: 	return "NotReached"; 	break ;
        case ALGO_OK: 			return "Ok"; 			break ;
        case ALGO_SKIPPED: 		return "Skipped" ; 		break ;
        case ALGO_FAILED: 		return "Failed"; 		break ;

        default: break ;
    }
	return "Unknown";
}