/* @(#)mallat.c	17.1.1.1 (ES0-DMD) 01/25/02 17:21:22 */
/*===========================================================================
  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: Jerome Bijaoui - Jean-Luc Starck
**
**    Date:  02/01/25
**    
**    File:  mallat.c
**
*******************************************************************************
**    DESCRIPTION  routines used for the wavelet transform algorithm 
**    -----------  defined by S. Mallat with biorthogal filters (Daubechies)
**
*******************************************************************************
** int ondelette_1d (nb_donnees, signal_entree, signal_sortie, detail_sortie)
** int nb_donnees;
**
** float   *signal_entree, *signal_sortie, *detail_sortie;
**
** one dimension wavelet transform 
**
*******************************************************************************
** int ondelette_inverse_1d (nb_donnees, signal_entree, 
**                          detail_entree, signal_sortie)
** int nb_donnees;
** float   *signal_entree, *detail_entree, *signal_sortie;
**
** one dimension inverse wavelet transform  
**
*******************************************************************************
** int ondelette_2d (n, m, niveau, image_entree, image_h0_h0, 
**                   image_h0_g0, image_g0_h0, image_g0_g0)
** int n, m, niveau;
** float *image_entree, *image_h0_h0, *image_h0_g0;
** float *image_g0_h0, *image_g0_g0;
**
** two dimensions wavelet transform
** 
*******************************************************************************
** int ondelette_inverse_2d (n, m, niveau, image_h1_h1, image_h1_g1, 
**                           image_g1_h1, image_g1_g1, image_sortie)
** int n, m, niveau;
** float *image_h1_h1, *image_h1_g1, *image_g1_h1;
** float *image_g1_g1, *image_sortie;
**
** two dimensions inverse wavelet transform 
**
*******************************************************************************
** mallat_2d_transform (image, Des_Wave, nb_lignes, nb_colonnes, Nbr_Plan)
** float *image;
** struct mallat_plan_des *Des_Wave;
** int nb_lignes, nb_colonnes, Nbr_Plan;
**
** computes the wavelet transform by Mallat's algorithm
** image = input image
** nb_lignes, nb_colonnes = number of lines and columns
** Nbr_Plan = number of scales
** Des_Wave = output: wavelet transform
**
*******************************************************************************
** mallat_2d_reconstruct (image, Des_Wave, nb_lignes, nb_colonnes, Nbr_Plan)
** float *image;
** struct mallat_plan_des *Des_Wave;
** int nb_lignes, nb_colonnes, Nbr_Plan;
**
** Reconstructs an image from its wavelet transform
** Des_Wave = wavelet transform
** nb_lignes, nb_colonnes = number of lines and columns
** Nbr_Plan = number of scales
** image = output: reconstructed image
**
****************************************************************************
** mallat_2d_extract_plan (Imag, Nl, Nc, Imag_Hor, Imag_Diag, Imag_Vert, 
**                         Des_Wave, Nbr_Plan)
** float *Imag, *Imag_Hor, *Imag_Diag, *Imag_Vert;
** struct mallat_plan_des *Des_Wave;
** int Nl, Nc, Nbr_Plan;
**
** Extract a plane from the wavelet transform
** Nl, Nc = number of lines and columns
** Imag_Hor, Imag_Diag, Imag_Vert = details images
** Des_Wave = wavelet transform
** Nbr_Plan = number of scales
** Imag = global representation at this scale
**
**
****************************************************************************
** 
** mallat_2d_enter_plan (Imag, Nl, Nc, Des_Wave, Nbr_Plan)
** float *Imag;
** struct mallat_plan_des *Des_Wave;
** int Nl, Nc, Nbr_Plan;
**
** enter the image in the wavelet transform
** 
** Nl, Nc = number of lines and columns
** Imag_Hor, Imag_Diag, Imag_Vert = details images
** Des_Wave = wavelet transform
** Nbr_Plan = number of scales
** Imag = global representation at this scale
**
***************************************************************************/ 

static char sccsid[] = "@(#)mallat.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"
#include "Def_Mallat.h"

extern float h0[9];
extern float g0[7];
extern float h1[7];
extern float g1[9];

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

