/* @(#)wa_comp.c	17.1.1.1 (ES0-DMD) 01/25/02 17:21:40 */
/*===========================================================================
  Copyright (C) 1995 European Southern Observatory (ESO)
 
  This program is free software; you can redistribute it and/or 
  modify it under the terms of the GNU General Public License as 
  published by the Free Software Foundation; either version 2 of 
  the License, or (at your option) any later version.
 
  This program is distributed in the hope that it will be useful,
  but WITHOUT ANY WARRANTY; without even the implied warranty of
  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  GNU General Public License for more details.
 
  You should have received a copy of the GNU General Public 
  License along with this program; if not, write to the Free 
  Software Foundation, Inc., 675 Massachusetss Ave, Cambridge, 
  MA 02139, USA.
 
  Corresponding concerning ESO-MIDAS should be addressed as follows:
	Internet e-mail: midas@eso.org
	Postal address: European Southern Observatory
			Data Management Division 
			Karl-Schwarzschild-Strasse 2
			D 85748 Garching bei Muenchen 
			GERMANY
===========================================================================*/

/******************************************************************************
**            Copyright (C) 1993 by European Southern Observatory
*******************************************************************************
**
**    UNIT
**
**    Version: 17.1.1.1
**
**    Author: Jean-Luc Starck
**
**    Date:  02/01/25
**    
**    File:  wa_comp.c
**
*******************************************************************************
**
**    DESCRIPTION   Computes the correlation coefficients and the signal to
**    -----------   noise ratio between two images.  The comparison is done 
**                  in the wavelet space. 
**                  At each scale we compute the SNR and the
**                  correlation. We get two curves, 
**                        . (coeff_correlation - frequence)
**                        . (SNR - frequence)
**
**    PARAMETRES    
**    ----------
**         image 1 (keyword:IN_A)
**
**         image 2 (keyword:IN_B)
**
**         Number of scales of the comparison (keyword:INPUTI[1])
**
**         Comparison parameter N_Sigma (keyword:INPUTR[1]):
**            if N_Sigma > 0, the comparison is not done with
**            the full image, but only with the structure in the
**            wavelet space > N_Sigma * sigma(scale)
**
**         Input-Output correlation table name (keyword:OUT_A)
**
**         Input-Output signal to noise ratio table name (keyword:OUT_B)
**
**    RESULTS two table is created which contains the result of the comparisons
**    ------- if the tables exists, columns are added.
**            the first line corresponds to comparison between the two images
**            the others line corresponds to comparison between the wavelet
**            planes.
******************************************************************************/

static char sccsid[] = "@(#)wa_comp.c 17.1.1.1 02/01/25 ESO 1993 @(#)";

#include <stdio.h>
#include <math.h>
#include <string.h>

#include<midas_def.h>

#include "Def_Math.h"
#include "Def_Mem.h"
#include "Def_Wavelet.h"

#define VISU_PARAM FALSE

#define TAB_MAX_SCALE 11

float N_Sigma;

/***************************************************************************/

