/* @(#)pyr_cf.c	17.1.1.1 (ES0-DMD) 01/25/02 17:21:23 */
/*===========================================================================
  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:  pyr_cf.c
**
*******************************************************************************
**
**    DESCRIPTION  routines for the pyramidal wavelet transform algorithm
**    -----------  which uses the FFT
**
*******************************************************************************
** 
** pyr_2d_cf_pyr_re (Pyr, Tab_Nl, Tab_Col, Tab_Pos, Nbr_Plan, Pyr_re)
** complex_float *Pyr;
** int *Tab_Nl, *Tab_Col, *Tab_Pos, Nbr_Plan;
** float *Pyr_re
**
** Copies the real part of the pyramid Pyr in Pyr_re
** 
****************************************************************************** 
**
** pyr_2d_cf_fft_2d (Pyr, Tab_Nl, Tab_Col, Tab_Pos, Nbr_Plan, Direction)
** complex_float *Pyr;
** int *Tab_Nl, *Tab_Col, *Tab_Pos, Nbr_Plan;
** int Direction;
**
** Computes the Fourier transform off all the planes
** if Direction = 1 then FFt
** if Direction = -1 then inverse  FFt
**
****************************************************************************** 
**
** pyr_2d_cf_normalize_fft (Pyr, Tab_Nl, Tab_Col, Tab_Pos, Nbr_Plan, Dir)
** complex_float *Pyr;
** int *Tab_Nl, *Tab_Col, *Tab_Pos, Nbr_Plan;
** int Dir;
**
**  Normalize the pyramid planes 
**  if Direction = 1 then normalisation Direct space -> Frequency
**  if Direction = -1 then normalisation Frequency space --> Direct
**
****************************************************************************** 
**
** pyr_2d_cf_mult_tab_imag (Pict, Result, Which_Filter, 
**                          Nl, Nc, Nl1, Nc1, Dep, Fc, Type_Transform)
** complex_float *Pict,*Result;
** int Which_Filter;
** int Nl, Nc, Nl1, Nc1;
** int Dep;
** float Fc;
** 
** Multiplies an image (Nl1xNc1) by the filter defined by Which filter.
** The size of the filter is NlxNc
** Dep is a power of  2
** Fc is the cut off-frequency
** Type_Transform is theselected wavelet transform algorithm
**  Result [i] = Pict [i] * Filter [Dep * i]
**
****************************************************************************** 
**
** pyr_2d_cf_resol_down (Im_Hf, Im_Bf, Tab_Nl, Tab_Col,  
**                       Num_Etap, Fc, Type_Transform)
** complex_float *Im_Hf,*Im_Bf;
** int *Tab_Nl, *Tab_Col;
** int Num_Etap, Type_Transform;
** float Fc;
**
** From the plane Im_Bf, we compute the wavelet coefficients in Im_Hf
** and the image at a lower resolution in Im_Bf. Im_Bf is undersampled
**
****************************************************************************** 
**
** pyr_2d_cf_tfo (Imag, Pyramid, Tab_Nl, Tab_Col, Tab_Pos, 
**                Nbr_Etap, Fc, Type_Transform)
** complex_float *Imag,*Pyramid;
** int *Tab_Nl, *Tab_Col, *Tab_Pos;
** int Nbr_Etap, Type_Transform;
** float Fc;
**
** Wavelet transform of an image. The result is put in a pyramid
** Pyramid will contains the Fourier transform of the wavelet pyramid
** Tab_Nl, Tab_Col, Tab_Pos defined the pyramid
** Nbr_Etap = Nbr_Plan - 1
** Fc = frequence de coupure
**
****************************************************************************** 
**
** pyr_2d_cf_shanon_interpolate (Imag, Imag_up, Nl, Nc, Nl_up, Nc_up)
** omplex_float *Imag, *Imag_up;
** int Nl, Nc, Nl_up, Nc_up;
**
** Interpolation by Schannon of an image (Nl,Nc) to an image (Nl_up,Nc_up)
** Imag and Imag_up are Fourier transform
**
****************************************************************************** 
**
** pyr_2d_cf_extra_wave (Pyramid, Tab_Nl, Tab_Col, Tab_Pos, 
**                       Nbr_Etap, Pyr_Up, Fc, Type_Transform)
** complex_float *Pyramid, *Pyr_Up;
** int *Tab_Nl, *Tab_Col, *Tab_Pos;
** int Nbr_Etap, Type_Transform;
** float Fc;
** 
** Computes the wavelet coefficient of each plane from the precedent plane
** (interpolation of the wavelet plane at a lower scale)
**
****************************************************************************** 
**
** pyr_2d_cf_resol_up (Im_Hf, Im_Bf, Im_Up, Tab_Nl, Tab_Col, 
**                     Num_Etap, Fc, Type_Transform)
** complex_float *Im_Hf,*Im_Bf,*Im_Up;
** int *Tab_Nl, *Tab_Col;
** int Num_Etap, Type_Transform;
** float Fc;
**
** Computes the plane a upper resolution from a plane at a lower resolution
** and from the wavelet plane.
** Im_Hf and Im_Bf are the high and low frequency planes
** Computes the new Im_Hf from Im_Hf and Im_Bf
** All the data are in Fourier space
**
****************************************************************************** 
**
** pyr_2d_cf_build  (Imag, Pyramid, Tab_Nl, Tab_Col, Tab_Pos, 
**                   Nbr_Etap, Fc, Type_Transform)
** complex_float *Pyramid;
** complex_float *Imag;
** int *Tab_Nl, *Tab_Col, *Tab_Pos;
** int Nbr_Etap, Type_Transform;
** float Fc;
** 
** Compute  the Fourier transform of an image  from the Fourier
** transform of its wavelet pyramid
**
****************************************************************************** 
**
** pyr_2d_cf_build_direct (Imag, Pyramid, Tab_Nl, Tab_Col, Tab_Pos, Nbr_Etap)
** complex_float *Pyramid;
** complex_float *Imag;
** int *Tab_Nl, *Tab_Col, *Tab_Pos;
** int Nbr_Etap;
**
** Computes the Fourier transform of an image  from the Fourier
** transform of its wavelet pyramid by a direct addition of the
** wavelet coefficient (case of the difference between
**  two interpolations without changing the wavelet coefficient
**
****************************************************************************** 
**
** pyr_2d_cf_transform (Imag, Pyr, Tab_Nl, Tab_Col, Tab_Pos,
**                      Nbr_Plan,Type_Transform,Fc)
** float *Imag, *Pyr;
** int *Tab_Nl, *Tab_Col, *Tab_Pos;
** int Nbr_Plan, Type_Transform;
** float Fc;
**
** Compute the pyramid in the real space from the image
**
****************************************************************************** 
**
** pyr_2d_cf_build_pict_from_pyr (Pyr, Imag, Tab_Nl, Tab_Col, Tab_Pos, 
**                               Nbr_Plan, Type_Transform, Fc, Build_Direct_Ok)
** float *Imag, *Pyr;
** int *Tab_Nl, *Tab_Col, *Tab_Pos;
** int Nbr_Plan, Type_Transform, Build_Direct_Ok;
** float Fc;
**
** Compute the image in the the real space from the pyramid 
** if Build_Direct_Ok = 1 (TRUE) then the reconstruction is done by a
** sommation of the wavelet coefficients in the Fourier space
** if Build_Direct_Ok = 0 (FALSE) then the reconstruction is done
** in a least mean square sens
**
****************************************************************************** 
**
** pyr_2d_cf_interp_plan (Pyr, Imag, Tab_Nl, Tab_Col, 
**                       Tab_Pos, Num_Plan, Num_Plan_Interp)
** float *Imag, *Pyr;
** int *Tab_Nl, *Tab_Col, *Tab_Pos;
** int Num_Plan, Num_Plan_Interp;
**
** Pyr and Imag are in the real space
** Interpolates the plan  Num_Plan of the pyramid Pyr to the
** resolution defined Num_Plan_Interp
** the result is placed in Imag
** 
*****************************************************************************/ 

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


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