int filtrer_h0 (nb_donnees, entree, sortie)
int nb_donnees;
float *entree, *sortie;
/* convolves the data with the filter h0 */
{
    int indice, indice_2, Index, position_masque;

    for (indice = 0; indice < nb_donnees; indice += 2)
    {
        indice_2 = indice>>1;
        sortie[indice_2] = 0;
        for (position_masque = 0; position_masque < 9; position_masque++)
	{
            Index = indice+position_masque-4;
            if (Index < 0 ) Index = -Index;
            if (Index >= nb_donnees) Index = ( (nb_donnees-1)<<1 ) - Index;
            sortie[indice_2] += entree[Index] * h0[position_masque];
	}
    }
    return (0);
}

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

int filtrer_g0 (nb_donnees, entree, sortie)
int nb_donnees;
float   *entree, *sortie;
/* convolves the data with the filter g0 */
{
    int indice, indice_2, Index, position_masque;

    for (indice = 1; indice < nb_donnees; indice += 2)
    {
        indice_2 = (indice-1) >> 1;
        sortie[indice_2] = 0;
        for (position_masque = 0; position_masque<7; position_masque++)
	{
            Index = indice+position_masque-3;
            if (Index < 0) Index = -Index;
            if (Index >= nb_donnees) Index = ((nb_donnees-1)<<1) - Index;
            sortie[indice_2] += entree[Index]*g0[position_masque];
	}
    }
    return (0);
}

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

int filtrer_h1 (nb_donnees, entree, sortie)
int nb_donnees;
float   *entree, *sortie;
/* convolves the data with the filter h1 */
{
    int indice, Index, position_masque;
    float   *temporaire;

    temporaire = (float *) calloc ((unsigned)nb_donnees, sizeof (float));

    for (indice = 0; indice < nb_donnees; indice += 2) 
                             temporaire[indice] = entree[indice>>1];

    for (indice = 0; indice < nb_donnees; indice++)
    {
        sortie[indice] = 0;
        for (position_masque = 0; position_masque < 7; position_masque++)
	{
            Index = indice+position_masque-3;
            if (Index< 0) Index = -Index;
            if (Index>=nb_donnees) Index = ((nb_donnees-1)<<1)-Index;
            sortie[indice] += temporaire[Index]*h1[position_masque];
	}
    }

    free ((char *) temporaire);
    return (0);
}

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

int filtrer_g1 (nb_donnees, entree, sortie)
int nb_donnees;
float   *entree, *sortie;
/* convolves the data with the filter g1 */
{
    int indice, Index, position_masque;
    float   *temporaire;

    temporaire = (float *) calloc ((unsigned)nb_donnees, sizeof (float));

    for (indice = 1; indice<nb_donnees; indice += 2) 
                                 temporaire[indice] = entree[(indice-1)>>1];
	
    for (indice = 0; indice<nb_donnees; indice++)
    {
        sortie[indice] = 0;
        for (position_masque = 0; position_masque < 9; position_masque++)
	{
            Index = indice+position_masque-4;
            if (Index < 0) Index = -Index;
            if (Index>=nb_donnees) Index = ((nb_donnees-1)<<1)-Index;
            sortie[indice] += temporaire[Index]*g1[position_masque];
	}
    }
    free ((char *) temporaire);
    return (0);
}

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

int ondelette_inverse_1d (nb_donnees, signal_entree, 
                          detail_entree, signal_sortie)
int nb_donnees;
float   *signal_entree, *detail_entree, *signal_sortie;
{
    int Index;
    float   *temporaire;

    temporaire = (float *)
              calloc ((unsigned)nb_donnees, (unsigned)sizeof (float));

    filtrer_h1 (nb_donnees, signal_entree, signal_sortie);
    filtrer_g1 (nb_donnees, detail_entree, temporaire);

    for (Index = 0; Index<nb_donnees; Index++) 
            signal_sortie[Index] = 2.*(signal_sortie[Index]+temporaire[Index]);
    free ((char *) temporaire);
    return (0);
}

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

int ondelette_1d (nb_donnees, signal_entree, signal_sortie, detail_sortie)
int nb_donnees;
float   *signal_entree, *signal_sortie, *detail_sortie;

{
    filtrer_h0 (nb_donnees, signal_entree, signal_sortie);
    filtrer_g0 (nb_donnees, signal_entree, detail_sortie);
    return (0);
}

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

int ondelette_2d (n, m, niveau, image_entree, image_h0_h0, 
                  image_h0_g0, image_g0_h0, image_g0_g0)
