/*----------------------------------------------------------------------------
   
   File name    :   distortion.c
   Author       :	Y. Jung
   Created on   :	Feb 2001
   Description  :	ISAAC distortion utilities 

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

/*
   $Id: distortion.c,v 1.13 2002/03/04 13:56:08 yjung Exp $
   $Author: yjung $
   $Date: 2002/03/04 13:56:08 $
   $Revision: 1.13 $
 */

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

#include "isaacp_lib.h"
#include "distortion.h"

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

#define SQR(x) ((x)*(x))
#define         SW_FWHM_YMIN        420
#define         SW_FWHM_YMAX        460
#define         LW_FWHM_YMIN        400
#define         LW_FWHM_YMAX        600
#define         LINE_HALF_LENGTH    10

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

/*-------------------------------------------------------------------------*/
/**
  @brief    Determine the distortion, correct, and wl calibrate
  @param    in                      image object
  @param    inimage_name            input image name
  @param    listnames               input files names
  @param    nb_coeffs               number of polynomial coefficients
  @param    rejected_ends           number of pix. to reject on image borders
  @param    rej_left                nb of pix to rej. on the left (lines id.)
  @param    rej_right               on the right ...
  @param    line_table              catalog
  @param    auto_dark_subtraction   flag to get rid of dark
  @param    file_name               output file name
  @param	mode					's' or 'l' for SW or LW mode
  @param    out_corrected           flag to output corrected frames
  @param    arcs_fwhm               array with arcs pos and fwhm 
  @return   3 columns with the 2d poly and the last column with the
            1D polynomial describing the dispersion
 */