pyr_2d_cf_pyr_re (Pyr, Tab_Nl, Tab_Col, Tab_Pos, Nbr_Plan, Pyr_re)
complex_float *Pyr;
int *Tab_Nl, *Tab_Col, *Tab_Pos, Nbr_Plan;
float *Pyr_re;
{
    int Size,Nbr_Etap,i;

    Nbr_Etap = Nbr_Plan - 1;
    Size = Tab_Pos [Nbr_Etap] + Tab_Nl[Nbr_Etap] * Tab_Col[Nbr_Etap];
    for (i = 0; i < Size; i++) Pyr_re [i] = Pyr [i].re;
}

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

pyr_2d_cf_fft_2d (Pyr, Tab_Nl, Tab_Col, Tab_Pos, Nbr_Plan, Direction)
complex_float *Pyr;
int *Tab_Nl, *Tab_Col, *Tab_Pos, Nbr_Plan;
int Direction;
{
    int i,j,Pos,Nl1,Nc1;
    complex_float *Ptr_Pyr;
    
    for (i = 0; i < Nbr_Plan; i++)
    {
        Nl1 = Tab_Nl [i];
        Nc1 = Tab_Col [i];
        Pos = Tab_Pos [i];
        Ptr_Pyr = Pyr + Pos;

        /* FFT of the plane */
        ft_cf_any_power_of_2 (Ptr_Pyr, Direction, Nl1);

        /* the imaginary part must be equal to 0 in the direct space */
        if (Direction == -1)
        {
            for (j = 0; j < Nl1 * Nc1; j++) Pyr [Pos+j].im = 0.;
        }
    }
}

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