int n, m, niveau;
float   *image_entree, *image_h0_h0, *image_h0_g0, *image_g0_h0, *image_g0_g0;
{
    int nb_lignes, nb_colonnes, nb_colonnes_2, nb_lignes_2;
    int colonne, ligne, Index;
    float *image_h0, *image_g0, *colonne_h0, *colonne_g0, *colonne_h0_h0;
    float *colonne_h0_g0, *colonne_g0_h0,*colonne_g0_g0;

    nb_lignes = m >> niveau;
    nb_colonnes   = n          >>niveau;
    nb_colonnes_2 = nb_colonnes>>     1;
    nb_lignes_2   = nb_lignes  >>     1;

    image_h0 = f_vector_alloc (nb_lignes*nb_colonnes_2);
    image_g0 = f_vector_alloc (nb_lignes*nb_colonnes_2);
    for (ligne = 0; ligne<nb_lignes; ligne++) 
         ondelette_1d (nb_colonnes, 
                      image_entree+nb_colonnes*ligne,
                      image_h0+nb_colonnes_2*ligne,
                      image_g0+nb_colonnes_2*ligne);

    colonne_h0    = f_vector_alloc (nb_lignes);
    colonne_g0    = f_vector_alloc (nb_lignes);
    colonne_h0_h0 = f_vector_alloc (nb_lignes_2);
    colonne_h0_g0 = f_vector_alloc (nb_lignes_2);
    colonne_g0_h0 = f_vector_alloc (nb_lignes_2);
    colonne_g0_g0 = f_vector_alloc (nb_lignes_2);

    for (colonne = 0; colonne<nb_colonnes_2; colonne++)
    {
        for (ligne = 0; ligne<nb_lignes; ligne++)
	{
            Index = nb_colonnes_2 * ligne + colonne;
            colonne_h0[ligne] = image_h0[Index];
            colonne_g0[ligne] = image_g0[Index];
	}
        ondelette_1d (nb_lignes, colonne_h0, colonne_h0_h0, colonne_h0_g0);
        ondelette_1d (nb_lignes, colonne_g0, colonne_g0_h0, colonne_g0_g0);
        for (ligne = 0; ligne<nb_lignes_2; ligne++)
	{
            Index = nb_colonnes_2 * ligne + colonne;
            image_h0_h0[Index] = colonne_h0_h0[ligne];
            image_h0_g0[Index] = colonne_h0_g0[ligne];
            image_g0_h0[Index] = colonne_g0_h0[ligne];
            image_g0_g0[Index] = colonne_g0_g0[ligne];
        }
    }
    free ((char *) image_h0);
    free ((char *) image_g0);
    free ((char *) colonne_h0);
    free ((char *) colonne_g0);
    free ((char *) colonne_h0_h0);
    free ((char *) colonne_h0_g0);
    free ((char *) colonne_g0_h0);
    free ((char *) colonne_g0_g0);
    return (0);
}			

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

int ondelette_inverse_2d (n, m, niveau, image_h1_h1, image_h1_g1, 
                          image_g1_h1, image_g1_g1, image_sortie)