compare_plan (Imag_1, Imag_2, Nl, Nc, Snr_Db, Correl)
float *Imag_1, *Imag_2, *Snr_Db, *Correl;
int Nl, Nc;
/* computes information between two images

   Snr = variance (Imag_1) / variance de (Imag_1 - Imag_2)
   Snr_Db = 10 log10 (Snr)
   Correl = correlation des images Imag_1 et Imag_2
*/
{
    int i;
    float Snr,S1,S2;
    float Error;
    float Sum_X2,Sum_Y2,Sum_XY;
    float Moy = 0., Ecart = 0.;
    float Moy_Err = 0., Ecart_Err = 0., Noise;
    int Size = 0;
    float lib_mat_ecart_type();

    Noise = N_Sigma * lib_mat_ecart_type(Imag_1, Nl, Nc);
    Sum_X2 = Sum_Y2 = Sum_XY = 0.;
    for (i = 0; i < Nl * Nc; i++)
    {
        if (fabs(Imag_1 [i]) > Noise)
	{
           /* Calcul Correlation */
           Sum_X2 += Imag_1 [i] * Imag_1 [i];
           Sum_Y2 += Imag_2 [i] * Imag_2 [i];
           Sum_XY += Imag_1 [i] * Imag_2 [i];

           /* Calcul of the standard deviation of Imag_1 */
           Moy += Imag_1 [i];
           Ecart += Imag_1 [i] * Imag_1 [i];

           /* Calcul of the standard deviation of the error */
           Error = Imag_1[i] - Imag_2[i];
           Moy_Err += Error;
           Ecart_Err += Error * Error;
           Size ++;
	 }
    }
    Snr = (float) Size / (float)(Nl*Nc) * 100.;
#if VISU_PARAM
    printf ("Signif = %5.2f %%\n", Snr);
#endif

    /* Correlation */
    *Correl = Sum_XY / sqrt (Sum_X2 * Sum_Y2);

    /* variance of Imag_1 */
    Moy /= (float) Size;
    Ecart /= (float) Size;
    S1 = Ecart - Moy*Moy;

    /* variance of the error */
    Moy_Err /= (float) Size;
    Ecart_Err /= (float) Size;
    S2 = Ecart_Err - Moy_Err*Moy_Err;

    /* signal to noise ration */
    Snr = S1 / S2;

    /* signal to noise ration (dB) */
    *Snr_Db = 10. * log10 (Snr);
}

/***************************************************************************/

compare_wavelet (Imag1, Imag2, Nl, Nc, Nbr_Plan, Tab_Snr, Tab_Correl)
float *Imag1, *Imag2;
int Nbr_Plan, Nl, Nc;
float *Tab_Snr, *Tab_Correl;
{
    int i, Nbr_Etap;
    int Type_Transform = 5;
    float Moy, Fc = 0.5, *Ptr1, *Ptr2;
    wave_transf_des Wavelet_1, Wavelet_2;

    compare_plan (Imag1, Imag2, Nl, Nc, Tab_Snr, Tab_Correl);

    wavelet_transform_data (Imag1, Nl, Nc, &Wavelet_1,
                            Type_Transform, Fc, Nbr_Plan);
    wavelet_transform_data (Imag2, Nl, Nc, &Wavelet_2,
                            Type_Transform, Fc, Nbr_Plan);

    Nbr_Etap = Nbr_Plan - 1;

    for (i = 0; i < Nbr_Etap; i++)
    {
       wavelet_pointer_plan (&Wavelet_1, &Ptr1, &Nl, &Nc, i+1, 0);
       wavelet_pointer_plan (&Wavelet_2, &Ptr2, &Nl, &Nc, i+1, 0);

       compare_plan (Ptr1, Ptr2, Nl, Nc, &(Tab_Snr[i+1]), &(Tab_Correl[i+1]));
    }
    wave_io_free (&Wavelet_1);
    wave_io_free (&Wavelet_2);
}

/****************************************************************************/

int test_tab_exit (File_Name_Step)
char *File_Name_Step;
{
   FILE *File_Des;
   string File_Name;
   int L, Nbr, Val_Ok = 0;
   char *Ptr;

    /* test if ".tbl" is in the file name step */
    strcpy (File_Name, File_Name_Step);
    L = strlen (File_Name_Step);
    Ptr = File_Name_Step;
    while ((L > 4) &&
      ((Ptr[0] != '.') || (Ptr[1] != 't') || (Ptr[2] != 'b') || (Ptr[3] != 'l')))
    {
        Ptr ++;
        L --;
    }

    /* if not, add it to the name */
    if ((Ptr[0] != '.') || (Ptr[1] != 't') || (Ptr[2] != 'b') || (Ptr[3] != 'l'))
    {
        strcat(File_Name, ".tbl");
    }

    /* test if the file exists */
    File_Des = fopen (File_Name, "r");
    if (File_Des == NULL) Val_Ok = -1;
    fclose (File_Des);

    /* return the result */
    return (Val_Ok);
}