static pyr_2d_cf_normalise_plan (Plan, Nl, Nl_0, Dir)
complex_float *Plan;
int Nl,Nl_0;
int Dir;
{
    int j;
    float Coef;

    /* Computes the normalisation term */
    if (Dir == -1) Coef = (float)(Nl * Nl) / (float)(Nl_0 * Nl_0);
    else Coef = (float)(Nl_0 * Nl_0) / (float)(Nl * Nl);

    for (j = 0; j < Nl * Nl; j++)
    {
        Plan [j].re *= Coef;
        Plan [j].im *= Coef;
    }
}

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


pyr_2d_cf_normalize_fft (Pyr, Tab_Nl, Tab_Col, Tab_Pos, Nbr_Plan, Dir)
complex_float *Pyr;
int *Tab_Nl, *Tab_Col, *Tab_Pos, Nbr_Plan;
int Dir;
{
    int i,Nl,Pos,Nl1;
    complex_float *Ptr_Pyr;

    Nl = Tab_Nl [0];
    for (i = 1; i < Nbr_Plan; i++)
    {
        Nl1 = Tab_Nl [i];
        Pos = Tab_Pos [i];
        Ptr_Pyr = Pyr + Pos;

        pyr_2d_cf_normalise_plan (Ptr_Pyr, Nl1, Nl, Dir);
    }
}

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

static reduction (Pict, Nl, Nc, Nl1, Nc1)
/* reduction of an image in the Fourier space,
   we cut the high frequency */
complex_float *Pict;
int Nl,Nc,Nl1,Nc1;
{
    int i,j,ind, u,v,ind1;
    
    for (i = 0; i < Nl1; i++)
    {
        u = i - Nl1 / 2;
        for (j = 0; j < Nc1; j++)
        {
            ind = i * Nc1 + j;
            v = j - Nc1 / 2;
            ind1 = (u + Nl/2) * Nc + v + Nc / 2;
            Pict [ind].re =  Pict [ind1].re;
            Pict [ind].im =  Pict [ind1].im;
        }
    }
}

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