int n, m, niveau;
float *image_h1_h1, *image_h1_g1, *image_g1_h1;
float *image_g1_g1, *image_sortie;
{
    int nb_lignes, nb_colonnes, nb_colonnes_2, nb_lignes_2;
    register int colonne, ligne, Index;
    float *image_h1, *image_g1, *colonne_h1_h1;
    float *colonne_h1_g1, *colonne_g1_h1, *colonne_g1_g1;
    float *colonne_h1, *colonne_g1;

    nb_lignes = m >> niveau;
    nb_colonnes = n >> niveau;
    nb_colonnes_2 = nb_colonnes >> 1;
    nb_lignes_2 = nb_lignes >> 1;

    /* Allocation  */
    image_h1 = f_vector_alloc (nb_lignes*nb_colonnes_2);
    image_g1 = f_vector_alloc (nb_lignes*nb_colonnes_2);

    colonne_h1 = f_vector_alloc (nb_lignes);
    colonne_g1 = f_vector_alloc (nb_lignes);
    colonne_h1_h1 = f_vector_alloc (nb_lignes_2);
    colonne_h1_g1 = f_vector_alloc(nb_lignes_2);
    colonne_g1_h1 = f_vector_alloc (nb_lignes_2);
    colonne_g1_g1 = f_vector_alloc(nb_lignes_2);

    /* Transform on the columns */
    for (colonne = 0; colonne<nb_colonnes_2; colonne++)
    {
        for (ligne = 0; ligne<nb_lignes_2; ligne++)
	{
            Index = nb_colonnes_2*ligne+colonne;
            colonne_h1_h1[ligne] = image_h1_h1[Index];
            colonne_h1_g1[ligne] = image_h1_g1[Index];
            colonne_g1_h1[ligne] = image_g1_h1[Index];
            colonne_g1_g1[ligne] = image_g1_g1[Index];
	}
        ondelette_inverse_1d (nb_lignes, colonne_h1_h1, 
                              colonne_h1_g1, colonne_h1);
        ondelette_inverse_1d (nb_lignes, colonne_g1_h1, 
                              colonne_g1_g1, colonne_g1);
        for (ligne = 0; ligne<nb_lignes; ligne++)
	{
            Index = nb_colonnes_2 * ligne + colonne;
            image_h1[Index] = colonne_h1[ligne];
            image_g1[Index] = colonne_g1[ligne];
	}
    }

    free ((char *) colonne_h1);
    free ((char *) colonne_g1);
    free ((char *) colonne_h1_h1);
    free ((char *) colonne_h1_g1);
    free ((char *) colonne_g1_h1);
    free ((char *) colonne_g1_g1);

    /* Transforms on the lines */
    for (ligne = 0; ligne<nb_lignes; ligne++) 
            ondelette_inverse_1d (nb_colonnes, 
                                  image_h1+nb_colonnes_2*ligne,
                                  image_g1+nb_colonnes_2*ligne,
                                  image_sortie+nb_colonnes*ligne);

    free ((char *) image_h1);
    free ((char *) image_g1);
    return (0);
}

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

mallat_2d_transform (image, Des_Wave, nb_lignes, nb_colonnes, Nbr_Plan)
float *image;
struct mallat_plan_des *Des_Wave;
int nb_lignes, nb_colonnes, Nbr_Plan;
{
    int i, profondeur;
    int dim_images_resultat, niveau, nb_lignes_2, nb_colonnes_2;
    float *image_resultat, *image_resultat_2;
    float *image_resultat_3, *image_resultat_4;
    register int Index;
    struct mallat_plan_des *Ptr_Mallat;
 
   /* Allocation  */
    image_resultat      = f_vector_alloc ((nb_colonnes>>1) * (nb_lignes>>1));
    image_resultat_2    = f_vector_alloc ((nb_colonnes>>1) * (nb_lignes>>1));
    image_resultat_3    = f_vector_alloc ((nb_colonnes>>1) * (nb_lignes>>1));
    image_resultat_4    = f_vector_alloc ((nb_colonnes>>1) * (nb_lignes>>1));

    nb_colonnes_2 = nb_colonnes;
    nb_lignes_2 = nb_lignes;
    profondeur = Nbr_Plan - 1;

    Ptr_Mallat = Des_Wave;

    /* Calcule the wavelet coefficients */
    for (niveau = 0; niveau < profondeur; niveau++)
    {
        nb_colonnes_2 = nb_colonnes_2>>1;
        nb_lignes_2 = nb_lignes_2  >>1;
        dim_images_resultat = nb_colonnes_2 * nb_lignes_2;
        ondelette_2d (nb_colonnes, nb_lignes, niveau, image, 
                      image_resultat, 
                      image_resultat_2, 
                      image_resultat_3, 
                      image_resultat_4);


	/* We keep the details */
	Ptr_Mallat -> Nl = nb_lignes_2;
	Ptr_Mallat -> Nc = nb_colonnes_2;

	for (i=0; i< dim_images_resultat; i++)
	{ 
	    (Ptr_Mallat -> Coef_Vert) [i] = image_resultat_2[i];
	    (Ptr_Mallat -> Coef_Horiz) [i]  = image_resultat_3[i];
	    (Ptr_Mallat -> Coef_Diag) [i]  = image_resultat_4[i];
	}

        for (Index = 0; Index < dim_images_resultat; Index++) 
                                image[Index] = image_resultat[Index];

        if (niveau < profondeur - 1) Ptr_Mallat = Ptr_Mallat -> Smooth_Imag;
    }

    for (i=0; i< dim_images_resultat; i++) 
                          (Ptr_Mallat -> Low_Resol)[i] = image_resultat[i];

    free ((char *)image_resultat  );
    free ((char *)image_resultat_2);
    free ((char *)image_resultat_3);
    free ((char *)image_resultat_4);
}

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

