/*----------------------------------------------------------------------------
   
   File name    :   slitpos.c
   Author       :	Y. Jung
   Created on   :	Feb. 2001
   Description  :	Slit position 

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

/*
	$Id: slitpos.c,v 1.39 2002/03/04 14:11:01 yjung Exp $
	$Author: yjung $
	$Date: 2002/03/04 14:11:01 $
	$Revision: 1.39 $
 */

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

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

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

#define MAX_NB_EROSIONS     1024
#define KERNEL_SIZE_Y       5

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

static int slitpos_engine(char *, char *, int) ;
static int slitpos_write_paffile(char *, char *, double, double, double) ;
static int slitpos_write_outfile(char *, int, double, double3 **, 
		framelist *, int) ;
static int slitpos_find_edges_one_line(image_t *, int, int *, int *) ;
static int slitpos_find_vert_slit_ends(image_t *, int, int *, int *) ;
static int slitpos_find_vert_pos(image_t *, int *) ;
static double3 ** slitpos_analysis(image_t *, int, double *, int *) ;

/*----------------------------------------------------------------------------
                            Main code
 ---------------------------------------------------------------------------*/
int isaac_slitpos_main(void * dict)
{
	dictionary	*	d ;
		
	int				slit_max_width ;
	
	char			argname[10] ;
	char    	*	name_i ;
    char    	*	name_o ;
    int     		nfiles ;

	int				errors ;
	int				i ;
	 
	d = (dictionary*)dict ;
   	/* Get options */
	slit_max_width   = dictionary_getint(d, "arg.max_width", 20) ;

	/* 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 += slitpos_engine(name_i,
                            		name_o,
									slit_max_width) ;
		free(name_o) ;
	}
	return errors ;
}


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


/*-------------------------------------------------------------------------*/
/**
  @name     slitpos_engine
  @memo     engine for the is_spec_slitpos recipe
  @param    inname  input file
  @param    outname output file name
  @param    slit_max_width  maximum slit width
  @return   int -1 in error case
  @doc      Find out if the first frame is a dark and subtract it from 
            slit images. Then detect the slit on each slit images
 */
/*--------------------------------------------------------------------------*/
static int slitpos_engine(
        char    *   inname,
        char    *   outname,
        int         slit_max_width)
{
	cube_t		*	images ;
	framelist	*	filenames ;
	char		*	current_file ;
	char			output_name[FILENAMESZ] ;
	char		*	mode ;
	int				dark_present ;
	image_t		*	dark ;
	int				first_frame_ind ;
	double				slit_angle ;
	int				slit_length ;
	double			xcenter,
					ycenter ;

	double3		**	out_table ;				
	
	int				i,
					j ;
	
	/* Initialize */
	dark = NULL ;
	dark_present = 0 ;
	first_frame_ind = 0 ;
	 
	/* Prepare the cube for reduction */

	/* If the input file is a FITS file */
	if (is_fits_file(inname)) {
		images = cube_load(inname) ;
		if (images == NULL) {
			e_error("cannot load FITS file [%s]", inname) ;
			return -1 ;
		}
	
		/* Fill filenames */
		filenames = framelist_new(1) ;
		filenames->name[0] = strdup(inname) ;
		
	/* If the input file is an ASCII file */
	} else {
	
		/* Read the input ASCII file  */
		filenames = framelist_load(inname) ;
		if (filenames == NULL) {
			e_error("cannot read the input ASCII file: [%s]", inname) ;
			return -1 ;
		}
	
		/* Find out if the first frame is a dark and load it */
		current_file = filenames->name[0] ;
		mode = isaac_get_mode(current_file) ;
		if (!strcmp(mode, "SW_DARK") || !strcmp(mode, "LW_DARK")) {
			e_comment(1, "dark present: %s", current_file) ;
			dark_present = 1 ;
			dark = image_load(current_file) ;
			first_frame_ind = 1 ;
			
			/* Test if there is not only the dark */
			if (filenames->n < 2) {
				e_error("only a DARK frame in the list - abort") ;
				framelist_del(filenames) ;
				image_del(dark) ;
				return -1 ;
			}
		} else {
			e_comment(1, "No dark present") ;
		}
		
		/* Load the cube */
		images = cube_load(inname) ;
		if (images == NULL) {
			e_error("cannot load ASCII file [%s]", inname) ;
			framelist_del(filenames) ;
			if (dark_present) image_del(dark) ;
			return -1 ;
		}

		/* Subtract the dark if present */
		if (dark_present) {
			if (cube_sub_im(images, dark) == -1) {
				e_error("cannot subtract dark - abort") ;
				image_del(dark) ;
				cube_del(images) ;
				framelist_del(filenames) ;
				return -1 ;
			}
			image_del(dark) ;
		}
	}
	
	/* Loop on all the slit images */
	for (i=first_frame_ind ; i<images->np ; i++) {
		e_comment(1, "Slit image no %d", i+1) ;
		/* Slit analysis */
		if ((out_table = slitpos_analysis(images->plane[i],
								slit_max_width,
								&slit_angle,
								&slit_length)) == NULL) {
			e_error("in slit position analysis: [%s]", filenames->name[i]) ;
		} else {
			sprintf(output_name, "%s_%d.tfits", outname, i+1) ;
			/* Write the output files (TFITS & PAF) */
			if (slitpos_write_outfile(output_name,
										slit_length,
										slit_angle,
										out_table,
										filenames, 
										i) == -1) {
				e_error("in writing the output FITS table") ;
				framelist_del(filenames) ;
				cube_del(images) ;
				for (j=0 ; j<3 ; j++) {
					double3_del(out_table[j]) ;
				}
				free(out_table) ;
				return -1 ;
			}
			
			xcenter = (out_table[1]->x[0]+out_table[1]->x[slit_length-1])/2.0 ;
			ycenter = (out_table[1]->y[0]+out_table[1]->y[slit_length-1])/2.0 ;
			for (j=0 ; j<3 ; j++) {
				double3_del(out_table[j]) ;
			}
			free(out_table) ;

			/* Write the output PAF file */
			sprintf(output_name, "%s.paf", get_rootname(output_name)) ;
			if (slitpos_write_paffile(output_name,
									filenames->name[i],
									xcenter,
									ycenter,
									slit_angle) == -1) {
				e_error("in writing the output PAF file") ;
				framelist_del(filenames) ;
				cube_del(images) ;
				return -1 ;
			}
		}
	}

	/* Free and return */
	framelist_del(filenames) ;
	cube_del(images) ;
	return 0 ;
}


/*-------------------------------------------------------------------------*/
/**
  @name		slitpos_write_paffile
  @memo		Write the output PAF file 
  @param	outname	name of the out file
  @param	inname	name of the input FITS file
  @param	xcenter	X coordinate of the slit center
  @param	ycenter Y coordinate of the slit center
  @param	slit_angle	slit angle with the horizontal in degrees
  @return	-1 in error case, 0 otherwise
 */
/*--------------------------------------------------------------------------*/
static int slitpos_write_paffile(
		char	*	outname,
		char	*	inname,
		double		xcenter,
		double		ycenter,
		double		slit_angle)
{
	FILE	*	paf ;

	char	*	mjd_obs ;
	char	*	strvar ;
	
    paf = paf_print_header( outname,
                            "ISAAC/slitpos",
                            "Slit position recipe results") ;
   if (paf == NULL) {
       e_warning("cannot output PAF file") ;
   } else {
	   fprintf(paf, "\n");
	   /* ARCFILE */
	   strvar = isaac_get_arcfile(inname) ;
	   if (strvar != NULL) {
		   fprintf(paf, "ARCFILE   \"%s\"    \n", strvar) ;
	   }
       /* MJD-OBS */
	   mjd_obs = isaac_get_mjdobs(inname) ;
       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(inname) ;
	   if (strvar != NULL) {
		   fprintf(paf, "INSTRUME \"%s\" \n", strvar) ; 
	   } 
	   /* TPL.ID  */
	   strvar = isaac_get_templateid(inname) ;
	   if (strvar != NULL) {
		   fprintf(paf, "TPL.ID  \"%s\" \n", strvar) ;
	   }
	   /* TPL.NEXP */
	   strvar = isaac_get_numbexp(inname) ;
	   if (strvar != NULL) {
		   fprintf(paf, "TPL.NEXP  %s \n", strvar) ;
	   } 
	   /* DPR.CATG */
	   strvar = isaac_get_dpr_catg(inname) ;
	   if (strvar != NULL) {
		   fprintf(paf, "DPR.CATG  \"%s\" \n", strvar) ;
	   }
	   /* DPR.TYPE */
	   strvar = isaac_get_dpr_type(inname) ;
	   if (strvar != NULL) {
		   fprintf(paf, "DPR.TYPE  \"%s\" \n", strvar) ;
	   }
	   /* DPR.TECH */
	   strvar = isaac_get_dpr_tech(inname) ;
	   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(isaac_spec_slitpos_qc)) ;
        /* Add the date */
        fprintf(paf, "DATE-OBS \"%s\" ;# Date\n",
                isaac_get_date_obs(inname)) ;
	   /* INS.OPTI1.ID */
	   strvar = isaac_get_optical_id(inname) ;
	   if (strvar != NULL) {
		   fprintf(paf, "INS.OPTI1.ID  \"%s\" \n", strvar) ;
	   }
	   /* QC.SLIT.XPOS  */
	   fprintf(paf, "QC.SLIT.XPOS  %g \n", xcenter) ;
	   
	   /* QC.SLIT.YPOS */
	   fprintf(paf, "QC.SLIT.YPOS  %g \n", ycenter) ;
	   
	   /* QC.SLIT.ANGLE */
	   fprintf(paf, "QC.SLIT.POSANG  %g \n", slit_angle) ;
	  
	   fclose(paf) ;
	   e_comment(0, "file [%s] produced", outname) ; 	
   }

   return 0 ;
}

