/* @(#)wa_pers.c	17.1.1.1 (ES0-DMD) 01/25/02 17:21:41 */
/*===========================================================================
  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) by European Southern Observatory
*******************************************************************************
**
**    UNIT
**
**    Version: 17.1.1.1
**
**    Author: Jean-Luc Starck 
**
**    Date:  02/01/25
**    
**    File: wa_pers.c
**
*******************************************************************************
**
**    DESCRIPTION creation of an image from a wavelet transform
**    ----------- 
**                 
**    PARAMETRES    
**    ----------   
** 
**      File_Name_Transform (keyword: IN_A):
**             File name of the input wavelet transform
**
**      File_Name_Imag (keyword: OUT_A):
**             File name of the output image
**
**      Increment (keyword: INPUTI(1)): (default is 1)
**             we keep one line on Increment
**
**      Visualisation Mode (keyword: IN_B): CO or BW (default is CO)
**               CO =  Color 
**               BW = black and white
**
**      Threshold value (keyword: INPUTR(1)): (default is 5)
**               all the wavelet coefficients > Treshold*Sigma_Scale
**               are set to Treshold*Sigma_Scale
**
**    RESULTS : create an image from a wavelet file   
**    -------   the image is a representation of the 
**              wavelet coefficicients with a 3 dimension 
**              visualisation 
******************************************************************************/

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

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

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

float lib_mat_ecart_type();


#define MAXI 255

#define COULEUR_1 1
#define COULEUR_2 1
#define BLACK_WHITE 0

int MAXX = 512;
int MAXY = 512;

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

ligne (x1,y1,x2,y2,Couleur,Nl,Nc,Tab_Caval)
/* 
   Draw a line in the image Tab_Caval between the two points:
    I1,J1 and I2,J2. The intensity of the line is defined by Couleur
*/
   
int x1,y1,x2,y2,Nl,Nc;
float Couleur;
float *Tab_Caval;
{
    int j,pasx,x,y,dx,dy,cumul;
    float pasy;
    
    if (x1 < x2) pasx = 1;
    else pasx = -1;
    if (y1 < y2) pasy = 1;
    else pasy = -1;
    dx = abs (x2 - x1);
    dy = abs (y2 - y1);
    x= x1;
    y = y1;
    Tab_Caval [x1*Nc+y1] = Couleur;
    
    if (dx > dy)
    {
        cumul = dx / 2;
        for (j = 0; j < dx; j++)
        {
            x += pasx;
            cumul += dy;
            if (cumul > dx)
            {
                cumul -= dx;
                y += pasy;
            }
            Tab_Caval [x*Nc+y] = Couleur;
        }
    }
    else
    {
        
        cumul = dy / 2;
        for (j = 0; j < dy; j++)
        {
            y += pasy;
            cumul += dx;
            if (cumul > dy)
            {
                cumul -= dy;
                x += pasx;
            }
            Tab_Caval [x*Nc+y] = Couleur;
        } 
    }
}  

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