/*--------------------------------------------------------------------------*/
double ** compute_arc_reduction(
        image_t     *   in,
        char        *   inimage_name,
        framelist   *   listnames,
        int         *   nb_coeffs,
        int             rejected_ends,
        int             rej_left,
        int             rej_right,
        char        *   line_table,
        int             auto_dark_subtraction,
        char        *   file_name,
		char			mode,
        int             out_corrected,
        double3     **  arcs_fwhm)
{
    image_t         *   corrected ;
    image_t         *   collapsed ;
    poly2d          *   coeffs ;
    poly2d          *   poly_id ;
    int                 xmin,
                        ymin,
                        xmax,
                        ymax ;
    double          *   disprel ;
    double              wave_min, wave_max ;
    int                 spec_length ;
    int                 discard_lo, discard_hi ;
    char                name[FILENAMESZ] ;
    qfits_header    *   fh ;
    int                 nb_arcs ;
    double          *   arcs ;
    pixelvalue      *   line ;
    int                 line_start,
                        line_stop ;
    int                 maxpos ;
    double          **  arc_array ;
	int					fwhm_ymin,
						fwhm_ymax ;
	pro_catg_key		pro_category ;					
	char				recipe_id[FILENAMESZ] ;
	
    int                 i, j ;

    /* Initialize */
    xmin = 0 ;
    xmax = in->lx - 1 ;
    ymin = rejected_ends ;
    ymax = in->ly - 1 - rejected_ends ;
	if (mode == 's') {
		fwhm_ymin = SW_FWHM_YMIN ;
		fwhm_ymax = SW_FWHM_YMAX ;
		pro_category = isaac_spec_sw_arc_corr ;
		sprintf(recipe_id, "spec_tec_arc") ;
	} else if (mode == 'l') {
		fwhm_ymin = LW_FWHM_YMIN ;
		fwhm_ymax = LW_FWHM_YMAX ;
		pro_category = isaac_spec_lw_arc_corr ;
		sprintf(recipe_id, "lw_spec_arc") ;
	} else {
		e_error("unrecognized mode") ;
		return NULL ;
	}
	
    /* Distortion estimation */
    e_comment(1, "estimate the distortion") ;
    spec_length = in->lx ;

    if ((coeffs=isaac_compute_distortion(in,
                                        xmin,
                                        ymin,
                                        xmax,
                                        ymax,
                                        auto_dark_subtraction,
                                        &nb_arcs,
                                        &arcs)) == NULL) {
        e_error("in compute distortion") ;
        return NULL ;
    }
    *nb_coeffs = coeffs->nc ;

    /* Correction of the distorsion */
    e_comment(1, "correct the distortion of the input image") ;
    if ((poly_id=poly2d_build_from_string("0 1 1.0")) == NULL) {
        e_error("in buiding a poly2d object") ;
        poly2d_free(coeffs) ;
        free(arcs) ;
        return NULL ;
    }
    if ((corrected=image_warp_generic(in,
                                "default",
                                coeffs,
                                poly_id)) == NULL) {
        e_error("in the correction of the distorsion") ;
        poly2d_free(coeffs) ;
        poly2d_free(poly_id) ;
        free(arcs) ;
        return NULL ;
    }

    /* Find out the FWHM of the used arcs */
    (*arcs_fwhm) = double3_new(nb_arcs) ;
    if ((collapsed = image_collapse_vig(corrected,
            1, fwhm_ymin, corrected->lx, fwhm_ymax, 0)) == NULL) {
        e_error("cannot create collapsed image") ;
        double3_del(*arcs_fwhm) ;
        *arcs_fwhm = NULL ;
        free(arcs) ;
    } else {
        for (i=0 ; i<nb_arcs ; i++) {
            (*arcs_fwhm)->x[i] = arcs[i] ;
            line = malloc((2*LINE_HALF_LENGTH+1)*sizeof(pixelvalue)) ;
            line_start = arcs[i] - LINE_HALF_LENGTH ;
            line_stop = arcs[i] + LINE_HALF_LENGTH ;
            maxpos = LINE_HALF_LENGTH ;
            if (line_start < 1) {
                line_start = 1 ;
                line_stop = 2*LINE_HALF_LENGTH + 1 ;
            }
            if (line_stop > collapsed->lx) {
                line_start = collapsed->lx - 2*LINE_HALF_LENGTH - 1 ;
                line_stop = collapsed->lx ;
            }
            for (j=0 ; j<2*LINE_HALF_LENGTH+1 ; j++) {
                line[j] = collapsed->data[line_start + j] ;
            }
            (*arcs_fwhm)->y[i] = function1d_get_fwhm(line,
                                            2*LINE_HALF_LENGTH+1,
                                            &maxpos,
                                            NULL) ;
            free(line) ;
        }
        free(arcs) ;
        image_del(collapsed) ;
    }

    /* Allocate the output array */
    arc_array = malloc(4*sizeof(double*)) ;
    for (i=0 ; i<4 ; i++) {
        arc_array[i] = malloc(*nb_coeffs*sizeof(double)) ;
    }

    /* Fill the 3 first columns of the output array */
    for (i=0 ; i<*nb_coeffs ; i++) {
        arc_array[0][i] = (double)coeffs->px[i] ;
        arc_array[1][i] = (double)coeffs->py[i] ;
        arc_array[2][i] = (double)coeffs->c[i] ;
    }
    poly2d_free(coeffs) ;
    poly2d_free(poly_id) ;

    /* Output the corrected images if required */
    if (out_corrected) {
        sprintf(name, "%s_corrected.fits", get_rootname(file_name)) ;
        /* Read the FITS header of the input file    */
        fh = qfits_header_read(inimage_name) ;
        isaac_header_for_image(fh) ;
        isaac_pro_fits(fh,
                   name,
                   "REDUCED",
                   NULL,
                   pro_category,
                   "OK",
                   recipe_id,
                   1,
                   listnames,
                   NULL) ;
        /* Write the used line table in the header */
        qfits_header_add(fh, "HIERARCH ESO PRO CATALOG", line_table,
                "Catalog used", NULL) ;
        /* Write HISTORY keywords in the header */
        if (isaac_add_files_history(fh, listnames) == -1) {
            e_warning("cannot write HISTORY keywords in out file") ;
        }
        image_save_fits_hdrdump(corrected, name, fh, BPP_DEFAULT) ;
        qfits_header_destroy(fh) ;
        e_comment(0, "Arc corrected image produced: [%s]", name) ;
    }

    /* Wavelength calibration */
    e_comment(1, "Wavelength calibration on the corrected image") ;
    /* First estimation using a physical model */
    if ((disprel=isaac_get_disprel_estimate(inimage_name,
                                            spec_length,
                                            2)) == NULL) {
        e_error("cannot estimate the dispersion relation") ;
        image_del(corrected) ;
        for (i=0 ; i<4 ; i++) {
            free(arc_array[i]) ;
        }
        free(arc_array) ;
        double3_del(*arcs_fwhm) ;
        return NULL ;
    }
    wave_min = disprel[0] + disprel[1] + disprel[2] ;
    wave_max = disprel[0] +
               disprel[1] * spec_length +
               disprel[2] * spec_length * spec_length ;
    free(disprel);

    discard_lo = discard_hi = rejected_ends ;
    if ((disprel=spectro_compute_disprel_from_table(corrected,
                                        discard_lo,
                                        discard_hi,
                                        rej_left,
                                        rej_right,
                                        0,  /* remove thermal */
                                        line_table,
                                        wave_min,
                                        wave_max)) == NULL) {
        e_warning("cannot compute the dispersion relation") ;
        /* Fill the 4th column of the output array with 0*/
        for (i=0 ; i<*nb_coeffs ; i++) {
            arc_array[3][i] = (double)0 ;
        }
    } else {
        /* Fill the 4th column of the output array */
        arc_array[3][0] = (double)disprel[0] ;
        arc_array[3][1] = (double)disprel[1] ;
        for (i=2 ; i<*nb_coeffs ; i++) {
            arc_array[3][i] = (double)0 ;
        }
    }

    /* Free and return */
    image_del(corrected) ;
    if (disprel != NULL) free(disprel) ;
    return arc_array ;
}