/****************************************************************************/

main()
{
    string File_Name_Imag_1, File_Name_Imag_2;
    string File_Tab_1, File_Tab_2;
    int Nbr_Plan,i;
    float *Imag_1, *Imag_2;
    int Nl, Nc, Nl_Tab, Nc_Tab;
    float *Tab_Correl, *Tab_Snr, Val;
    int Unit;
    int Null, Actvals, Buffer_Int, Maxvals, Felem;
    int Stat, Dummy, Des_Tab_1, Des_Tab_2;
    double Exp;
    int Min,temp;
    int nrow, ncol, ocol;

    /* Initialisation */  
    SCSPRO("wave_2d_comp_ima");

    /* read the image 1 file name */
    Felem = 1;
    Maxvals = 60;
    Stat = SCKGETC("IN_A", Felem, Maxvals, &Actvals, File_Name_Imag_1);

    /* read the image 2 file name */
    Felem = 1;
    Maxvals = 60;
    Stat = SCKGETC("IN_B", Felem, Maxvals, &Actvals, File_Name_Imag_2);

    /* read the  number of scales */
    Felem = 1;
    Maxvals = 1;
    Stat = SCKRDI("INPUTI",Felem, Maxvals, &Actvals, &Buffer_Int, &Unit, &Null);
    Nbr_Plan = Buffer_Int + 1;

    /* read the N_Sigma_Parameter */
    Felem = 1;
    Maxvals = 1;
    Stat = SCKRDR("INPUTR", Felem, Maxvals, &Actvals, &N_Sigma, &Unit, &Null);


    /* read the table 1 file name */
    Felem = 1;
    Maxvals = 60;
    Stat = SCKGETC("OUT_A", Felem, Maxvals, &Actvals, File_Tab_1);

    /* read the table 2 file name */
    Felem = 1;
    Maxvals = 60;
    Stat = SCKGETC("OUT_B", Felem, Maxvals, &Actvals, File_Tab_2);

    Tab_Snr = f_vector_alloc (Nbr_Plan);
    Tab_Correl = f_vector_alloc (Nbr_Plan);

    /* Lecture des images d'entree */
    io_read_file_to_pict_f (File_Name_Imag_1, &Imag_1, &Nl, &Nc);
    io_read_file_to_pict_f (File_Name_Imag_2, &Imag_2, &Nl_Tab, &Nc_Tab);

    if ((Nl != Nl_Tab) || (Nc != Nc_Tab))
       SCETER(10, "frames 1 and 2 must have the same size");

    /* test if the number of planes is not too high */
    Min = (Nl < Nc) ? Nl : Nc;
    Exp = (double) Nbr_Plan + 2.;
    temp = pow(2., Exp) + 0.5;
    if (Min < temp) io_err_message_exit (ERR_NUMBER_OF_PLANES, " ");

    /* Comparaison des images */
    compare_wavelet (Imag_1, Imag_2, Nl, Nc, Nbr_Plan, Tab_Snr, Tab_Correl);

#if VISU_PARAM
    {
    char Send[100];

    sprintf(Send," File_Name_Imag1 = %s\n", File_Name_Imag_1);
    SCTPUT(Send);
    sprintf(Send," File_Name_Imag2 = %s\n", File_Name_Imag_2);
    SCTPUT(Send);
    sprintf(Send," File_Name_Tab1 = %s\n", File_Tab_1);
    SCTPUT(Send);
    sprintf(Send," File_Name_Tab2 = %s\n", File_Tab_2);
    SCTPUT(Send);
    sprintf(Send," Nbr_Plan = %d\n", Nbr_Plan);
    SCTPUT(Send);
    for (i = 0;  i < Nbr_Plan; i++)
    {
        sprintf(Send,"%d: correl = %f, snr = %f\n",i,Tab_Correl[i],Tab_Snr[i]);
    SCTPUT(Send);
    }
  }
#endif

    /* the results are stored in Midas table */
    
    Stat = test_tab_exit (File_Tab_1);
    if (Stat == 0)
    {
        /* open a table */
        Stat = TCTOPN (File_Tab_1, F_IO_MODE, &Des_Tab_1);

        /* read information */
        Stat = TCIGET (Des_Tab_1, &ncol, &nrow, &Dummy,&Dummy,&Dummy);
        Nc_Tab = ncol;
        Nl_Tab = nrow;
    }
    else
    {
        nrow = TAB_MAX_SCALE;
        ncol = 2;
        
        /* create a table */
        Stat = TCTINI(File_Tab_1, F_TRANS, F_O_MODE, 
                      ncol, nrow, &Des_Tab_1);   
        Nc_Tab = 1;
        Nl_Tab = nrow;
        
        /* init a column for the scale */
        Stat = TCCINI(Des_Tab_1, D_R4_FORMAT, 1, "F12.4", " ", "Scale", &ocol);

        /* write the data */
        for (nrow = 1; nrow <= Nbr_Plan; nrow++)
        {
           Val = (float) nrow - 1.;
           ncol = 1;
           TCEWRR (Des_Tab_1, nrow, ncol, &Val);
        }
     }
     ncol = Nc_Tab + 1;
     /* init a column for the correlation */
     Stat = TCCINI(Des_Tab_1, D_R4_FORMAT, 1, "F12.4", 
            " ", File_Name_Imag_2, &ocol);   
     Nc_Tab = ncol;
     
     /* write the data */
     for (nrow = 1; nrow <= Nbr_Plan; nrow++)
     {
         i = nrow - 1;
         TCEWRR (Des_Tab_1, nrow, ncol, &(Tab_Correl[i]));
     }


    /* Signal to noise ration */

    Stat = test_tab_exit (File_Tab_2);
    if (Stat == 0)
    {
        Stat = TCTOPN (File_Tab_2, F_IO_MODE, &Des_Tab_2);
        Stat = TCIGET (Des_Tab_2, &ncol, &nrow, &Dummy,&Dummy,&Dummy);
        Nc_Tab = ncol;
        Nl_Tab = nrow;
    }
    else
    {
        nrow = TAB_MAX_SCALE;
        ncol = 2;
        
        Stat = TCTINI(File_Tab_2, F_TRANS, F_O_MODE, 
                      ncol, nrow, &Des_Tab_2);   
        Nc_Tab = 1;
        Nl_Tab = nrow;

        Stat = TCCINI(Des_Tab_2, D_R4_FORMAT, 1, "F12.4", " ", "Scale", &ocol);
        for (nrow = 1; nrow <= Nbr_Plan; nrow++)
        {
           Val = (float) nrow - 1.;
           ncol = 1;
           TCEWRR (Des_Tab_2, nrow, ncol, &Val);
        }
     }
     
     ncol = Nc_Tab + 1;
     Stat = TCCINI(Des_Tab_2, D_R4_FORMAT, 1, "F12.4", 
            " ", File_Name_Imag_2, &ocol);   
     Nc_Tab = ncol;
     
     for (nrow = 1; nrow <= Nbr_Plan; nrow++)
     {
         i = nrow - 1;
         TCEWRR (Des_Tab_2, nrow, ncol, &(Tab_Snr[i]));
     }

    free ((char *) Imag_1);
    free ((char *) Imag_2);
    free ((char *) Tab_Snr);
    free ((char *) Tab_Correl);

    TCTCLO(Des_Tab_1);
    TCTCLO(Des_Tab_2);

    SCSEPI();
}

/***************************************************************************/