lib_caval_pict (Nl, Nc, Maxy, Maxx, Tab_Img, Tab_Caval, Inc, Mod_Visu)
int Nl,Nc,Maxy,Maxx,Inc;
float *Tab_Img,*Tab_Caval; 
int Mod_Visu;
{
    int j,Size,Indi,Indj;
    float i,Seuil,Inclinaison,X_Pas,Increment,Nb_Line;
    float *Tab_Max;
    int Borne;
    int X1,Y1,X2,Y2;
    int Couleur;
    float Val,H_Facteur,L_Facteur,u,Interp;
    int Angle;
    float Level;

    Size = Maxx * Maxy;    
    Tab_Max = f_vector_alloc (Maxx);
    for (j = 0; j < Maxx; j++) Tab_Max [j] = 0.;
    for (j = 0; j < Size; j++) Tab_Caval [j] = 0.;

    Increment  = (float) Inc;
    {
        float Min, Max;
 
        Min = Max = Tab_Img[0];
        for (j = 0; j < Nl*Nc; j++)
        {
            if (Min > Tab_Img[j]) Min = Tab_Img[j];
            if (Max < Tab_Img[j]) Max = Tab_Img[j];
        }
        for (j = 0; j < Nl*Nc; j++) Tab_Img[j] = (Tab_Img[j]-Min) / (Max - Min);
    }

    /* Color threshold  */
    Seuil = 1. / 5.;

    /* scale parameter realted to the visualisation angle*/
    H_Facteur = (float) Maxy / 2.;
    L_Facteur = (float) Maxx / (float) Nc;
    X_Pas = 1. / L_Facteur;
    Nb_Line = (float) Nl / ((float) Maxy / 2.);
    Borne = (Maxy /2. > Nl)? Nl : Maxy / 2.;
    
    for (i = 0; i < Borne; i += Increment)
    {
        /* First point position */
        u = i * Nb_Line;
        Indi = (int) u * Nc;

        if (Indi > Nl*Nc)
        {
             printf ("(%f,%d)",Nb_Line,Indi);
        }
        Val = H_Facteur * Tab_Img [Indi] +  i;
        if (Val > Tab_Max [0]) Tab_Max [0] = Val;
        X1 = 0;
        Y1 = Val;
        u = 0.;
        
        for (j = 1; j < Maxx; j ++)
        {
            /* ordinate position */
            u += X_Pas;
            u = ((float) Nc / (float) Maxx) * (float) j;
            Indj = u;
            
            if ((L_Facteur > 1.) && ((Indj + 1) < Nc))
            {
                Interp = u - Indj;
                Val  = Tab_Img [Indi+Indj];
                Val = H_Facteur * (Val + (Tab_Img [Indi+Indj+1] 
                                             - Val) * Interp) +  i;
                if ((Val >= Maxy) || (Val < 0.))
                {
                    printf ("zoom val (%d,%d)=%f ",Indi,Indj,Val);
                    exit (0);
                }
            }
            else
            {
                Val = H_Facteur * Tab_Img [Indi+Indj] + i;
                if ((Val >= Maxy) || (Val < 0.))
                {
                    printf ("calcul val %5.1f:(%d,%d)=%f ",i,Indi,Indj,Val);
                    exit (0);
                }
            }
                
            /* Tab_Max = Z-array: contains the max for each columns   */
            if (Val > Tab_Max [j])
            {
                Tab_Max [j] = Val;
                X2 = j;
               /* Y2 = Maxy - 1 - (Tab_Max [j] + 0.5);*/
                Y2 = Val;
                if (Y2 < 0)
		{
                    printf ("Y2 < 0 : %d",Y2);
                    exit (0);
                    Y2 = 0;
		}
		if (Tab_Img [Indi+Indj] < Seuil) Couleur = COULEUR_1;
                else Couleur = COULEUR_2;
                
                if ((Couleur == COULEUR_1) 
                            && ((Tab_Caval [Y1*Maxx+X1] == COULEUR_2)
                                || (Tab_Caval [Y2*Maxx+X2] == COULEUR_2)))
                {
                    X1 = X1+1;
                    Y1 = Y2;
                }
                else
                {
                    Level = Tab_Img [Indi+Indj];
                    if (Y1 < 0)
		    {
                        printf ("Y1(%d,%5.1f)",j,Tab_Max [j]);
                        exit (-1);
                        Y1 = 0;
		    }
		    if (Mod_Visu == BLACK_WHITE)
		    {
		        Level = Couleur;
		    }
                    ligne (Y1,X1,Y2,X2,Level,Maxy,Maxx,Tab_Caval);
                    X1 = X2;
                    Y1 = Y2;
                }
            }
            else
            {
                X1 = j;
               /*  Y1 = Maxy - Tab_Max [j];*/
                Y1 = Tab_Max [j];
                if (Y1 < 0)
		{
                    printf ("Y1(%d,%f)",j,Tab_Max [j]);
                    exit (-1);
                    Y1 = 0;
		}
            }
        }
    }
    
    free ((char *)Tab_Max);
}

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