/*-------------------------------------------------------------------------*/
/**
  @name		slitpos_write_outfile
  @memo		Write the output FITS table		
  @param	outname	name of out file	
  @param	slit_length	length of the out table
  @param	slit_angle	angle of the slit
  @param	out_table	result of slitpos_analysis()
  @param	filenames	list of input file names (charmatrix)
  @param	int			file id 
  @return	-1 in error case, 0 otherwise
 */
/*--------------------------------------------------------------------------*/
static int slitpos_write_outfile(
	char		*	outname,
	int				slit_length,
	double			slit_angle,
	double3		**	out_table,
	framelist	*	filenames,
	int				file_id)
{
	qfits_header	*	fh ;
	qfits_table		*	table ;
	qfits_col		*	col ;
	int					i, j ;
	double			**	data ;
	
	/* Write the output qfits_table table (informations) */
	table = qfits_table_new(outname, QFITS_BINTABLE, -1, 4, slit_length) ;
    col = table->col ;
    for (i=0 ; i<table->nc ; i++) {
    	qfits_col_fill(col, 1, sizeof(double), TFITS_BIN_TYPE_D, "pixel", " ", 
                " ", " ", 0, 0.0, 0, 1.0, i*sizeof(double)) ;
		col++ ;
    }
	col = table->col ;
	strcpy(col->tlabel, "Y") ;
	col++ ;
	strcpy(col->tlabel, "LEFT_POSITION") ;
	col++ ;
	strcpy(col->tlabel, "CENTER_POSITION") ;
	col++ ;
	strcpy(col->tlabel, "RIGHT_POSITION") ;

    /* Put out_table in data */
    data = malloc(table->nc * sizeof(double*)) ;
    for (i=0 ; i<table->nc ; i++) {
        data[i] = malloc(table->nr * sizeof(double)) ;
    }
	for (j=0 ; j<table->nr ; j++) {
		data[0][j] = out_table[0]->y[j] ;
		data[1][j] = out_table[0]->x[j] ;
		data[2][j] = out_table[1]->x[j] ;
		data[3][j] = out_table[2]->x[j] ;
	}

	/* WRITE THE OUTPUT FILE */
	/* Read the input header */
	if ((fh = qfits_header_read(filenames->name[file_id])) == NULL) {
        e_error("in writing the output fits file") ;
		for (i=0 ; i<table->nc ; i++) free(data[i]) ;
        free(data) ;
		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") ;
		for (i=0 ; i<table->nc ; i++) free(data[i]) ;
        free(data) ;
        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,
                        isaac_spec_slitpos_table,
                        "OK",
                        "img_tec_slitposition",
                        filenames->n,
						filenames,
						NULL) == -1) {
        e_error("in writing PRO keywords in output file") ;
		for (i=0 ; i<table->nc ; i++) free(data[i]) ;
        free(data) ;
        qfits_header_destroy(fh) ;
		qfits_table_close(table) ;
        return -1 ;
    }

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

	/* Write the file on disk */
    if (qfits_save_table_hdrdump((void**)data, table, fh) == -1) {
        e_error("cannot write file: %s", outname) ;
        qfits_header_destroy(fh) ;
		qfits_table_close(table) ;
		for (i=0 ; i<table->nc ; i++) free(data[i]) ;
        free(data) ;
        return -1 ;
    }
	qfits_table_close(table) ;
    qfits_header_destroy(fh) ;
	for (i=0 ; i<table->nc ; i++) free(data[i]) ;
    free(data) ;

    e_comment(0, "File [%s] produced", outname) ;

    /* Write slit_angle slit_length and the slit center coordinates */ 
	/* on the screen */
	e_comment(0,"Slit angle with horizontal (in deg): %g\n", slit_angle) ;
	e_comment(0,"Slit center coordinates: (%g, %g)\n", 
					(out_table[1]->x[0]+out_table[1]->x[slit_length-1])/2.0,
				(out_table[1]->y[0]+out_table[1]->y[slit_length-1])/2.0) ;
	e_comment(0,"Slit length in pixels: %d\n", slit_length) ;
	e_comment(0,"file [%s] produced", outname) ;

	return 0 ;
}