pyr_2d_cf_mult_tab_imag (Pict, Result, Which_Filter, 
                         Nl, Nc, Nl1, Nc1, Dep, Fc, Type_Transform)
complex_float *Pict,*Result;
int Which_Filter;
int Nl, Nc, Nl1, Nc1;
int Dep, Type_Transform;
float Fc;
{
    int i, j, ind;
    float  u,v,Coef;
    int Nl_2, Nc_2;
    float pyr_2d_cf_filter ();

    Nl_2 = Nl1 / 2.;
    Nc_2 = Nc1 / 2.;
    
    for (i = 0; i < Nl1; i++)
    {
        u = (float) (Dep * (i - Nl_2)); 
        for (j = 0; j < Nc1; j++)
        {
            ind = i * Nc1 + j;
            v = (float) (Dep * (j - Nc_2));
            if ((u < - Nl_2 ) || (u >= Nl_2) ||
                 (v < - Nc_2) || (v >= Nc_2))
            {
                switch (Which_Filter)
                {
                    case FILTER_H:
                    case FILTER_H_TILDE:
                             Coef = 0.;
                             break;
                    case FILTER_G:
                    case FILTER_G_TILDE:
                             Coef = 1.;
                             break;
                    default:
                             printf ("Pb: Unknown Filter\n");
                             Coef = 0.;
                             break;
                }
            }
            else Coef = pyr_2d_cf_filter (Which_Filter, u, v, Fc, 
                                         Nl, Nc, Type_Transform);
            Result [ind].re = Pict [ind].re * Coef;
            Result [ind].im = Pict [ind].im * Coef;
        }
    }
} 

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

pyr_2d_cf_resol_down (Im_Hf, Im_Bf, Tab_Nl, Tab_Col,  Num_Etap, 
                      Fc, Type_Transform)
complex_float *Im_Hf,*Im_Bf;
int *Tab_Nl, *Tab_Col;
int Num_Etap, Type_Transform;
float Fc;
{
    int Nl,Nc,Nl1,Nc1;
    int Dep;
  
    Nl = Tab_Nl [0];
    Nc = Tab_Col [0];
    Nl1 = Tab_Nl [Num_Etap];
    Nc1 = Tab_Col [Num_Etap];
    
    /* the frequencies Im_Hf = Im_Bf * FILTER_G */
    INT_POW (2,Num_Etap,Dep);
    pyr_2d_cf_mult_tab_imag (Im_Bf, Im_Hf, FILTER_G, Nl, Nc, 
                             Nl1, Nc1, Dep, Fc, Type_Transform);

    /* the low frequencies Im_Bf = Im_Bf * FILTER_H */
    pyr_2d_cf_mult_tab_imag (Im_Bf, Im_Bf, FILTER_H, Nl, Nc, 
                             Nl1, Nc1, Dep, Fc, Type_Transform);

    /* unedersampling of the low frequencies */
    reduction (Im_Bf, Nl1, Nc1, Tab_Nl [Num_Etap+1], Tab_Col [Num_Etap+1]);
}

/***************************************************************************/    
pyr_2d_cf_tfo (Imag, Pyramid, Tab_Nl, Tab_Col, Tab_Pos, 
               Nbr_Etap, Fc, Type_Transform)
complex_float *Imag,*Pyramid;
int *Tab_Nl, *Tab_Col, *Tab_Pos;
int Nbr_Etap, Type_Transform;
float Fc;
{
    int i,j,Nl1,Nc1;
    int Pos;
    complex_float *Im_Hf;
        
    for (i = 0; i < Nbr_Etap; i++)
    {
        Im_Hf = Pyramid + Tab_Pos [i];
        /* Compute the wavelet coefficients Im_Hf and the image at a lower 
           resolution 
        */
        pyr_2d_cf_resol_down (Im_Hf, Imag, Tab_Nl, Tab_Col, 
                              i, Fc, Type_Transform);
    }
    
    Nl1 = Tab_Nl [Nbr_Etap];
    Nc1 = Tab_Col [Nbr_Etap];
    Pos = Tab_Pos [Nbr_Etap];

    /* copy in the pyramid the image at a lower resolution */
    for (j = 0; j < Nl1*Nc1; j++)
    {
        Pyramid [Pos+j].re = Imag [j].re;
        Pyramid [Pos+j].im = Imag [j].im;
    }
}


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