/*-------------------------------------------------------------------------*/
/**
  @brief    Write the out fits table 
  @param    outname         outfile name    
  @param    nb_coeffs       nb of coeffs
  @param    out_table       result of engine
  @param    inimage_name    one input image name
  @param	mode			'l' or 's' for LW or SW mode
  @param    lnames          input files list
  @param    lines_type      used lines type ("argon", "xenon", argon+xenon")
  @param    arcs_fwhm       array with arcs positions and fwhm
  @return   error code: 0 ok, -1 in error case
 */
/*--------------------------------------------------------------------------*/
int arc_write_outfile(
        char        *   outname,
        int             nb_coeffs,
        double      **  out_table,
        char        *   inimage_name,
		char			mode,
        framelist   *   lnames,
        char        *   lines_type,
        double3     *   arcs_fwhm)
{
    qfits_header    *   fh ;
    qfits_table	    *   table ;
    qfits_col       *   col ;
    FILE            *   paf ;
    char            *   mjd_obs ;
    char                pafname[FILENAMESZ] ;
    char            *   strvar ;
    double              fwhm_med ;
	double          *   fwhm_valid ;
	int                 nb_valid ;
	char				recipe_id[FILENAMESZ] ;
	pro_catg_key        pro_category_tab ;
	pro_catg_key        pro_category_qc ;
    int                 i ;

    /* Initialize */
    if (mode == 's') {
        pro_category_tab = isaac_spec_sw_arc_coef ;
		pro_category_qc  = isaac_spec_sw_arc_qc ;
        sprintf(recipe_id, "spec_tec_arc") ;
    } else if (mode == 'l') {
        pro_category_tab = isaac_spec_lw_arc_coef ;
        pro_category_qc  = isaac_spec_lw_arc_qc ;
        sprintf(recipe_id, "lw_spec_arc") ;
    } else {
        e_error("unrecognized mode") ;
        return -1 ;
    }

   /* Write the output qfits_table table (informations) */
	table = qfits_table_new(outname, QFITS_BINTABLE, -1, 4, nb_coeffs) ;
    col = table->col ;
    for (i=0 ; i<table->nc ; i++) {
		qfits_col_fill(col, 1, sizeof(double), TFITS_BIN_TYPE_D, " ", " ", 
                " ", " ", 0, 0.0, 0, 1.0, i*sizeof(double)) ;
		col++ ;
	}
    /* Update the columns labels */
	col = table->col ;
    sprintf(col->tlabel, "Degree_of_x") ;
    col++ ;
    sprintf(col->tlabel, "Degree_of_y") ;
    col++ ;
    sprintf(col->tlabel, "poly2d_coef") ;
    col++ ;
    sprintf(col->tlabel, "WL_coefficients") ;

    /* WRITE THE OUTPUT FITS TABLE */
    /* Read the input header */
    if ((fh = qfits_header_read(inimage_name)) == NULL) {
        e_error("in writing the output fits file") ;
        qfits_table_close(table) ;
        return -1 ;
    }

    /* Prepare it for table output */
    if (isaac_header_for_table(fh) == -1) {
        e_error("in writing the output fits file") ;
        qfits_header_destroy(fh) ;
        qfits_table_close(table) ;
        return -1 ;
    }

    /* Write the PRO keywords in the header */
    if (isaac_pro_fits(fh,
                        outname,
                        "REDUCED",
                        NULL,
						pro_category_tab,
                        "OK",
                        recipe_id,
                        lnames->n,
                        lnames,
                        NULL) == -1) {
        e_error("in writing PRO keywords in output file") ;
        qfits_header_destroy(fh) ;
        qfits_table_close(table) ;
        return -1 ;
    }

    /* Write the used CATALOG in the header as PRO keyword*/
    qfits_header_add(fh, "HIERARCH ESO PRO CATALOG", lines_type, "lines", NULL);

    /* Write the HISTORY keywords with the input file names */
    if (isaac_add_files_history(fh, lnames) == -1) {
        e_warning("cannot write HISTORY keywords in out file") ;
    }

    /* Write the file on disk */
    if (qfits_save_table_hdrdump((void**)out_table, table, fh) == -1) {
        e_error("cannot write file: %s", outname) ;
        qfits_header_destroy(fh) ;
        qfits_table_close(table) ;
        return -1 ;
    }
    qfits_table_close(table) ;
    qfits_header_destroy(fh) ;

    /* WRITE THE OUTPUT PAF FILE */
    sprintf(pafname, "%s.paf", get_rootname(outname)) ;
    paf = paf_print_header( pafname,
                            "ISAAC/arcs",
                            "Arc recipe results") ;
    if (paf == NULL) {
        e_warning("cannot output PAF file") ;
    } else {
        fprintf(paf, "\n");
		/* ARCFILE */
		strvar = isaac_get_arcfile(inimage_name) ;
		if (strvar != NULL) {
			fprintf(paf, "ARCFILE \"%s\" \n", strvar) ;
		}
        /* MJD-OBS */
        mjd_obs = isaac_get_mjdobs(inimage_name) ;
        if (mjd_obs!=NULL) {
            fprintf(paf, "MJD-OBS  %s; # Obs start\n\n", mjd_obs);
        } else {
            fprintf(paf, "MJD-OBS  0.0; # Obs start unknown\n\n");
        }
        /* INSTRUME keyword  */
        strvar = isaac_get_instrument(inimage_name) ;
        if (strvar != NULL) {
            fprintf(paf, "INSTRUME \"%s\" \n", strvar) ;
        }
        /* TPL.ID  */
        strvar = isaac_get_templateid(inimage_name) ;
        if (strvar != NULL) {
            fprintf(paf, "TPL.ID  \"%s\" \n", strvar) ;
        }
        /* TPL.NEXP */
        strvar = isaac_get_numbexp(inimage_name) ;
        if (strvar != NULL) {
            fprintf(paf, "TPL.NEXP  %s \n", strvar) ;
        }
        /* DPR.CATG */
        strvar = isaac_get_dpr_catg(inimage_name) ;
        if (strvar != NULL) {
            fprintf(paf, "DPR.CATG  \"%s\" \n", strvar) ;
        }
        /* DPR.TYPE */
        strvar = isaac_get_dpr_type(inimage_name) ;
        if (strvar != NULL) {
            fprintf(paf, "DPR.TYPE  \"%s\" \n", strvar) ;
        }
        /* DPR.TECH */
        strvar = isaac_get_dpr_tech(inimage_name) ;
        if (strvar != NULL) {
            fprintf(paf, "DPR.TECH  \"%s\" \n", strvar) ;
        }
        /* Add PRO.CATG */
        fprintf(paf, "PRO.CATG  \"%s\" ;# Product category\n",
                isaac_get_pro_catg_value(pro_category_qc)) ;
        /* Add the date */
        fprintf(paf, "DATE-OBS  \"%s\" ;# Date\n",
                isaac_get_date_obs(inimage_name)) ;
        /* INS.GRAT.NAME */
        strvar = isaac_get_resolution(inimage_name) ;
        if (strvar != NULL) {
            fprintf(paf, "INS.GRAT.NAME  \"%s\" \n", strvar) ;
        }
        /* INS.GRAT.WLEN */
        fprintf(paf, "INS.GRAT.WLEN  %g \n", 
				isaac_get_central_wavelength(inimage_name)) ;
        /* INS.GRAT.ORDER */
        fprintf(paf, "INS.GRAT.ORDER %s \n", isaac_get_order(inimage_name)) ;
        /* INS.MODE */
        fprintf(paf, "INS.MODE       %s \n", isaac_get_mode(inimage_name)) ;
        /* INS.OPTI1.ID */
        fprintf(paf, "INS.OPTI1.ID   %s \n", 
				isaac_get_optical_id(inimage_name)) ;
        /* QC.LAMP */
        fprintf(paf, "QC.LAMP  \"%s\" \n", lines_type) ;
        /* QC.CENWL */
        fprintf(paf, "QC.WLEN     %g \n", 
				(double)(out_table[3][0]+out_table[3][1]*(1024/2)));
        /* QC.DISPCO1 */
        fprintf(paf, "QC.DISPCO1  %g \n", (double)out_table[3][0]);
        /* QC.DISPCO2 */
        fprintf(paf, "QC.DISPCO2  %g \n", (double)out_table[3][1]);
        /* QC.DISPCO3 */
        fprintf(paf, "QC.DISPCO3  %g \n", (double)out_table[3][2]);
        /* QC.DIST1 */
        fprintf(paf, "QC.DIST1   %g \n", (double)out_table[2][0]) ;
        /* QC.DISTX */
        fprintf(paf, "QC.DISTX   %g \n", (double)out_table[2][1]) ;
        /* QC.DISTY */
        fprintf(paf, "QC.DISTY   %g \n", (double)out_table[2][2]) ;
        /* QC.DISTXY */
        fprintf(paf, "QC.DISTXY  %g \n", (double)out_table[2][3]) ;
        /* QC.DISTXX */
        fprintf(paf, "QC.DISTXX  %g \n", (double)out_table[2][4]) ;
        /* QC.DISTYY */
        fprintf(paf, "QC.DISTYY  %g \n", (double)out_table[2][5]) ;
        /* QC.FILTER.OBS */
        strvar = isaac_get_filter(inimage_name) ;
        if (strvar != NULL) {
            fprintf(paf, "QC.FILTER.OBS        \"%s\" ;\n", strvar) ;
        }

        if (arcs_fwhm != NULL) {
            /* ARCS FWHM */
            fprintf(paf, "# ARCS FWHM\n") ;
            for(i=0 ; i<arcs_fwhm->n ; i++) {
                fprintf(paf, "QC.ARCS%d.XPOS  %.1f \n", i+1, arcs_fwhm->x[i]) ;
                fprintf(paf, "QC.ARCS%d.FWHM  %.2f \n\n", i+1, arcs_fwhm->y[i]);
            }
			/* Purge the -1 values   */
            nb_valid = 0 ;
            for (i=0 ; i<arcs_fwhm->n ; i++) 
				if (arcs_fwhm->y[i] > 0) nb_valid++ ;
            fwhm_valid = malloc(nb_valid*sizeof(double)) ;
            nb_valid = 0 ;
            for (i=0 ; i<arcs_fwhm->n ; i++)
                if (arcs_fwhm->y[i] > 0) {
                    fwhm_valid[nb_valid] = arcs_fwhm->y[i] ;
                    nb_valid++ ;
                }
			fwhm_med = double_median(fwhm_valid, nb_valid) ;
            fprintf(paf, "QC.FWHM.MED  %.2f \n", fwhm_med) ;
			free(fwhm_valid) ;
        }
        fclose(paf) ;
    }
    return 0 ;
}


