/*-------------------------------------------------------------------------*/
/**
   @file	spectro_wave.c
   @author	N.Devillard
   @date	October 1999
   @version	$Revision: 1.21 $
   @brief	spectroscopy routines
*/
/*--------------------------------------------------------------------------*/

/*
	$Id: spectro_wave.c,v 1.21 2001/11/07 16:14:38 yjung Exp $
	$Author: yjung $
	$Date: 2001/11/07 16:14:38 $
	$Revision: 1.21 $
*/

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

#include "spectro_wave.h"
#include "spectral_lines.h"

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

static double * spectro_refine_solution(pixelvalue *, int, char	*, double, 
		double, int) ; 

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


/*-------------------------------------------------------------------------*/
/**
  @brief	compute a dispersion relation
  @param  inimage     Allocated spectroscopic image
  @param  discard_lo  number of pixels to discard at the bottom
  @param  discard_hi  (at the top)
  @param  discard_le  Number of pixels to discard on the left
  @param  discard_ri  Number of pixels to discard on the left
  @param  remove_thermal  Flag to force thermal background removal.
  @param  table_name  spectral table name (see below)
  @param  wl_min      wavelength range as [wl_min..wl_max] in angstroms
  @param  wl_max
  @return two linear coefficients such as wave = c[0] + c[1]*pix

  Compute a dispersion relation from a spectroscopic image showing some
  strong emission lines. A vital assumption is that strong emission lines
  can be seen in the image.
  
  The spectral table name is a character string. Possible values are:
  
  \begin{tabular}{ll}
  "oh"            &   OH lines \\
  "Xe"            &   Xenon lines \\
  "Ar"            &   Argon lines \\
  "Xe+Ar"         &   Xenon and Argon lines \\
  "/path/file"    &   Full pathname of an ASCII table
  \end{tabular}
  
  The latter specifies an ASCII table containing the lines you want to
  use for spectral calibration. Notice that this table must respect the
  format described in spectral_lines.c.
  
  Algorithm:
  \begin{itemize}
  \item Spectrum extraction along the spectrum direction (horizontal).
  For each column the lower and higher
  pixels are discarded, then a median value of the remaining pixels in
  the column is returned. This forms a 1d signal of same size as the
  image in the spectrum direction.
  \item If a thermal background should be removed, this is done at this
  point.
  \item If some values must be discarded (set to zero) in the input
  spectrum, they are zeroed at that point.
  \item Spectral lines are retrieved from a catalog (internal).
  \item Cross-correlation of the extracted signal and the spectral lines
  catalog at various scales. The best scale (highest cross-correlation
  factor) is kept, and the displacement giving the best cross-correlation
  factor is computed.
  \item The offset and scale are converted to wavelengths, using the
  catalog.
  \item A linear relation is computed and returned.
  \end{itemize}

  Setting parts of the input spectrum to zero enables a correct
  wavelength calibration in the thermal regime. You should provide -1 and
  -1 if you want to let this function decide for you about a correct
  zeroing interval, (0,0) if you do not want to zero the spectrum at all,
  and any other values depending on which interval you want to set to
  zero on each side of the spectrum.
  Notice that setting these values to (10,20) will set to zero the 10
  left pixels and 20 rightmost pixels of the input spectrum before
  throwing it into the cross-correlation procedure.
*/
/*---------------------------------------------------------------------------*/
/* Beginning of thermal regime in angstroms */
#define THERMAL_START		20000.0
/* Number of scales to compare */
#define XCORR_HNSCALES_PASS1        50
/* Number of pixel steps on each side to compare */
#define XCORR_WIDTH_PIX				50

#define SCALE_BEG       1.00
#define SCALE_END       1.025