pyr_2d_cf_shanon_interpolate (Imag, Imag_up, Nl, Nc, Nl_up, Nc_up)
complex_float *Imag, *Imag_up;
int Nl, Nc, Nl_up, Nc_up;
{
    int i,j,u,v,ind;
   
    /* initialisation */
    for (i = 0; i < Nl_up * Nc_up; i++) Imag_up [i].re = Imag_up [i].im = 0.;
    
    /* Copy imag in the center of Imag */
    for (i = 0; i < Nl; i++)
    {
        u = i - Nl / 2;
        for (j = 0; j < Nc; j++)
        {
            v = j - Nc / 2;
            ind = (u  + Nl_up / 2) * Nc_up + v + Nc_up / 2;
            Imag_up [ind].re = Imag [i * Nc + j].re;
            Imag_up [ind].im = Imag [i * Nc + j].im;
        }
    }
}

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

pyr_2d_cf_extra_wave (Pyramid, Tab_Nl, Tab_Col, Tab_Pos, Nbr_Etap, 
                      Pyr_Up, Fc, Type_Transform)
complex_float *Pyramid, *Pyr_Up;
int *Tab_Nl, *Tab_Col, *Tab_Pos;
int Nbr_Etap, Type_Transform;
float Fc;
{
    int i,Nl,Nc,Nl1,Nc1;
    complex_float *Imag, *Imag_Up, *Im_Hf, *Im_Bf;
    
    /* Interpolation on the wavelet planes
       The first plane of Pyr_Up = second plane of Pyramid ... */
    for (i = 1; i < Nbr_Etap; i++)
    {
        Nl1 = Tab_Nl [i];
        Nc1 = Tab_Col [i];
        Nl = Tab_Nl [i-1];
        Nc = Tab_Col [i-1];
        Imag_Up = Pyr_Up + Tab_Pos [i-1];
        Imag = Pyramid + Tab_Pos [i];
        pyr_2d_cf_shanon_interpolate (Imag, Imag_Up, Nl1, Nc1, Nl, Nc);
    }

    /* Compute the wavelet plane associated to the low resolution plane */

    Nl1 = Tab_Nl [Nbr_Etap];
    Nc1 = Tab_Col [Nbr_Etap];
    Nl = Tab_Nl [Nbr_Etap-1];
    Nc = Tab_Col [Nbr_Etap-1];

    Im_Hf = cf_vector_alloc (Nl * Nc);
    Im_Bf = cf_vector_alloc (Nl * Nc);
    
    Imag = Pyramid + Tab_Pos [Nbr_Etap];
    Imag_Up = Pyr_Up + Tab_Pos [Nbr_Etap];
    for (i = 1; i < Nl1*Nc1; i++)
    {
        Im_Bf [i].re = Imag [i].re;
        Im_Bf [i].im = Imag [i].im;
        Imag_Up [i].re = Imag [i].re;
        Imag_Up [i].im = Imag [i].im;
    }
    Imag_Up = Pyr_Up + Tab_Pos [Nbr_Etap-1];
    
    pyr_2d_cf_resol_down (Im_Hf, Im_Bf, Tab_Nl, Tab_Col, 
                          Nbr_Etap, Fc, Type_Transform);
    pyr_2d_cf_shanon_interpolate (Im_Hf, Imag_Up, Nl1, Nc1, Nl, Nc);
    
    free ((char *) Im_Hf); 
    free ((char *) Im_Bf);  
}

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

pyr_2d_cf_resol_up (Im_Hf, Im_Bf, Im_Up, Tab_Nl, Tab_Col, 
                    Num_Etap, Fc, Type_Transform)