/*-------------------------------------------------------------------------*/
/**
  @name		slitpos_find_edges_one_line
  @memo		return the first pixel higher than avg starting from the
  			left and from the right of the line	
  @param	inimage	input image
  @param	line_pos	line position
  @param	left_pos	pointer to left position
  @param	right_pos	pointer to right position
  @return	-1 in error case, 0 otherwise
  @doc	Positions in C coordinates (first pixel position is 0)
 */
/*--------------------------------------------------------------------------*/
static int slitpos_find_edges_one_line(
        image_t    *   inimage,
        int             line_pos,
        int         *   left_pos,
        int         *   right_pos)
{
    image_stats *   statistics ;
    pixelvalue      threshold ;
	int				zone[4] ;

    int             i ;

    /* Find the threshold */
	zone[0] = 1 ;
	zone[1] = inimage->lx ;
	zone[2] = line_pos+1 ;
	zone[3] = line_pos+1 ;
	
    if ((statistics = image_getstats_opts(inimage,
										   NULL,
										   NULL,
										   zone))==NULL) {
        e_error("cannot get the statistics of [%d]th line",
                line_pos+1);
        return -1 ;
    }
    threshold = statistics->avg_pix ;
    free(statistics) ;

    /* Detect the left edge */
    i = 0 ;
    while ((inimage->data[line_pos*inimage->lx+i] < threshold)
            && (i < inimage->lx)) i++ ;
    *left_pos = i ;

    /* Detect the right edge */
    i = inimage->lx - 1 ;
    while ((inimage->data[line_pos*inimage->lx+i] < threshold)
            && (i >= 0 )) i-- ;
    *right_pos = i ;

    return 0 ;
}