/*-------------------------------------------------------------------------*/
/**
  @name     isaac_compute_distortion
  @memo     computes the distortion
  @param    in  input image
  @param	xmin
  @param	ymin	The region of interest
  @param	xmax
  @param	ymax
  @param	auto_dark_sub	Flag to automatically subtract the dark
  @param	arcs	array with arcs positions
  @return   a 2-D polynomial of size 6 double
  @doc      This function is Isaac specific. It attempts to detect a
  dark ramp and subtract it if found. See compute_distortion for a
  generic version.
 */
/*--------------------------------------------------------------------------*/
poly2d * isaac_compute_distortion(
		image_t    *   in,
		int				xmin,
		int				ymin,
		int				xmax,
		int				ymax,
		int				auto_dark_sub,
		int			*	nb_arcs,
		double		**	arcs)
{
    poly2d      	*	poly_u ;
    int         		err=0;
    int					dark_ramp_present ;
    double      		rampslope ;
    image_t    		* 	loc ;

	/* Initialize */
	poly_u = NULL ;
	dark_ramp_present = 0 ;
	rampslope = 0.0 ;
	
	/* Local copy */
    loc = image_copy(in) ;

	if (auto_dark_sub == 1) {
    	dark_ramp_present = isaac_detect_dark_ramp(loc, &rampslope);
    	if (dark_ramp_present)
    	    err=isaac_level_dark(rampslope, loc);
    	if (!err)
    	    poly_u=compute_distortion(loc, xmin, ymin, xmax, ymax,
					nb_arcs, arcs) ;
	} else if (auto_dark_sub == 0) {
		poly_u=compute_distortion(loc, xmin, ymin, xmax, ymax, nb_arcs, arcs) ;
	}

	/* Free and return */
	image_del(loc) ;
    return poly_u;
}