#define ZEROPIX_LE		100
#define ZEROPIX_RI		100
double * spectro_compute_disprel_from_table(
	    image_t *   in,
	    int         discard_lo,
	    int         discard_hi,
		int			discard_le,
		int			discard_ri,
		int			remove_thermal,
		char	*	table_name,
	    double      wl_min,
	    double      wl_max)
{
    spectral_table	*   spt ;
    double      	*   disprel ;
	image_t			*	collapsed ;
    pixelvalue  	*   line_i ;
	pixelvalue		*	line_s ;
    pixelvalue  	*   line_t ;
    double      	*   d_t ;
    int         	    i ;
    int        		    width_i ;
    int         	    width_t ;
    double      	*   xcorr ;
    double      	*   offs ;
    double      	    xcorr_max ;
    int         	    maxpos ;
    double      	    delta ;
    int         	    nscales ;
    double      	*   scale ;
    double      	    s_min, s_max, s_step ;
    double      	    wl_from, wl_to, wl_central ;
    double      	    wl_offset, wl_scale ;
    double      	    best_scale ;
    double      	    best_offset ;

	/* Sanity checks */
	if (in==NULL) return NULL ;
	if (table_name==NULL) return NULL ;

	/* Check acceptable wavelength range */
	if ((wl_min > wl_max) || (wl_min < MIN_WAVELENGTH) ||
		(wl_max < MIN_WAVELENGTH) || (wl_min > MAX_WAVELENGTH) ||
		(wl_max > MAX_WAVELENGTH)) {
		e_error("in provided wavelength range: [%g %g]", wl_min, wl_max);
		return NULL ;
	}

    /* Get a list of lines in the image by median-collapsing it over */
    /* the horizontal direction. With a little help from the median, */
    /* should get rid of cosmics and other spiky defects. */
    collapsed = image_collapse_median(in, 0, discard_lo, discard_hi);
	line_i  = collapsed->data ;
    width_i =  in->lx ;

	if (collapsed==NULL) {
		e_error("collapsing input image: aborting wavelength calibration");
		return NULL ;
	}

	/* Remove thermal background contributions above 2 microns */
	if ((remove_thermal!=0) && (wl_max > THERMAL_START)) {
		e_comment(1, "removing thermal background");
		line_s = function1d_remove_thermalbg(line_i, width_i);
		memcpy(line_i, line_s, width_i * sizeof(pixelvalue));
		free(line_s);
	} 

	/* See if default zeroing widths have been requested */
	if ((discard_le<0) && (discard_ri<0)) {
		if (!strcmp(table_name, "oh")) {
			/* Calibrate using OH lines */
			if (wl_max > THERMAL_START) {
				/* Thermal regime: discard the right side of the signal */
				discard_le = 10 ;
				discard_ri = 512 ;
			} else {
				/* Other bands: get the default values */
				discard_le = ZEROPIX_LE ;
				discard_ri = ZEROPIX_RI ;
			}
		} else {
			/* Calibrate using standard lamps */
			discard_le = ZEROPIX_LE;
			discard_ri = ZEROPIX_RI;
		}
	}

	/* Zero the signal where requested */
	if ((discard_le>0) || (discard_ri>0)) {
		e_comment(1, "zeroing the input signal on the specified range");
		if (discard_le>0) {
			e_comment(2, "zeroing pixels [1-%d]", discard_le);
			for (i=0 ; i<discard_le ; i++) line_i[i]=0 ;
		}
		if (discard_ri>0) {
			e_comment(2, "zeroing pixels [%d-%d]", 
					width_i-discard_ri+1, width_i);
			for (i=0 ; i<discard_ri ; i++) line_i[width_i-i-1]=0 ;
		}
	}

    /* Get a list of spectral lines from official sources */
    if ((spt = spectral_table_init(table_name)) == NULL) {
		e_error("cannot initialize table: [%s]", table_name);
		image_del(collapsed);
		return NULL ;
	}

    /* Create a signal corresponding to the list of lines at a given */
    /* scale. Loop over all scales, first in a coarse pass, then in a */
    /* finer way on the second pass. */
    nscales = 2*XCORR_HNSCALES_PASS1+1 ;
    s_min   = SCALE_BEG ;
    s_max   = SCALE_END ;
    s_step  = (s_max - s_min) / (double)(nscales-1) ;
    scale = malloc(nscales * sizeof(double));
    for (i=0 ; i<nscales ; i++) scale[i] = s_min + (double)i * s_step ;

    xcorr = malloc(nscales * sizeof(double)) ;
    offs  = malloc(nscales * sizeof(double)) ;

    width_t = width_i ;
    wl_central = (wl_min + wl_max) * 0.5 ;

    for (i=0 ; i<nscales ; i++) {
        compute_status("1d x-correlation", i, nscales, 1);

        wl_from = wl_central - scale[i] * (wl_central - wl_min);
        wl_to   = wl_central + scale[i] * (wl_max - wl_central);

        if ((d_t = spectral_table_build_signal(spt, 
											wl_from, 
											wl_to, 
											width_t)) == NULL) {
			e_error("generating 1d signal from line table: aborting");
			spectral_table_destroy(spt);
			image_del(collapsed);
			free(scale);
			free(xcorr);
			free(offs);
			return NULL ;
		}
        line_t = double2pixel_array(d_t, width_t);
        free(d_t);

        xcorr[i] = function1d_xcorrelate(line_t, width_t,
										 line_i, width_i,
										 XCORR_WIDTH_PIX,
										 &delta);
        offs[i]  = delta ;
        free(line_t);
    }

    spectral_table_destroy(spt);

    /* Look for best correlation point */
    xcorr_max = xcorr[0] ;
    maxpos    = 0 ;
    for (i=1 ; i<nscales ; i++) {
        if (xcorr[i]>xcorr_max) {
            maxpos = i ;
            xcorr_max = xcorr[i] ;
        }
    }
    best_scale = scale[maxpos];
    best_offset = offs[maxpos];
    free(scale) ;
    free(xcorr);
    free(offs);

	if ((best_offset==XCORR_WIDTH_PIX) || (best_offset==-XCORR_WIDTH_PIX)) {
		e_error("cannot cross-correlate signal: aborting");
		image_del(collapsed);
		return NULL ;
	}
	
    wl_from = wl_central - best_scale * (wl_central - wl_min);
    wl_to   = wl_central + best_scale * (wl_max - wl_central);
    e_comment(0, "wave in [%g %g]", wl_from, wl_to);

    /* Change offset and scale units to wavelength */
	wl_scale = best_scale * (wl_max - wl_min) / (double)(width_i-1);
	wl_offset = wl_from - wl_scale * (1 - best_offset);

    /* Allocate results */
    disprel    = malloc(2 * sizeof(double)) ;
    disprel[0] = wl_offset;
    disprel[1] = wl_scale ;

	/* FIXME Refining is still experimental ... */
	if (debug_active()>2) {
		spectro_refine_solution(line_i, width_i,
								table_name,
								disprel[0],
								disprel[1],
								0);
	}

	image_del(collapsed);
    return disprel ;
}