/*-------------------------------------------------------------------------*/
/**
  @name		slitpos_find_vert_slit_ends
  @memo		Find the ends of a vertical slit (y coordinates in FITS
  convention)
  @param	in	input image
  @param	kernel_size	vertical kernel size
  @param	bot_slit_y	bottom slit y position
  @param	top_slit_y	top slit y position
  @return	-1 in error case, 0 otherwise
  @doc		The input image as to be as thin as possible to contain
  'only' the slit
 */
/*--------------------------------------------------------------------------*/
static int slitpos_find_vert_slit_ends(
        image_t    *   in,
        int             kernel_size,
        int         *   bot_slit_y,
        int         *   top_slit_y)
{
    image_stats     *   statistics ;
    pixelmap       	*  	binary ;
    pixelmap       	*  	kernel ;
    intimage	    *  	label_image ;
	int					nobj ;
	int                 erosions_nb ;
    int                 i ;


    /* Threshold to have a binary image */
    if ((statistics = image_getstats(in)) == NULL) {
        e_error("unable to get the image stats") ;
        return -1 ;
    }
    binary = image_threshold2pixelmap(in,
									 (double)statistics->avg_pix,
									 (double)statistics->max_pix);
	free(statistics) ;
	if (binary==NULL) {
        e_error("failed while binarizing the image") ;
        return -1 ;
    }

    /* Define the kernel for morpho operations */
    if ((kernel = pixelmap_new(1,kernel_size)) == NULL) {
        e_error("cannot create the kernel") ;
        pixelmap_del(binary) ;
        return -1 ;
    }

    /* Erode until there is 1 object left in the image */
    i = 0 ;
    label_image = intimage_labelize_pixelmap(binary, &nobj) ;
	intimage_del(label_image);
    while (nobj>1) {
		if (pixelmap_morpho_erosion_k(binary, kernel) != 0) {
            e_error("cannot erode") ;
            pixelmap_del(binary) ;
            pixelmap_del(kernel) ;
            return -1 ;
        }
        i++ ;
        if (i >= MAX_NB_EROSIONS) {
            e_error("max number of erosions reached : [%d] - aborting",
                    i) ;
            pixelmap_del(binary) ;
            pixelmap_del(kernel) ;
            return -1 ;
        }
        label_image = intimage_labelize_pixelmap(binary, &nobj) ;
		intimage_del(label_image);
    }
    if (nobj<=0) {
        e_error("no detected slit") ;
        pixelmap_del(binary) ;
        pixelmap_del(kernel) ;
        return -1 ;
    }
    erosions_nb = i ;

    /* Reconstruct the slit with dilatations */
    if (erosions_nb > 0) {
        pixelmap_del(kernel) ;
        if ((kernel = pixelmap_new(1,
                            (kernel_size-1)*erosions_nb+1)) == NULL) {
            e_error("cannot create the kernel") ;
            pixelmap_del(binary) ;
            return -1 ;
        }
        if (pixelmap_morpho_dilation_k(binary, kernel) != 0) {
            e_error("cannot dilate") ;
            pixelmap_del(binary) ;
            pixelmap_del(kernel) ;
            return -1 ;
        }
    }
    pixelmap_del(kernel) ;

    /* Find the ends of the slit */
    i = 0 ;
    while (binary->data[i] != 1) i++ ;
    *bot_slit_y = (int)(i/binary->lx + 1) ;
    i = (binary->lx * binary->ly) - 1 ;
    while (binary->data[i] != 1) i-- ;
    *top_slit_y = (int)(i/binary->lx + 1) ;

    /* Free and return */
    pixelmap_del(binary) ;
    return 0 ;
}



