/* @(#)pave_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:  pave_cf.c
**
*******************************************************************************
**
**    DESCRIPTION  routines for the wavelet transform algorithm using
**    -----------  the FFT without reduction of sampling.
**
*******************************************************************************
**
** pave_2d_cf_tfo (Imag, Pave, Nl, Nc, Nbr_Plan, Tab_h, Tab_g)
** complex_float *Imag, *Pave;
** int Nl, Nc, Nbr_Plan;
** float *Tab_h, *Tab_g;
**
** Computes the wavelet transform of an image by using the FFT and does
** not reduce the sampling
**
** Imag = INPUT: Fourier transform of the image
** Pave = OUTPUT: Fourier transform of all the planes of the wavelet transform
** Nbr_Plan = INPUT: Scale number
** Tab_h, Tab_g = INPUT: Filters used for the wavelet transform
**
*******************************************************************************
**
** pave_2d_cf_down (I0, I1, C1, H, G, Etap, Nl, Nc)
** complex_float *I0, *I1, *C1;
** float *H, *G;
** int Etap, Nl, Nc;
**
** passage from a resolution to the next one:
** We computes from I0 an image I1 at a lower resolution and the
** wavelet coefficients C1 related to this resolution by using
** the filter H and G
**
** I0 = INPUT: Fourier transform of the image
** I1 = OUTPUT: Fourier transform of the image at the lower resolution
** C1 = OUTPUT: Fourier transform of the wavelet coefficients
** Etap = INPUT: Resolution number
** Nl, Nc = INPUT: lines and columns number 
**
*******************************************************************************
**
** pave_2d_cf_fft (Pave, Nl, Nc, Dir, Nbr_Plan)
** complex_float *Pave;
** int Nl, Nc;
** int Dir, Nbr_Plan;
**
** Computes the Fourier transform of all the planes of the wavelet transform
** Dir = 1  ==> direct Fourier transform
** Dir = -1 ==>  inverse Fourier transform 
**
** Pave = INPUT, OUTPUT: wavelet transform
** Nl, Nc = INPUT: lines and columns number 
** Dir = INPUT: direction
** Nbr_Plan = INPUT: Scale number
**
*******************************************************************************
**
** pave_2d_cf_transform (Imag, Pave, Nl, Nc, Nbr_Plan, Type_Wavelet, Fc)
** float *Imag, *Pyr;
** int *Tab_Nl, *Tab_Col, *Tab_Pos;
** int Nbr_Plan, Type_Wavelet;
** float Fc;
**
** Computes the wavelet transform  of an image
** 
** Image = INPUT: image
** Pave = OUTPUT: wavelet transform of the image
** Nl, Nc = INPUT: lines and columns number 
** Nbr_Plan = INPUT: Scale number
** Type_Wavelet = INPUT: Wavelet algorithm type (ex TO_PAVE_BSPLINE_FFT)
** Fc = INPUT: frequency cut off
** 
*******************************************************************************
** 
** pave_2d_cf_build (Pave, Imag, Nl, Nc, Nbr_Plan)
** float *Imag, *Pave;
** int Nl, Nc, Nbr_Plan, Type_Wavelet;
** float Fc;
**
** Reconstructs an image from its wavelet transform
**
** Image = OUTPUT: reconstructed image
** Pave = INPUT: wavelet transform of the image
** Nl, Nc = INPUT: lines and columns number 
** Nbr_Plan = INPUT: Scale number
**
*******************************************************************************
** 
** pave_2d_cf_extract_plan (Pave, Imag, Nl, Nc, Num_Plan)
** float *Imag, *Pave;
** int Nl, Nc, Num_Plan;
** 
** Extracts a plane from the wavelet transform
**
** Pave = INPUT: wavelet transform 
** Imag = OUTPUT: extracted plane
** Nl, Nc = INPUT: lines and columns number 
** Num_Plan = INPUT: plane number
**
******************************************************************************/ 

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

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

pave_2d_cf_down (I0, I1, C1, H, G, Etap, Nl, Nc)
complex_float *I0, *I1, *C1;
float *H, *G;
int Nl,Nc;
{
    register int i,j,u,v,Nl_2,Nc_2;
    register int indi,indu,indv;
    int Dep;
    float Coef_h, Coef_g;

    Dep = pow(2., (float) Etap) + 0.5;
    Nl_2 = Nl / 2.;
    Nc_2 = Nc / 2.;
    indi = 0;
    for (i = 0; i < Nl; i++)
    {
         u = (Dep * (i - Nl_2));
         indu = (u + Nl_2) * Nc;
         for (j = 0; j < Nc; j++,indi++)
         {
             v = Dep * (j - Nc_2);
             indv = indu + v + Nc_2;
             if ((u < - Nl_2 ) || (u >= Nl_2) ||
                 (v < - Nc_2) || (v >= Nc_2))
             {
                 Coef_h = 0.;
                 Coef_g = 1.;
             }
             else
             {
                 Coef_h = H[indv];
                 Coef_g = G[indv];
             }
             I1 [indi].re = I0[indi].re * Coef_h;
             I1 [indi].im = I0[indi].im * Coef_h;
             C1 [indi].re = I0[indi].re * Coef_g;
             C1 [indi].im = I0[indi].im * Coef_g;
         }
    }
}

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

