/*----------------------------------------------------------------------------
   
   File name    :   correltest.c
   Author       :	Y. Jung
   Created on   :	July 2001
   Description  :	CONICA correlationrecipe test.

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

/*
   $Id: correltest.c,v 1.10 2002/01/23 14:15:25 ndevilla Exp $
   $Author: ndevilla $
   $Date: 2002/01/23 14:15:25 $
   $Revision $
 */

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

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include "eclipse.h"
#include "conicap_lib.h"

/*----------------------------------------------------------------------------
                                Define
 ---------------------------------------------------------------------------*/

#define KAPPA_MIN           2
#define KAPPA_MAX           50
#define KAPPA_NBTESTS       10
#define OFFSETS_PRECISION   1

/*----------------------------------------------------------------------------
                            Private functions
 ---------------------------------------------------------------------------*/

static double3 * find_reference_point_center(image_t *, double, double3 *,
		int, int, int, int) ;
static double determine_kappa(image_t *, image_t *, double, double, int) ;
static int conica_correltest_engine(char *, char *, int) ;

/*----------------------------------------------------------------------------
                            Main code
 ---------------------------------------------------------------------------*/

int conica_correltest_main(void * dict)	
{
	dictionary  *   d ;

	char            argname[10] ;
	char        *   name_i ;
	char        *   name_o ;
	int				nfiles ;
    int             refine ;

	int				errors ;
	int				i ;
	
	d = (dictionary*)dict ;
	/* Get options */
    refine = 1 ;
    
	/* 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 += conica_correltest_engine(name_i, name_o, refine) ;
        free(name_o) ;
    }
    return errors ;
}


static int conica_correltest_engine(
		char 	* 	name_i, 
		char 	* 	name_o,
        int         refine)
{
	cube_t		*	incube ;
    double          kappa ; 
    double3     *   offsets ;
    double3     *   fine_offsets ;
    image_t     *   stacked ;
    double3     *   positions ;    
    
	int				i ;

	/* Initialize */
	if ((incube = cube_load(name_i)) == NULL) {
		e_error("cannot load frames") ;
		return -1 ;
	}

    /* Determine KAPPA for star detection using the 2 first images */
    e_comment(0, "Find kappa for detection") ;
    if (fabs((kappa = determine_kappa(incube->plane[0], 
                            incube->plane[1], 
                            KAPPA_MIN,
                            KAPPA_MAX,
                            KAPPA_NBTESTS)) + 1.0) < 1e-3) {
		e_error("cannot compute kappa") ;
		cube_del(incube) ;
		return -1 ;
	}
    /* Determine offsets */ 
    e_comment(0, "Look for offsets") ;
    offsets = double3_new(incube->np) ;
    offsets->x[0] = offsets->y[0] = 0.0 ;
    for (i=1 ; i<incube->np ; i++) {
        offsets_estimates(incube->plane[0],
                          incube->plane[i],
                          &(offsets->x[i]),
                          &(offsets->y[i]),
                          kappa) ;
    	e_comment(0, "Rough offsets: %g\t%g\n", 
				offsets->x[i], offsets->y[i]) ;
    }

    if (refine) {
        e_comment(0, "Refine offsets") ;
       	positions = find_reference_point_center(incube->plane[0], 
				kappa, offsets, 15, 15, 55, 55) ;
		fine_offsets = xcorr_with_objs(incube, incube->plane[0], offsets, 
				positions, 15, 15, 55, 55) ;
        double3_del(positions) ;
        for (i=1 ; i<incube->np ; i++) {
			e_comment(0, "Fine offsets: %g\t%g\n", 
					fine_offsets->x[i], fine_offsets->y[i]) ;
            fine_offsets->x[i] = -fine_offsets->x[i] ;
            fine_offsets->y[i] = -fine_offsets->y[i] ;
        }
        stacked = cube_shiftandadd(incube, fine_offsets, "default", 0, 0) ;
        image_save_fits(stacked, "fine_stacked.fits", BPP_DEFAULT) ;
        image_del(stacked) ;
    }

    /* Shift and add */
    for (i=1 ; i<incube->np ; i++) {
        offsets->x[i] = -offsets->x[i] ;
        offsets->y[i] = -offsets->y[i] ;
    }
    stacked = cube_shiftandadd(incube, offsets, "default", 0, 0) ;
    image_save_fits(stacked, "stacked.fits", BPP_DEFAULT) ;
    image_del(stacked) ;
    
    cube_del(incube) ;
    double3_del(fine_offsets) ; 
    double3_del(offsets) ;
    return 0 ;
}