/*-------------------------------------------------------------------------*/
/**
  @name		slitpos_find_vert_pos
  @memo		Find a vertical slit position (x coordinate of the slit)
  @param	in	input image
  @param	slit_pos	pointre to the searched position
  @return	-1 in error case, 0 otherwise
  @doc		Coordinate given in FITS convention (ll is (1,1))
 */
/*--------------------------------------------------------------------------*/
static int slitpos_find_vert_pos(
        image_t    *   in,
        int         *   slit_pos)
{
    image_t    *   filtered ;
    image_t    *   image1D ;
    image_stats *   statistics ;

    /* Filter the input image to 'erase' the bad pixels */
    if ((filtered = image_filter_median(in)) == NULL) {
        e_error("unable to filter the image") ;
        return -1 ;
    }

    /* Collapse the image to a horizontal 1D image */
    if ((image1D = image_collapse(filtered, 0)) == NULL) {
        e_error("unable to collapse the image") ;
        image_del(filtered) ;
        return -1 ;
    }

    /* Search the max of the 1D image to identify the slit position */
    if ((statistics = image_getstats(image1D)) == NULL) {
        e_error("unable to get image statistics") ;
        image_del(filtered) ;
        image_del(image1D) ;
        return -1 ;
    }
    *slit_pos = statistics->max_x + 1;

    /* Free and return */
    free(statistics) ;
    image_del(image1D) ;
    image_del(filtered) ;
    return 0 ;
}


/*-------------------------------------------------------------------------*/
/**
  @name		slitpos_analysis	
  @memo		detect the slit position, detect its ends, extract a thin
  image containing only the slit and find its edges
  @param	inimage	input image
  @param	slit_max_width maximum slit width
  @param	slit_angle	pointer to the angle horizontal-slit
  @param	slit_length	pointer to the slit length
  @return	ptr to 3 double3 objects.
  @doc

  This function can be used for horizontal or vertical slits.

  This function returns 3 double3 objects:
  \begin{itemize}
  \item Left or Lower edge of the slit
  \item Center of the slit
  \item Right or Upper edge of the slit
  \end{itemize}

  NB: Coordinates use FITS convention.
 */