/* Build in the Fourier space the image at a resolution
   Num_Etap from the TF of the wavelet plane Im_Hf and the
   TF of the plane at a lower resolution */
complex_float *Im_Hf,*Im_Bf,*Im_Up;
int *Tab_Nl, *Tab_Col;
int Num_Etap, Type_Transform;
float Fc;
{
    int Nl,Nc,Nl1,Nc1,i,j,u,v,ind,ind1;
    int Nl2,Nc2;
    int Dep;

    Nl = Tab_Nl [0];
    Nc = Tab_Col [0];

    Nl1 = Tab_Nl [Num_Etap];
    Nc1 = Tab_Col [Num_Etap];
    
    Nl2 = Tab_Nl [Num_Etap+1];
    Nc2 = Tab_Col [Num_Etap+1];
    
    /* Calcule Im_Hf * FILTER_G_TILDE + Im_Bf * FILTER_H_TILDE */
    INT_POW (2,Num_Etap,Dep);

    /* Im_Hf * FILTER_G_TILDE */
    pyr_2d_cf_mult_tab_imag (Im_Hf, Im_Up, FILTER_G_TILDE, 
                                    Nl, Nc, Nl1, Nc1, Dep, Fc, Type_Transform);

    /* Im_Bf * FILTER_H_TILDE */
    pyr_2d_cf_mult_tab_imag (Im_Bf, Im_Bf, FILTER_H_TILDE, 
                                    Nl, Nc, Nl2, Nc2, Dep, Fc, Type_Transform);
    
    /* addition */
    for (i = 0; i < Nl2; i++)
    {
        u = i - Nl2 / 2;
        for (j = 0; j < Nl2; j++)
        {
            v = j - Nl2 / 2;
            ind = i * Nc2 + j;
            ind1 = (u + Nl1 /2) * Nc1 + v + Nc1 / 2;
            Im_Up [ind1].re +=  Im_Bf [ind].re;
            Im_Up [ind1].im +=  Im_Bf [ind].im;
        }
    }
}

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

pyr_2d_cf_build  (Imag, Pyramid, Tab_Nl, Tab_Col, Tab_Pos, 
                  Nbr_Etap, Fc, Type_Transform)
complex_float *Pyramid;
complex_float *Imag;
int *Tab_Nl, *Tab_Col, *Tab_Pos;
int Nbr_Etap, Type_Transform;
float Fc;
{
    int i,j,Nl1,Nc1;
    int Pos;
    complex_float *Im_Hf, *Imag_Up;
        
    i = Nbr_Etap;
    Nl1 = Tab_Nl [i];
    Nc1 = Tab_Col [i];
    Pos = Tab_Pos [i];
    Imag_Up = cf_vector_alloc (Tab_Nl [0] * Tab_Col [0]);

    /* copy the low resolution plane in the image */
    for (j = 0; j < Nl1*Nc1; j++)
    {
        Imag [j].re = Pyramid [Pos+j].re;
        Imag [j].im = Pyramid [Pos+j].im;
    }

    /* Iterative Reconstruction */
    for (i = Nbr_Etap-1; i >= 0 ; i--)
    {
        Nl1 = Tab_Nl [i];
        Nc1 = Tab_Col [i];
        Pos = Tab_Pos [i];
        Im_Hf = Pyramid + Pos;

        pyr_2d_cf_resol_up (Im_Hf, Imag, Imag_Up, Tab_Nl, 
                            Tab_Col, i, Fc, Type_Transform);

        for (j = 0; j < Nl1*Nc1; j++)
        {
            Imag [j].re = Imag_Up [j].re;
            Imag [j].im = Imag_Up [j].im;
        }
    }
    free ((char *) Imag_Up);
}

/***************************************************************************/
           