static double3 * find_reference_point_center(
		image_t		*	in,
		double			kappa,
		double3		*	offsets,
		int				sx,
		int				sy,
		int				hx,
		int				hy)
{
	image_t	*	extracted ;
	int			xmin, xmax, ymin, ymax ;
	double		offsetxmin, offsetxmax, offsetymin, offsetymax ;
	double3	*	positions ;

	int			i ;
	
	/* Determine central common zone */
	offsetxmin = offsetxmax = offsetymin = offsetymax = 0.0 ;
	for (i=0 ; i<offsets->n ; i++) {
		if (offsets->x[i] < offsetxmin) offsetxmin = offsets->x[i] ;
		if (offsets->y[i] < offsetymin) offsetymin = offsets->y[i] ;
		if (offsets->x[i] > offsetxmax) offsetxmax = offsets->x[i] ;
		if (offsets->y[i] > offsetymax) offsetymax = offsets->y[i] ;
	}
	xmin = - offsetxmin ;
	ymin = - offsetymin ;
	xmax = in->lx - 1 - offsetxmax ;
	ymax = in->ly - 1 - offsetymax ;
	
	/* Detect in the center zone */
	extracted = image_getvig(in, xmin+1, ymin+1, xmax+1, ymax+1) ;
	positions = get_points_engine(extracted, kappa, hx+sx, hy+sy, 1, 5) ;
	image_del(extracted) ;

	for (i=0 ; i<positions->n ; i++) {
		positions->x[i] += xmin ;
		positions->y[i] += ymin ;
	}
	
	return positions ;	
}

static double determine_kappa(
        image_t *   im1,
        image_t *   im2,
        double      kappa_min,
        double      kappa_max,
        int         nbtests)
{
    double      kappa ; 
    double  *   offsetsx ;
    double  *   offsetsy ;
    int         ind ;
  	int			tests_done ;
	
    int         i ;
    
    tests_done = nbtests ;
	
	/* Allocate offsets */    
    offsetsx = malloc(nbtests * sizeof(double)) ;
    offsetsy = malloc(nbtests * sizeof(double)) ;
    
    /* Compute offsets for each kappa */
    for (i=0 ; i<nbtests ; i++) {
        kappa = kappa_min + i * (kappa_max - kappa_min) / nbtests ;
        if (offsets_estimates(im1, 
								im2, 
								&offsetsx[i], 
								&offsetsy[i], 
								kappa) == -1) {
			tests_done = i ;
			break ;
		}
    }
    
    /* Determine ind */
    ind = 0 ;
    while (((fabs(offsetsx[ind]-offsetsx[ind+1]) > OFFSETS_PRECISION) || 
			(fabs(offsetsy[ind]-offsetsy[ind+1]) > OFFSETS_PRECISION)) && 
			(ind < tests_done)) ind++ ;
    
    if (ind == tests_done) {
        e_warning("Offsets do not converge") ;
        free(offsetsx) ;
        free(offsetsy) ;
        return -1.0 ;
    }
    
    free(offsetsx) ;
    free(offsetsy) ;
    kappa = kappa_min + ind * (kappa_max - kappa_min) / nbtests ;
    
    return kappa ;
}