mallat_2d_reconstruct (image, Des_Wave, nb_lignes, nb_colonnes, Nbr_Plan)
float *image;
struct mallat_plan_des *Des_Wave;
int nb_lignes, nb_colonnes, Nbr_Plan;
{
    int profondeur, dim_images_entree, niveau;
    register int Index,i;
    float *image_h1_h1;
    float *image_h1_g1, *image_g1_h1, *image_g1_g1;
    struct mallat_plan_des *Ptr_Mallat;

    profondeur = Nbr_Plan - 1;

    /* Allocation */
    image_h1_h1 = f_vector_alloc (nb_colonnes*nb_lignes);

    /* input outut array size definition */
    dim_images_entree = nb_colonnes*nb_lignes>>(profondeur<<1);

    /* Position on the last plan  */
    Ptr_Mallat = Des_Wave;
    for (i = 1; i < profondeur; i++) Ptr_Mallat = Ptr_Mallat -> Smooth_Imag;

     /* initial image constructuction : image_h1_h1. */
    for (Index = 0; Index < (Ptr_Mallat -> Nl * Ptr_Mallat -> Nc); Index++) 
                     image_h1_h1[Index] = (Ptr_Mallat -> Low_Resol)[Index];
    
    for (niveau = profondeur; niveau > 0; niveau--)
    {
	Ptr_Mallat = Des_Wave;
	for (i = 1; i < niveau; i++) Ptr_Mallat = Ptr_Mallat -> Smooth_Imag; 

	image_h1_g1 = Ptr_Mallat -> Coef_Vert;
	image_g1_h1 = Ptr_Mallat -> Coef_Horiz;
	image_g1_g1 = Ptr_Mallat -> Coef_Diag;

        ondelette_inverse_2d (nb_colonnes, nb_lignes, niveau-1, 
                              image_h1_h1, image_h1_g1, image_g1_h1,
                              image_g1_g1, image);

        /* New dimensions */
        dim_images_entree <<= 2;

        /* Next iteration */
        for (Index = 0; Index < dim_images_entree; Index++) 
                                      image_h1_h1[Index] = image[Index];
    }

    free((char *) image_h1_h1);
}

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

mallat_2d_extract_plan (Imag, Nl, Nc, Imag_Hor, Imag_Diag, Imag_Vert, 
                        Des_Wave, Nbr_Plan)
float *Imag, *Imag_Hor, *Imag_Diag, *Imag_Vert;
struct mallat_plan_des *Des_Wave;
int Nl, Nc, Nbr_Plan;
{
    int Nl_2, Nc_2, i, j, ind, Nbr_Etap = Nbr_Plan - 1;
    struct mallat_plan_des *Ptr_Mallat;
    int Dep_H_Nl,Dep_H_Nc, Dep_V_Nl, Dep_V_Nc;
    int Dep_D_Nl, Dep_D_Nc, Dep_L_Nl, Dep_L_Nc;
    int Num_Etap;

    Nl_2 = Des_Wave -> Nl;
    Nc_2 = Des_Wave -> Nc;

    /* copy the details */
    for (i = 0; i < Nl_2*Nc_2; i++)
    {
        Imag_Hor[i] = (Des_Wave -> Coef_Horiz) [i];
        Imag_Diag[i] = (Des_Wave -> Coef_Diag) [i];
        Imag_Vert[i] = (Des_Wave -> Coef_Vert) [i];
    }

     /* Position of each detail image in the big image */
    Dep_H_Nl = Nl_2;
    Dep_H_Nc = Nc_2;
    Dep_V_Nl = 0;
    Dep_V_Nc = 0;
    Dep_D_Nl = 0;
    Dep_D_Nc = Nc_2;
    Dep_L_Nl = Nl - Nl_2;
    Dep_L_Nc = 0;

    Ptr_Mallat = Des_Wave;
    for (Num_Etap = 1; Num_Etap <= Nbr_Etap; Num_Etap++)
    {
        /* detail plane size */
        Nl_2 = Ptr_Mallat -> Nl;
        Nc_2 = Ptr_Mallat -> Nc;

        /* copy the details in the big image */
        for (i = 0; i < Nl_2; i++)
        {
            for (j = 0; j < Nc_2; j++)
            {
                ind = (i + Dep_H_Nl) * Nc + Dep_H_Nc + j;
                Imag [ind]  = (Ptr_Mallat -> Coef_Horiz) [i*Nc_2 + j];
                ind = (i + Dep_V_Nl) * Nc + Dep_V_Nc + j;
                Imag [ind]  = (Ptr_Mallat -> Coef_Vert) [i*Nc_2 + j];
                ind = (i + Dep_D_Nl) * Nc + Dep_D_Nc + j;
                Imag [ind]  = (Ptr_Mallat -> Coef_Diag) [i*Nc_2 + j];
                if (Num_Etap == Nbr_Etap)
                {
                    ind = (i + Dep_L_Nl) * Nc + Dep_L_Nc + j;
                    Imag [ind] = (Ptr_Mallat -> Low_Resol) [i*Nc_2 + j];
                }
            }
        }
        Dep_H_Nl += Nl_2 / 2;
        Dep_V_Nl += Nl_2;
        Dep_D_Nl += Nl_2;
        Dep_L_Nl += Nl_2 / 2;

        Dep_H_Nc -= Nc_2 / 2;
        Dep_V_Nc -= 0.;
        Dep_D_Nc -= Nc_2 / 2;
        Dep_L_Nc -= 0;
        /* Next iteration */
        Ptr_Mallat = Ptr_Mallat -> Smooth_Imag;
    }
}

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