pyr_2d_cf_build_direct (Imag, Pyramid, Tab_Nl, Tab_Col, Tab_Pos, Nbr_Etap)
complex_float *Pyramid;
complex_float *Imag;
int *Tab_Nl, *Tab_Col, *Tab_Pos;
int Nbr_Etap;
{
    int i,j,Nl1,Nc1,Nl2,Nc2;
    int Pos;
    complex_float *Ptr_Pyr, *Imag_Up;
        
    i = Nbr_Etap;
    Nl1 = Tab_Nl [i];
    Nc1 = Tab_Col [i];
    Pos = Tab_Pos [i];
    Imag_Up = cf_vector_alloc (Tab_Nl [0] * Tab_Col [0]);

   /* copy the low resolution plane in the image */
    for (j = 0; j < Nl1*Nc1; j++)
    {
        Imag [j].re = Pyramid [Pos+j].re;
        Imag [j].im = Pyramid [Pos+j].im;
    }

    for (i = Nbr_Etap-1; i >= 0 ; i--)
    {
        Nl2 = Tab_Nl [i+1];
        Nc2 = Tab_Col [i+1];

        Nl1 = Tab_Nl [i];
        Nc1 = Tab_Col [i];
        Pos = Tab_Pos [i];
        Ptr_Pyr = Pyramid + Pos;

        /* We give the same size to two images */
        pyr_2d_cf_shanon_interpolate (Imag, Imag_Up, Nl2, Nc2, Nl1, Nc1);

        /* addition */
        for (j = 0; j < Nl1*Nc1; j++)
        {
            Imag [j].re = Imag_Up [j].re + Ptr_Pyr[j].re;
            Imag [j].im = Imag_Up [j].im + Ptr_Pyr[j].im;
        }
    }
    free ((char *) Imag_Up);
}

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

pyr_2d_cf_transform (Imag, Pyr, Tab_Nl, Tab_Col, Tab_Pos,
                     Nbr_Plan,Type_Transform,Fc)
float *Imag, *Pyr;
int *Tab_Nl, *Tab_Col, *Tab_Pos;
int Nbr_Plan, Type_Transform;
float Fc;
{
    int Size_Ima, Size_Pyr;
    int Nl,Nc;
    int Nbr_Etap = Nbr_Plan - 1;
    complex_float *Imag_cf, *Pyr_cf;

    Nl = Tab_Nl [0];
    Nc = Tab_Col[0];

    /* image FFT */
    Size_Ima = Nl*Nc;
    Imag_cf = cf_vector_alloc (Size_Ima);
    prepare_fft_real (Imag, Imag_cf, Nl);
    ft_cf_any_power_of_2 (Imag_cf, 1,  Nl);

    Size_Pyr = Tab_Pos [Nbr_Etap] + Tab_Nl[Nbr_Etap] * Tab_Col[Nbr_Etap];
    Pyr_cf = cf_vector_alloc (Size_Pyr);

    /* wavelet transform of the image */
    pyr_2d_cf_tfo (Imag_cf, Pyr_cf, Tab_Nl, Tab_Col, Tab_Pos, 
                   Nbr_Etap, Fc, Type_Transform);

    /* inverse Fourier transform of the pyramid */
    pyr_2d_cf_fft_2d (Pyr_cf, Tab_Nl, Tab_Col, Tab_Pos, Nbr_Plan, -1);

    /* Normalisation of the  pyramide */
    pyr_2d_cf_normalize_fft (Pyr_cf, Tab_Nl, Tab_Col, Tab_Pos, Nbr_Plan, -1);

    /* Copy the  result in the  pyramid */
    pyr_2d_cf_pyr_re (Pyr_cf, Tab_Nl, Tab_Col, Tab_Pos, Nbr_Plan, Pyr);

    free ((char *) Imag_cf);
    free ((char *) Pyr_cf);
}

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

pyr_2d_cf_build_pict_from_pyr (Pyr, Imag, Tab_Nl, Tab_Col, Tab_Pos, 
                               Nbr_Plan, Type_Transform, Fc, Build_Direct_Ok)