/*-------------------------------------------------------------------------*/
/**
  @name		isaac_distortion_estim_corr
  @memo		Estimate and correct the distortion for a batch of frames
  @param	inputs	Input images to correct	
  @param	nb_inputs	nb of input images
  @param	estimation_frame	Frame used for slit curv estimation
  @param	arc_flag	Flag to use the arc file
  @param	arc_name	Arc table that gives the slit curve distortion
  @param	startrace_flag	Flag to use the startrace file
  @param	startrace_name	Startrace table that gives the sttr distortion
  @param	xmin
  @param	ymin	Region used for disto. estimation
  @param	xmax	(C coord std : 0 -> 1023)
  @param	ymax
  @param	dark_auto_sub	Flag to automat. remove the dark ramp
  @param	status_slit_curv	slit curv. status. MODIFIED
  @param	status_startrace	startrace status. MODIFIED
  @return	corrected images or NULL if unsuccessfull correction
 */
/*--------------------------------------------------------------------------*/
image_t ** isaac_distortion_estim_corr(
		image_t		**	inputs,
		int				nb_inputs,
		char		*	estimation_frame,
		int				arc_flag,
		char		*	arc_name,
		int				startrace_flag,
		char		*	startrace_name,
		int				xmin,
		int				ymin,
		int				xmax,
		int				ymax,
		int				dark_auto_sub,
		int			*	status_slit_curv,
		int			*	status_startrace)
{
	image_t        	**  corrected ;
	poly2d          *   correct_arc ;
	poly2d          *   correct_sttr ;
	image_t        	* 	estim_image ;

	int					i, j ;

	/* Initialize */
	correct_arc = NULL ;
	correct_sttr = NULL ;

	/* FIND THE COEFFICIENTS OF THE DISTORTION */
	/* Test if arc calibr. table is available */
	if (arc_flag) {
		correct_arc = read_poly2d_from_table(arc_name) ;
		if (correct_arc == NULL) {
			e_error("cannot read the arc table") ;
			*status_slit_curv = ALGO_FAILED ;
		} else {
			*status_slit_curv = ALGO_OK ;
		}
	} else if (estimation_frame != NULL) {
		/* Estimate the disto. with the sky provided frame */
		estim_image = image_load(estimation_frame) ;
		e_comment(2, "computing distortion coefficients...");
		if ((correct_arc=isaac_compute_distortion(estim_image,
												xmin,
												ymin,
												xmax,
												ymax,
												dark_auto_sub,
												NULL,
												NULL)) == NULL) {
			e_error("in compute distortion") ;
			*status_slit_curv = ALGO_FAILED ;
		} else {
			*status_slit_curv = ALGO_OK ;
		}
		image_del(estim_image) ;
	}
	/* Test if startrace cal. table is available */
	if (startrace_flag) {
		correct_sttr = read_poly2d_from_table(startrace_name) ;
		if (correct_sttr == NULL) {
			e_error("cannot read the startrace table") ;
			*status_startrace = ALGO_FAILED ;
		} else {
			*status_startrace = ALGO_OK ;
		}
	} else {
		*status_startrace = ALGO_SKIPPED ;
	}


    /* CORRECTION OF THE DISTORTION ON ALL THE DIFFERENCE FRAMES */
    /* Use the identity polynomial in no distortion cases */
    if (*status_startrace != ALGO_OK) {
		/* Polynomial f(x,y) = y */
		correct_sttr = poly2d_build_from_string("0 1 1.0") ;
    }
    if (*status_slit_curv != ALGO_OK) {
		/* Polynomial f(x,y) = x */
		correct_arc = poly2d_build_from_string("1 0 1.0") ;
    }

    if ((*status_slit_curv == ALGO_OK) || (*status_startrace == ALGO_OK)) {
		/* Apply the transformation */
		corrected = malloc(nb_inputs*sizeof(image_t*)) ;

        /* For each difference image */
        for (i=0 ; i<nb_inputs ; i++) {
            compute_status("warping images", i, nb_inputs, 1) ;
            /* Apply the transformation on the images of the current cube */
			if ((corrected[i] = image_warp_generic(inputs[i],
					"default",
					correct_arc,
					correct_sttr)) == NULL) {
                e_error("in the correction of the distorsion") ;
                poly2d_free(correct_arc) ;
                poly2d_free(correct_sttr) ;
                for (j=0 ; j<i ; j++) {
                    image_del(corrected[j]) ;
                }
                free(corrected) ;
                if (*status_startrace == ALGO_OK) 
					*status_startrace = ALGO_FAILED ;
                if (*status_slit_curv == ALGO_OK) 
					*status_slit_curv = ALGO_FAILED ;
                return NULL ;
            }
        }
    } else {
    	e_error("cannot estimate the distortion") ;
    	poly2d_free(correct_arc) ;
    	poly2d_free(correct_sttr) ;
		return NULL ;
	}

	/* Free and return  */
	poly2d_free(correct_arc) ;
    poly2d_free(correct_sttr) ;
	return corrected ;
}