pave_2d_cf_tfo (Imag, Pave, Nl, Nc, Nbr_Plan, Tab_h, Tab_g)
complex_float *Imag, *Pave;
int Nl,Nc,Nbr_Plan;
float *Tab_h, *Tab_g;
{
    int i;
    complex_float *Ptr_I, *Ptr_Ip, *Ptr_Cp;

    /* Copy the image */
    for (i = 0; i < Nl * Nc; i++)
    {
        Pave[i].re = Imag[i].re;
        Pave[i].im = Imag[i].im;
    }

    /* wavelet transfrom scale by scale */
    for (i = 0; i < Nbr_Plan - 1; i++)
    {
        Ptr_I = Pave + i * Nl * Nc;
        Ptr_Ip = Pave + (i+1) * Nl * Nc;
        Ptr_Cp = Ptr_I;
        pave_2d_cf_down (Ptr_I, Ptr_Ip, Ptr_Cp, Tab_h, Tab_g, i, Nl, Nc);
    }
}

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

pave_2d_cf_fft (Pave, Nl, Nc, Dir, Nbr_Plan)
complex_float *Pave;
int Nl, Nc;
int Dir, Nbr_Plan;
{
    int i;
    complex_float *Ptr_Ci;

    for (i = 0; i < Nbr_Plan; i++)
    {
        Ptr_Ci = Pave + Nl * Nc * i;
        ft_cf_any_power_of_2 (Ptr_Ci, Dir,  Nl);
    }
}

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

pave_2d_cf_transform (Imag, Pave, Nl, Nc, Nbr_Plan, Type_Wavelet, Fc)
float *Imag, *Pave;
int Nl, Nc;
int Nbr_Plan, Type_Wavelet;
float Fc;
{
    int Size_Ima, Size_Pave, i;
    complex_float *Imag_cf, *Pave_cf;
    float *Tab_h, *Tab_g;

    Size_Ima = Nl*Nc;
    Tab_h = f_vector_alloc (Size_Ima);
    Tab_g = f_vector_alloc (Size_Ima);

    /* Creation of the filters H and G */
    pyr_2d_cf_create_filter (Fc, Nl, Nc, Tab_h, FILTER_H, Type_Wavelet);
    pyr_2d_cf_create_filter (Fc, Nl, Nc, Tab_g, FILTER_G, Type_Wavelet);

    /* Fourier transform of the image */
    Imag_cf = cf_vector_alloc (Size_Ima);
    prepare_fft_real (Imag, Imag_cf, Nl);
    ft_cf_any_power_of_2 (Imag_cf, 1,  Nl);

    /* memory allocation for the wavelet transform */
    Size_Pave = Nl * Nc * Nbr_Plan;
    Pave_cf = cf_vector_alloc (Size_Pave);

    /* wavelet transform of the image in the wavelet space */
    pave_2d_cf_tfo (Imag_cf, Pave_cf, Nl, Nc, Nbr_Plan, Tab_h, Tab_g);

    /* inverse Fourier transform of the wavelet planes */
    pave_2d_cf_fft (Pave_cf, Nl, Nc, -1, Nbr_Plan);

    /* Resuts in Pave */
    for (i = 0; i < Size_Pave; i++) Pave[i] = Pave_cf[i].re;

    free ((char *) Imag_cf);
    free ((char *) Pave_cf);
    free ((char *) Tab_h);
    free ((char *) Tab_g);
}

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

pave_2d_cf_build (Pave, Imag, Nl, Nc, Nbr_Plan)
float *Imag, *Pave;
int Nl, Nc, Nbr_Plan;
{
    int Num_Plan,i,Pos;
    float *Plan;

    for (i = 0; i < Nl*Nc; i++) Imag [i] = 0.;

    for (Num_Plan = Nbr_Plan - 1; Num_Plan >= 0; Num_Plan--)
    {
          Pos = Nl * Nc * Num_Plan;
          Plan = Pave + Pos;
          /* recopie de l'image dans le pave */
          for (i = 0; i < Nl*Nc; i++) Imag [i] += Plan [i];
    }
}

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

pave_2d_cf_extract_plan (Pave, Imag, Nl, Nc, Num_Plan)
float *Imag, *Pave;
int Nl, Nc, Num_Plan;
{
    int i,Pos;
    float *Plan;

    Pos = Nl * Nc * Num_Plan;
    Plan = Pave + Pos;

    /* Copy the plane to the image */
    for (i = 0; i < Nl*Nc; i++) Imag [i] = Plan [i];
}

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