float *Imag, *Pyr;
int *Tab_Nl, *Tab_Col, *Tab_Pos;
int Nbr_Plan, Type_Transform, Build_Direct_Ok;
float Fc;
{
    int Size_Ima, Size_Pyr;
    int Nl,Nc,i;
    int Nbr_Etap = Nbr_Plan - 1;
    complex_float *Imag_cf, *Pyr_cf;

    Nl = Tab_Nl [0];
    Nc = Tab_Col [0];
    Size_Ima = Nl*Nc;

    Size_Pyr = Tab_Pos [Nbr_Etap] + Tab_Nl[Nbr_Etap] * Tab_Col[Nbr_Etap];
    Pyr_cf = cf_vector_alloc (Size_Pyr);
    
    /* Copy the pyramid to a complex  structure */
    for (i = 0; i < Size_Pyr; i++)
    {
        Pyr_cf [i].re = Pyr [i];
        Pyr_cf [i].im = 0.;
    }

    /* Normalisation */
    pyr_2d_cf_normalize_fft (Pyr_cf, Tab_Nl, Tab_Col, Tab_Pos, Nbr_Plan, 1);

    /* Pyramid FFT   */
    pyr_2d_cf_fft_2d (Pyr_cf, Tab_Nl, Tab_Col, Tab_Pos, Nbr_Plan, 1);

    /* Reconstruction  */
    Imag_cf = cf_vector_alloc (Size_Ima);
    if (Build_Direct_Ok == TRUE) 
            pyr_2d_cf_build_direct (Imag_cf, Pyr_cf, Tab_Nl, 
                                    Tab_Col, Tab_Pos, Nbr_Etap);
    else    pyr_2d_cf_build (Imag_cf, Pyr_cf, Tab_Nl, Tab_Col, Tab_Pos, 
                             Nbr_Etap, Fc, Type_Transform);

    /* inverse FFT of the image  */
    ft_cf_any_power_of_2 (Imag_cf, -1,  Nl);

    /* Copy of the result in Imag */
    for (i = 0; i < Size_Ima; i++) Imag[i] = Imag_cf[i].re;

    free ((char *) Imag_cf);
    free ((char *) Pyr_cf);
}

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

pyr_2d_cf_interp_plan (Pyr, Imag, Tab_Nl, Tab_Col, 
                       Tab_Pos, Num_Plan, Num_Plan_Interp)
float *Imag, *Pyr;
int *Tab_Nl, *Tab_Col, *Tab_Pos;
int Num_Plan, Num_Plan_Interp;
{
    int Nl,Nc,Nl_i,Nc_i,i;
    complex_float *Plan, *Plan_up;
    float *Ptr_Pyr;
 
    Nl = Tab_Nl[Num_Plan];
    Nc = Tab_Col[Num_Plan];
    Nl_i = Tab_Nl[Num_Plan_Interp];
    Nc_i = Tab_Col[Num_Plan_Interp];
    
    Plan = cf_vector_alloc (Nl*Nc);
    Plan_up = cf_vector_alloc (Nl_i*Nc_i);
    
    Ptr_Pyr = Pyr + Tab_Pos[Num_Plan];
    for (i = 0; i < Nl*Nc; i++)
    {
        Plan [i].re = Ptr_Pyr [i];
        Plan [i].im = 0.;
    }

    /* Fourier transform */
    pyr_2d_cf_normalise_plan (Plan, Nl, Tab_Nl[0], 1);
    ft_cf_any_power_of_2 (Plan, 1,  Nl);
    
    /* Interpolation */
    pyr_2d_cf_shanon_interpolate (Plan, Plan_up, Nl, Nc, Nl_i, Nc_i);
    
    /* inverse Fourier transform */
    ft_cf_any_power_of_2 (Plan_up, -1, Nl_i);
    pyr_2d_cf_normalise_plan (Plan_up, Nl_i, Tab_Nl[0], -1);

    /* Copy of the interpolated map in Imag */
    for (i = 0; i < Nl_i*Nc_i; i++)
    {
        Imag[i] = Plan_up [i].re;
    }
    
    free ((char *) Plan);
    free ((char *) Plan_up);
}