/*-------------------------------------------------------------------------*/
/**
  @name     isaac_detect_dark_ramp
  @memo     detects a dark ramp in an image
  @param    in  input image 
  @param    slope slope 
  @return   a flag signalling whether the ramp has been detected or
  not
  @doc      If the the ramp in the input image is found either the
  dark image should be subtracted or (if no dark is available)
  level_dark() may be used. By convention lo has lower y values than
  hi.
 */
/*--------------------------------------------------------------------------*/
#define IS_NB_TESTPOINTS    8
#define IS_MIN_SLOPE        0.01
#define IS_MAX_SLOPE_DIF    0.075
#define IS_MAX_FIT_EDGE_DIF 0.05
#define IS_MIN_RAMP         10.0
#define IS_MAX_MNERR        13.0
#define IS_MAX_MNERR_DIF    8.0
#define IS_MAX_INTER_DIF    20.0
#define IS_SKIPZONE         2.5

int isaac_detect_dark_ramp(image_t * in, double * slope)
{
    int         	y,
					yhi,
					ylo;
    pixelvalue  *	tmpline ;
	double3		*	testpointlo,
				*	testpointhi ;
    int         	i,
					spacing;
    double      	rampdif,
					fitslope;
    double      *	pol_coefhi,
				*	pol_coeflo ;
    double      	median[IS_NB_TESTPOINTS],
					medianerrlo,
					medianerrhi;

    if (slope==NULL || in==NULL) return 0;

    if (in->ly<IS_SKIPZONE*IS_NB_TESTPOINTS){
        e_error("detect_dark_ramp: image has %d lines, min=%d ",
                in->ly,(int)(IS_SKIPZONE*IS_NB_TESTPOINTS*2));
        return 0;
    }

    *slope=0.0;
    tmpline=malloc(in->lx*sizeof(pixelvalue));
    spacing= in->ly/(IS_SKIPZONE*IS_NB_TESTPOINTS);
    yhi= in->ly/2;
    ylo= yhi-1;
	testpointhi = double3_new(IS_NB_TESTPOINTS);
	testpointlo = double3_new(IS_NB_TESTPOINTS);

    for (i=0; i<IS_NB_TESTPOINTS; i++){
        y = yhi + i * spacing;
        memcpy(tmpline,&(in->data[y*in->lx]), in->lx*sizeof(pixelvalue));
        testpointhi->y[i] = median_pixelvalue(tmpline, in->lx);
        testpointhi->x[i] = y - in->ly/2;
        y = ylo - i * spacing;
        memcpy(tmpline,&(in->data[y*in->lx]), in->lx*sizeof(pixelvalue));
        testpointlo->y[IS_NB_TESTPOINTS-i-1]=median_pixelvalue(tmpline, in->lx);
        testpointlo->x[IS_NB_TESTPOINTS-i-1]=y;
    }
    free(tmpline);

    pol_coefhi = fit_slope_robust(testpointhi);
    pol_coeflo = fit_slope_robust(testpointlo);

    for (i=0; i<IS_NB_TESTPOINTS; i++) {
        median[i]=SQR(testpointhi->y[i]
                - pol_coefhi[0] - pol_coefhi[1] * testpointhi->x[i]);
	}

    medianerrhi = double_median(median,IS_NB_TESTPOINTS);
    for (i=0; i<IS_NB_TESTPOINTS; i++) {
        median[i]=SQR(testpointlo->y[i]
                - pol_coeflo[0] - pol_coeflo[1] * testpointlo->x[i]);
	}
    medianerrlo = double_median(median,IS_NB_TESTPOINTS);
    rampdif = testpointlo->y[IS_NB_TESTPOINTS-1] - testpointhi->y[0];

	double3_del(testpointlo);
	double3_del(testpointhi);

    if (fabs(rampdif)<IS_MIN_RAMP) {
        free(pol_coeflo);
        free(pol_coefhi);
        return 0;
    }
	
    if ((fabs(pol_coefhi[1]) < IS_MIN_SLOPE) ||
        (fabs(pol_coeflo[1]) < IS_MIN_SLOPE) ||
        (pol_coefhi[1]/pol_coeflo[1]<0.5)   ||
        (pol_coefhi[1]/pol_coeflo[1]>2.0)   ||
        (fabs(pol_coefhi[1]-pol_coeflo[1])>IS_MAX_SLOPE_DIF)) {
        free(pol_coeflo);
        free(pol_coefhi);
        return 0;
    }
    if (fabs(pol_coefhi[0]-pol_coeflo[0]) > IS_MAX_INTER_DIF){
        free(pol_coeflo);
        free(pol_coefhi);
        return 0;
    }
    if ((medianerrlo> IS_MAX_MNERR) || (medianerrhi> IS_MAX_MNERR) ||
        (fabs(medianerrlo-medianerrhi) >IS_MAX_MNERR_DIF)){
        free(pol_coeflo);
        free(pol_coefhi);
        return 0;
    }
    /* the ramp is most precisely defined at the center, include ramdif 
       in final estimate with equal weight as the two fit estimates */
    fitslope=(pol_coefhi[1]+pol_coeflo[1])/2.0;
    *slope=rampdif/(in->ly/2.0);
    if ((fabs(*slope-fitslope) > IS_MAX_FIT_EDGE_DIF) ||
        (*slope/fitslope<0.5)   ||
        (*slope/fitslope>2.0)){
        free(pol_coeflo);
        free(pol_coefhi);
        return 0;
    }
    free(pol_coeflo);
    free(pol_coefhi);
    return 1;
}

/*-------------------------------------------------------------------------*/
/**
  @name     level_dark
  @memo     subtracts a first order model of the dark current
  @param    slope   slope
  @param    in  Non dark subtracted image with two ramps in Y
  direction centered on exactly in->ly/2
  @return   0 if OK -1 if not. The ramp subtracted image is written on
  disk
 */
/*--------------------------------------------------------------------------*/
int isaac_level_dark( double    slope, image_t *in)
{
    int         i,j;
    pixelvalue  val;

    if ((in==NULL)){
        e_error("level_dark: NULL pointer provided");
        return -1;
    }
    for (j=0;j<in->ly/2;j++){
        val = slope * (j-in->ly/2) ;
        for (i=0;i<in->lx;i++)
            in->data[i+j*in->lx] -= val;
    }
    for (j=in->ly/2;j<in->ly;j++){
        val = slope * (j-in->ly);
        for (i=0;i<in->lx;i++)
            in->data[i+j*in->lx] -= val;
    }
    return 0;
}