/*--------------------------------------------------------------------------*/
static double3 ** slitpos_analysis(
        image_t	*   inimage,
        int        		slit_max_width,
        double  	*  	slit_angle,
        int     	*   slit_length)
{
    image_t        *   filtered ;
    image_t        *   extr_fits ;
    int                 slit_pos ;
    int                 spec_size ;
    int                 slit_size ;
    int                 slit_top_y ;
    int                 slit_bot_y ;
	double3			**	out ;
	double3			*	slit_l,
					*	slit_c,
					*	slit_r ;
    int                 right_pos ;
    int                 left_pos ;
    double          *   coeff_r ;
    double          *   coeff_l ;

    int                 i ;

    /* Flip the image or not to have a vertical slit */
	spec_size = inimage->lx ;
	slit_size = inimage->ly ;

    /* Find the position of the slit */
    if ((slitpos_find_vert_pos(inimage, &slit_pos)) == -1) {
        e_error("cannot find the slit position") ;
        return NULL ;
    }

    /* Filter the input image to 'erase' the bad pixels */
    if ((filtered = image_filter_median(inimage)) == NULL) {
        e_error("unable to filter the image") ;
        return NULL ;
    }

    /* Extract a thin image containing the slit */
    if ((extr_fits = image_getvig(filtered,
                            (int)(slit_pos-(int)(slit_max_width/2)),
                            1,
                            (int)(slit_pos+(int)(slit_max_width/2)),
                            (int)(slit_size))) == NULL) {
        e_error("unable to extract the thin image") ;
        image_del(filtered) ;
        return NULL ;
    }

    /* Find the ends of the slit */
    if ((slitpos_find_vert_slit_ends(extr_fits,
                            KERNEL_SIZE_Y,
                            &slit_bot_y,
                            &slit_top_y)) == -1) {
        e_error("cannot find the ends of the slit") ;
        image_del(filtered) ;
        image_del(extr_fits) ;
        return NULL ;
    }
    image_del(extr_fits) ;
    *slit_length = slit_top_y - slit_bot_y ;

    /* Extract an image with exactly the slit */
    if ((extr_fits = image_getvig(filtered,
                            (int)(slit_pos-(int)(slit_max_width/2)),
                            slit_bot_y,
                            (int)(slit_pos+(int)(slit_max_width/2)),
                            slit_top_y)) == NULL) {
        e_error("cannot extract a thin image") ;
        image_del(filtered) ;
        return NULL ;
    }
    image_del(filtered) ;

    /* Allocate the slit position arrays */
    slit_l = double3_new(*slit_length);
    slit_c = double3_new(*slit_length);
    slit_r = double3_new(*slit_length);

    /* Find the edges of the slit */
    for (i=0 ; i<*slit_length ; i++) {
        if ((slitpos_find_edges_one_line(extr_fits,
                                i,
                                &left_pos,
                                &right_pos)) == -1) {
            e_error("cannot find the edges of the [%d]th line", i+1) ;
            image_del(extr_fits) ;
			double3_del(slit_l);
			double3_del(slit_c);
			double3_del(slit_r);
            return NULL ;
        }
        /* Store the edges as horizontal lines for the fit */

        slit_l->y[i] = (double)left_pos ;
        slit_l->x[i] = (double)(i+slit_bot_y-1) ;

        slit_r->y[i] = (double)right_pos ;
        slit_r->x[i] = (double)(i+slit_bot_y-1) ;
    }
    image_del(extr_fits) ;

    /* Linear regression to find the edges */
    coeff_l = fit_slope_robust(slit_l);
    coeff_r = fit_slope_robust(slit_r);

    /* Rewrite the edges in the out table, and write the center */
    for (i=0 ; i<*slit_length ; i++) {
		slit_l->y[i] = slit_c->y[i] = slit_r->y[i] = (double)(i+slit_bot_y);
		slit_l->x[i] =
			coeff_l[0]+coeff_l[1]*slit_l->y[i]+(slit_pos-slit_max_width/2);
		slit_r->x[i] =
			coeff_r[0]+coeff_r[1]*slit_r->y[i]+(slit_pos-slit_max_width/2);
		slit_c->x[i] = (double)((slit_l->x[i] + slit_r->x[i])/2.) ;
    }
    free(coeff_r) ;
    free(coeff_l) ;

    /* Find the slit angle in degrees with the horizontal axis   */
    *slit_angle = (double)(atan((slit_c->y[*slit_length-1] -
                        slit_c->y[0]) / (slit_c->x[*slit_length-1] -
                        slit_c->x[0]))) ;
    *slit_angle *= 180.0/M_PI ;
    
	/* Free and return   */
    out = malloc(3*sizeof(double3*)) ;
	out[0] = slit_l ;
	out[1] = slit_c ;
	out[2] = slit_r ;
    return out ;
}