mallat_2d_visu (Wavelet, Imag, Nl_Plan, Nc_Plan)
float **Imag;
wave_transf_des *Wavelet;
int *Nl_Plan, *Nc_Plan;
{
    struct mallat_plan_des *Ptr_Mallat;
    int Nbr_Plan;

    Nbr_Plan = Wavelet -> Nbr_Plan;
    Ptr_Mallat = &(Wavelet -> Mallat);

    /* Image size of the image containing this plane + the next */
    *Nl_Plan = Ptr_Mallat -> Nl * 2;
    *Nc_Plan = Ptr_Mallat -> Nc * 2;
    *Imag = f_vector_alloc (*Nl_Plan * *Nc_Plan);

     /* Extraction of the planes */
     mallat_2d_norm (*Imag, *Nl_Plan, *Nc_Plan, Ptr_Mallat, Nbr_Plan);
}

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

mallat_2d_norm (Imag,Nl,Nc,Des_Wave,Nbr_Plan)
float *Imag;
struct mallat_plan_des *Des_Wave;
int Nl, Nc, Nbr_Plan;
{
    int Nl_2, Nc_2, i, j, ind, Nbr_Etap = Nbr_Plan - 1;
    struct mallat_plan_des *Ptr_Mallat;
    int Dep_H_Nl,Dep_H_Nc, Dep_V_Nl, Dep_V_Nc;
    int Dep_D_Nl, Dep_D_Nc, Dep_L_Nl, Dep_L_Nc;
    int Num_Etap;

    Nl_2 = Des_Wave -> Nl;
    Nc_2 = Des_Wave -> Nc;

     /* Position of each detail image in the big image */
    Dep_H_Nl = Nl_2;
    Dep_H_Nc = Nc_2;
    Dep_V_Nl = 0;
    Dep_V_Nc = 0;
    Dep_D_Nl = 0;
    Dep_D_Nc = Nc_2;
    Dep_L_Nl = Nl - Nl_2;
    Dep_L_Nc = 0;

    Ptr_Mallat = Des_Wave;
    for (Num_Etap = 1; Num_Etap <= Nbr_Etap; Num_Etap++)
    {
        /* taille de chaque plan detail */
        Nl_2 = Ptr_Mallat -> Nl;
        Nc_2 = Ptr_Mallat -> Nc;

        NORM_TO_1(Ptr_Mallat -> Coef_Horiz, Nl_2 * Nc_2);
        NORM_TO_1(Ptr_Mallat -> Coef_Vert, Nl_2 * Nc_2);
        NORM_TO_1(Ptr_Mallat -> Coef_Diag, Nl_2 * Nc_2);
        if (Num_Etap == Nbr_Etap)
        {
            NORM_TO_1(Ptr_Mallat -> Low_Resol, Nl_2 * Nc_2);
        }

        /* copy the details in the big image */
        for (i = 0; i < Nl_2; i++)
        {
            for (j = 0; j < Nc_2; j++)
            {
                ind = (i + Dep_H_Nl) * Nc + Dep_H_Nc + j;
                Imag [ind]  = (Ptr_Mallat -> Coef_Horiz) [i*Nc_2 + j];
                ind = (i + Dep_V_Nl) * Nc + Dep_V_Nc + j;
                Imag [ind]  = (Ptr_Mallat -> Coef_Vert) [i*Nc_2 + j];
                ind = (i + Dep_D_Nl) * Nc + Dep_D_Nc + j;
                Imag [ind]  = (Ptr_Mallat -> Coef_Diag) [i*Nc_2 + j];
                if (Num_Etap == Nbr_Etap)
                {
                    ind = (i + Dep_L_Nl) * Nc + Dep_L_Nc + j;
                    Imag [ind] = (Ptr_Mallat -> Low_Resol) [i*Nc_2 + j];
                }
            }
        }
        Dep_H_Nl += Nl_2 / 2;
        Dep_V_Nl += Nl_2;
        Dep_D_Nl += Nl_2;
        Dep_L_Nl += Nl_2 / 2;

        Dep_H_Nc -= Nc_2 / 2;
        Dep_V_Nc -= 0.;
        Dep_D_Nc -= Nc_2 / 2;
        Dep_L_Nc -= 0;

        Nc_2 *= 2;
        Nl_2 *= 2;
        for (i = 0; i < Nl_2; i++)
        {
            Imag [(Nl-i-1)*Nc+Nc_2/2] = 1.;
            Imag [(Nl-Nl_2/2-1)*Nc+i] = 1.;
        }

        /* Next iteration */
        Ptr_Mallat = Ptr_Mallat -> Smooth_Imag;
    }
}

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