#undef THERMAL_START

	
static double * spectro_refine_solution(
	pixelvalue	*	sig,
	int				npix,
	char		*	table_name,
	double			disp0,
	double			disp1,
	int				poly_deg)
{
    spectral_table  *   spt ;
	double				wave, wave_min, wave_max ;
	int					i, j;
	int					i_min, i_max ;
	double			*	locmax ;
	double			*	lines ;
	double			*	locmax_f ;
	int					n_ok ;
	int					flag_ok ;
	int					nlines ;
	int					hs = 8 ;

    /*
     * Get a list of spectral lines from official sources
     */
    spt = spectral_table_init(table_name);
    if (spt==NULL) {
        e_error("cannot initialize table: [%s]", table_name);
        return NULL ;
    }

	/* Locate indexes of min and max wavelengths */
	wave_min = disp0 + disp1 ;
	wave_max = disp0 + (double)npix * disp1 ;

	i=0 ;
	while (((spt->lines[i]).wavel<wave_min) && (i<spt->nlines)) i++ ;
	i_min = i ;
	while (((spt->lines[i]).wavel<wave_max) && (i<spt->nlines)) i++ ;
	i_max = i ;
	if (i_min==spt->nlines) {
		e_error("wavelength range [%g %g] outside known bounds",
				wave_min, wave_max);
		return NULL ;
	}
	nlines = i_max - i_min + 1 ;

	/* Store lines in a double array */
	lines = malloc(nlines * sizeof(double));
	for (i=i_min ; i<=i_max ; i++) {
		lines[i-i_min] = (double)(spt->lines[i]).wavel ;
	}
	spectral_table_destroy(spt);

	/* Loop on all lines that should be present in the signal */
	locmax = malloc(nlines * sizeof(double)) ;
	for (i=0 ; i<nlines ; i++) {
		wave = lines[i];
		/* Look for index of corresponding wavelength in signal */
		j = (int)((wave-disp0)/disp1 - 1.0);
		if ((j>hs) && (j<(npix-hs-1))) {
			locmax[i] = function1d_find_locmax(sig, npix, j, hs);
		} else {
			locmax[i] = -1.0 ;
		}
	}
	/* Remove double lines */
	n_ok = 0 ;
	for (i=0 ; i<nlines ; i++) {
		flag_ok = 1 ;
		if (locmax[i]<0.0) {
			flag_ok = 0 ;
		} else if (i>0) {
			if (fabs(locmax[i]-locmax[i-1])<1e-4) {
				flag_ok=0;
			}
		}
		if (flag_ok) n_ok ++ ;
	}
	locmax_f = malloc(n_ok * sizeof(double));
	n_ok = 0 ;
	for (i=0 ; i<nlines ; i++) {
		flag_ok = 1 ;
		if (locmax[i]<0.0) {
			flag_ok = 0 ;
		} else if (i>0) {
			if (fabs(locmax[i]-locmax[i-1])<1e-4) {
				flag_ok=0;
			}
		}
		if (flag_ok) {
			locmax_f[n_ok] = locmax[i];
			n_ok ++ ;
		}
	}
	free(locmax);
	/* Print out new list */
	printf("-------------\n");
	for (i=0 ; i<n_ok ; i++) {
		/* printf("%g\t%g\n", lines[i], disp0 + disp1 *(1+locmax_f[i])); */
		printf("%g\t%g\n", locmax_f[i], lines[i]);
	}
	printf("-------------\n");

	free(lines);
	free(locmax_f);
	return NULL ;
}