static copy_plan_to_pict(Tab_Caval, Nl1, Nc1, Pict_Caval, Maxy, Maxx, Num_Etap)
float *Tab_Caval,*Pict_Caval;
int Nl1, Nc1, Maxy, Maxx;
int Num_Etap;
{
    int Offsx,Offsy;
    int i,j;
    int ind;
    
    Offsx = (Maxx - Nc1) / 2;
    Offsy = Nl1 * Num_Etap * Maxx;
    
    for (i = 0; i < Nl1; i++)
    {
        for (j = 0; j < Nc1; j++)
        {
            ind = Offsy + Offsx + i * Maxx  + j;
            Pict_Caval [ind] = Tab_Caval [i * Nc1 + j] * MAXI;
        }
    }
}

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

main()
{
    int i,j;
    int Nl,Nc;
    wave_transf_des Wavelet;
    float *ImaSynt, *Imag, Sigma;
    string File_Name_Imag, File_Name_Transform;
    int Actvals, Maxvals, Felem, Unit, Null;
    int Stat, Buffer_Int;
    int Nl1,Nc1,Size;
    int Inc = 1;
    int Mod_Visu = 1;
    float Zoom, *Pict_Caval, K_Sigma;
    string Buffer_Visu;
    int Size_Plan;

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

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

    /* read the image output file name */
    Stat = SCKGETC("OUT_A", Felem, Maxvals, &Actvals, File_Name_Imag);

    /* read the increment */
    Maxvals = 1;
    Stat = SCKRDI("INPUTI",Felem, Maxvals, &Actvals, &Buffer_Int, &Unit, &Null);
    Inc = Buffer_Int;

    /* read the visualisation mode (Black and White (BW) or Color (CO))*/
    Maxvals = 60;
    Stat = SCKGETC("IN_B", Felem, Maxvals, &Actvals, Buffer_Visu);
    if (   ((Buffer_Visu[0] == 'B') || (Buffer_Visu[0] == 'b'))
        && ((Buffer_Visu[1] == 'W') || (Buffer_Visu[1] == 'w'))) Mod_Visu = 0;

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

   /* read the wavelet */
    wave_io_read (File_Name_Transform, &Wavelet);

    Nl = Wavelet.Nbr_Ligne;
    Nc = Wavelet.Nbr_Col;
    Zoom = (float) MAXX / (float) Nc;
    Size_Plan = (float) MAXY / (float) Wavelet.Nbr_Plan;
    Size = MAXX * MAXY;    
    Pict_Caval = f_vector_alloc (Size);
    ImaSynt = f_vector_alloc (Size);

    for (i = 0; i < Size; i++) ImaSynt[i] = 0.;

    for (i = 0; i <Wavelet.Nbr_Plan ; i++)
    {
        wavelet_extract_plan (&Wavelet, &Imag, &Nl, &Nc, i+1);
        Sigma = K_Sigma * lib_mat_ecart_type(Imag,Nl,Nc);

        Nl1 = Size_Plan;
        Nc1 = (((float) Nc * Zoom) > MAXX) ? MAXX : ((float)Nc * Zoom);
        if (i != Wavelet.Nbr_Plan-1)
        for (j = 0; j < Nl*Nc; j++)
        {
            if (Imag [j] > Sigma) Imag [j] = Sigma;
            if (Imag [j] < - Sigma) Imag [j] = - Sigma;
        }
       lib_caval_pict  (Nl, Nc, Nl1, Nc1, Imag, Pict_Caval, Inc, Mod_Visu);
       copy_plan_to_pict(Pict_Caval,Nl1,Nc1,ImaSynt,MAXY,MAXX, i);
       free ((char *) Imag); 
    }

   io_write_pict_f_to_file (File_Name_Imag, ImaSynt, MAXY, MAXX);

   /* End */
   SCSEPI();
}