mallat_2d_enter_plan (Imag, Nl, Nc, Des_Wave, Nbr_Plan)
float *Imag;
struct mallat_plan_des *Des_Wave;
int Nl, Nc, Nbr_Plan;
{
    int Nl_2, Nc_2, i, j, ind, Nbr_Etap = Nbr_Plan - 1;
    struct mallat_plan_des *Ptr_Mallat;
    int Dep_H_Nl,Dep_H_Nc, Dep_V_Nl, Dep_V_Nc;
    int Dep_D_Nl, Dep_D_Nc, Dep_L_Nl, Dep_L_Nc;
    int Num_Etap;

    Nl_2 = Des_Wave -> Nl;
    Nc_2 = Des_Wave -> Nc;

     /* Position of each detail image in the big image */
    Dep_H_Nl = Nl_2;
    Dep_H_Nc = Nc_2;
    Dep_V_Nl = 0;
    Dep_V_Nc = 0;
    Dep_D_Nl = 0;
    Dep_D_Nc = Nc_2;
    Dep_L_Nl = Nl - Nl_2;
    Dep_L_Nc = 0;

    Ptr_Mallat = Des_Wave;
    for (Num_Etap = 1; Num_Etap <= Nbr_Etap; Num_Etap++)
    {
        /* taille de chaque plan detail */
        Nl_2 = Ptr_Mallat -> Nl;
        Nc_2 = Ptr_Mallat -> Nc;

        /* copy the details in the big image */
        for (i = 0; i < Nl_2; i++)
        {
            for (j = 0; j < Nc_2; j++)
            {
                ind = (i + Dep_H_Nl) * Nc + Dep_H_Nc + j;
                (Ptr_Mallat -> Coef_Horiz) [i*Nc_2 + j] = Imag [ind];
                ind = (i + Dep_V_Nl) * Nc + Dep_V_Nc + j;
                (Ptr_Mallat -> Coef_Vert) [i*Nc_2 + j] = Imag [ind];
                ind = (i + Dep_D_Nl) * Nc + Dep_D_Nc + j;
                (Ptr_Mallat -> Coef_Diag) [i*Nc_2 + j] = Imag [ind];
                if (Num_Etap == Nbr_Etap)
                {
                    ind = (i + Dep_L_Nl) * Nc + Dep_L_Nc + j;
                    (Ptr_Mallat -> Low_Resol) [i*Nc_2 + j] = Imag [ind];
                }
            }
        }
        Dep_H_Nl += Nl_2 / 2;
        Dep_V_Nl += Nl_2;
        Dep_D_Nl += Nl_2;
        Dep_L_Nl += Nl_2 / 2;

        Dep_H_Nc -= Nc_2 / 2;
        Dep_V_Nc -= 0.;
        Dep_D_Nc -= Nc_2 / 2;
        Dep_L_Nc -= 0;

        /* Next iteration */
        Ptr_Mallat = Ptr_Mallat -> Smooth_Imag;
    }
}

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