/* @(#)pyr.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.c
**
*******************************************************************************
**
**    DESCRIPTION  Pyramidal algorithm routines
**    -----------    
*******************************************************************************
**
** pyr_2d_difference_pict (Pict, Pyr, Nl, Nc, Ind_Pyramid)
** float *Pict, *Pyr
** int Nl, Nc, Ind_Pyramid
**
** Difference between  an image and a pyramid plane. The result is backup in the
** pyramid plane
**         Pyramid(Plan) = Pyramid(Plan) - Image
** Pict = Image
** Pyr = Pyramid
** Nl = Line number of the image and the pyramid plane 
** Nc = Column number of the image and the pyramid plane 
** Ind_Pyramid = position of the first pixel plane in the pyramid
**
*******************************************************************************
**
** pyr_2d_move_plan_pyr (Pyramid, Tab_Nl, Tab_Col, Tab_Pos, Num_Plan)
** float *Pyramid
** int *Tab_Nl, *Tab_Col, *Tab_Pos, Num_Plan
**
** Copies a pyramid plane to the upper pyramid plane with data compression
** Memory must be allowed before this call 
** Pyramid = Pyramid
** Tab_Nl = line number array of each plane 
** Tab_Col = column number array of each plane 
** Tab_Pos = position array of each plane 
** Num_Plan = Plane number to copy (Num_Plan = 0 .. Nbr_Plan - 2)
**
*******************************************************************************
**
** pyr_2d_copy_pyr_to_pict (Pyr, Pict, Nl, Nc, Pos)
** float *Pyr, *Pict
** int Nl, Nc, Pos
**
** Copies a pyramid plane in an image 
** Pyr  = Pyramid
** Pict = image (output data)
** Nl = Line number of the image and the pyramid plane 
** Nc = Column number of the image and the pyramid plane 
** Pos =  position of the first pixel plane in the pyramid
**
*******************************************************************************
**
** pyr_2d_copy_pict_to_pyr (Pyr, Pict, Nl, Nc, Pos)
** float *Pyr, *Pict
** int Nl, Nc, Pos
**
** Copies an image in a pyramid plane
** Pyr  = Pyramid
** Pict = image
** Nl = Line number of the image and the pyramid plane 
** Nc = Column number of the image and the pyramid plane 
** Pos = position of the first pixel plane in the pyramid
**
*******************************************************************************
**
** pyr_2d_smoothing_picture (Imag, Smooth, Nl, Nc, Ind_Pyramid)
** float *Imag, *Smooth
** int Nl, Nc, Ind_Pyramid
**
** Smoothes an image. The smoothed image is placed in the pyramid
** Pict = image 
** Smooth = Pyramid where is placed the smoothed image 
** Nl = Line number of the image and the pyramid plane 
** Nc = Column number of the image and the pyramid plane 
** Ind_Pyramid = position of the first pixel plane in the pyramid
**
*******************************************************************************
**
** pyr_2d_position_ind (Tab_Nl, Tab_Col, Tab_Pos, Nl, Nc, Nbr_Etap)
** int *Tab_Nl, *Tab_Col, *Tab_Pos, Nbr_Etap, Nl, Nc
**
** Computes the arrays which contain the line, column and index number
** of each plane from an image of size NlxNc and for pyramid with a
** number of plane equal to  Nbr_Etap + 1
**
** Tab_Nl = line number array of each plane 
** Tab_Col = column number array of each plane 
** Tab_Pos = position array of each plane
** Nl = Image line number
** Nc = Imag column number
** Nbr_Etap =  Number of Planes - 1
**
*******************************************************************************
**
** pyr_2d_pyramid_build (Pict, Tab_Nl, Tab_Col, Tab_Pos, Pyramid, 
**                       Nbr_Etap, Type_Transform)
** float *Pict, *Pyramid
** int *Tab_Nl, *Tab_Col, *Tab_Pos, Nbr_Etap, Type_Transform
** 
** Builds a pyramid from an image, The pyramid memory must be allowed
**
** Pict = image
** Tab_Nl = line number array of each plane 
** Tab_Col = column number array of each plane 
** Tab_Pos = position array of each plane
** Pyramid = pyramid created.
** Nbr_Etap =  Number of planes - 1
** Type_Transform = TO_PYR_LINEAR ==>  Linear scaling function
** Type_Transform = TO_PYR_BSPLINE ==> B3-spline scaling function
**
*******************************************************************************
**
** pyr_2d_add_pict_pyr (Image, Pyramid, Nl, Nc, Pos)
** float *Image, *Pyramid
** int Nl, Nc, Pos
**
** Adds a pyramid plane to an image
**              Image = Image + Pyramid plane
** Image = Image 
** Pyramid = Pyramid
** Nl = Line number of the image and the pyramid plane 
** Nc = Column number of the image and the pyramid plane 
** Pos = position of the first pixel plane in the pyramid
**
*******************************************************************************
**
** pyr_2d_interpol_pyr (Pict, Pyramid, Nl, Nc, Nl2, Nc2, Pos)
** float *Pict, *Pyramid
** int Nl, Nc, Nl2, Nc2, Pos
**
** Interpolates a pyramid plane and plane it in Pict
**
** Pict = image
** Pyramid  = Pyramid
** Nl = Line number after interpolation
** Nc = Column number after interpoplation
** Nl2 = Line number of the pyramid plane
** Nc2 = Column number of the pyramid plane
** Pos = position of the first pixel plane in the pyramid
**
*******************************************************************************
**
** pyr_2d_build_pict_from_pyr(Pyramid, Pict, Tab_Nl, Tab_Col, 
**                            Tab_Pos, Nbr_Etap, Type_Transform)
** float *Pict, *Pyramid
** int *Tab_Nl, *Tab_Col, *Tab_Pos, Nbr_Etap, Type_Transform
**
** Reconstructs an image from its pyramid. Non iterative reconstruction
**
** Pyramid = Pyramid
** Pict = reconstructed image
** Tab_Nl = line number array of each plane 
** Tab_Col = column number array of each plane 
** Tab_Pos = position array of each plane
** Nbr_Etap =  Number of planes  - 1
** Type_Transform = TO_PYR_LINEAR ==>  Linear scaling function
** Type_Transform = TO_PYR_BSPLINE ==> B3-spline scaling function
**
*******************************************************************************
**
** int pyr_2d_size_pyramide (Nl, Nc, Nbr_Plan)
** int Nl, Nc, Nbr_Plan
**
** Computes the size of a pyramid 
**
** Nl = line number of the image from which the pyramid results
** Nc = column number of the image from which the pyramid results
** Nbr_Plan = number of scales = number of planes
**
*******************************************************************************
**
** pyr_2d_dist_pyr (Pyr1, Pyr2, Size, Distance)
** float *Pyr1, *Pyr2
** int Size
** float *Distance
**
** Computes the distance betwwen two pyramids
**
** Pyr1 = first pyramid
** Pyr2 = second pyramid
** Size = size of the two pyramids
** Distance = calculated distance
**
*******************************************************************************
**
** pyr_2d_build_iter_pict_from_pyr (Pyramid, Pict, Tab_Nl, Tab_Col, 
**                                  Tab_Pos, Nbr_Etap, Nbr_Iter,Type_Transform)
** float *Pict, *Pyramid
** int *Tab_Nl, *Tab_Col, *Tab_Pos, Nbr_Etap,Nbr_Iter
**
** Reconstructs an image from its pyramid. Iterative reconstruction
**
** Pyramid = Pyramid
** Pict = reconstructed image
** Tab_Nl = line number array of each plane 
** Tab_Col = column number array of each plane 
** Tab_Pos = position array of each plane
** Nbr_Etap = Number of planes - 1
** Nbr_Iter = number of iterations
** Type_Transform = TO_PYR_LINEAR ==>  Linear scaling function
** Type_Transform = TO_PYR_BSPLINE ==> B3-spline scaling function
**
*******************************************************************************
** 
** pyr_2d_create_file_pict_from_pyramid (File_Name_In, File_Name_Out, Nbr_Iter)
** char *File_Name_In, *File_Name_Out
** int Nbr_Iter
**
** Creates a file wich contains an image from the file wich contains
** its pyramid.
**
** File_Name_In = pyramid file name
** File_Name_Out = image file name
** Nbr_Iter = number of iterations in the iterative reconstruction
**
*******************************************************************************
** 
** pyr_2d_info (Pyr, Nl, Nc, Pos, Min, Max, Moy, Sig)
** float *Pyr, *Min, *Max, *Moy, *Sig
** int Nl, Nc, Pos
**
** Computes the min,max,average and standard deviation of a pyramid plane
**
** Pyr = pyramid
** Nl = Line number of the pyramid plane 
** Nc = Column number of the pyramid plane 
** Pos = position of the first pixel plane in the pyramid
** Min = minimum value of the plane 
** Max = maximum value of the plane
** Moy = average
** Sig = standard  deviation
**
*******************************************************************************
** 
** pyr_2d_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;
**
** Interpolates a pyramid plane to a given resolution and places it in an
** image
**
** Pyr = pyramid
** Tab_Nl = line number array of each plane 
** Tab_Col = column number array of each plane 
** Tab_Pos = position array of each plane
** Num_Plan = plane number to interpolate
** Num_Plan_Interp = plane number of the resolution to reach
** Imag = image
**
******************************************************************************/ 

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

/* B3-spline coefficient values */

#define VAL_1_SUR_16 1./16.
#define VAL_1_SUR_4 1./4.
#define VAL_3_SUR_8 3./8.

#define VISU_PRINT FALSE

/* coefficients array size */
#define SIZE_TAB_COEF_INTER 10 
float Tab_Coef_Inter[SIZE_TAB_COEF_INTER] = {.00245,
                                             -0.00915,
                                             0.03414,
                                             -0.12720,
                                             0.60048,
                                             0.60048,
                                             -0.12740,
                                             0.03414,
                                             -0.00915,
                                             0.00245};

float Coef_Bspline[5] = {VAL_1_SUR_16, 
                         VAL_1_SUR_4, 
                         VAL_3_SUR_8, 
                         VAL_1_SUR_4, 
                         VAL_1_SUR_16};


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

pyr_2d_difference_pict (Pict1, Pict2, Nl, Nc, Ind_Pyramid)
float *Pict1, *Pict2;
int Nl, Nc, Ind_Pyramid;
{
    int i,j;
    int indi = 0;
    int Ind_Pyr = Ind_Pyramid;
    
    for (i = 0; i < Nl; i++)
    {
        for (j = 0; j < Nc; j++)
        {
            Pict2 [Ind_Pyr] = Pict1 [indi] - Pict2 [Ind_Pyr];
            indi++;
            Ind_Pyr++;
        }
    }
}

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

pyr_2d_move_plan_pyr (Pyramid, Tab_Nl, Tab_Col, Tab_Pos, Etap)
float *Pyramid;
int *Tab_Nl, *Tab_Col, *Tab_Pos, Etap;
{
    int i,j,Pos_Etap1,Pos_Etap2;
    
    Pos_Etap2 = Tab_Pos [Etap+1];
    for (i = 0; i < Tab_Nl [Etap]; i += 2)
    {
       Pos_Etap1 = Tab_Pos [Etap] + Tab_Col [Etap] * i;
       for (j = 0; j < Tab_Col [Etap]; j += 2)
       {
           Pyramid [Pos_Etap2] = Pyramid [Pos_Etap1];
           Pos_Etap1 += 2;
           Pos_Etap2 ++;
        }
    }
}

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

pyr_2d_copy_pyr_to_pict (Pyr, Pict, Nl, Nc, Pos)
float *Pyr, *Pict;
int Nl, Nc, Pos;
{
    int i,j; 
    int ind1 = Pos;
    int ind2 = 0;
       
    for (i = 0; i < Nl; i++)
    {
        for (j = 0; j < Nc; j++)
        {
            Pict [ind2] = Pyr [ind1];
            ind1 ++;
            ind2 ++;
        }
    }
}
/***************************************************************************/

pyr_2d_copy_pict_to_pyr (Pyr, Pict, Nl, Nc, Pos)
float *Pyr, *Pict;
int Nl, Nc, Pos;
{
    int i,j; 
    int ind1 = Pos;
    int ind2 = 0;
       
    for (i = 0; i < Nl; i++)
    {
        for (j = 0; j < Nc; j++)
        {
            Pyr [ind1] = Pict [ind2];
            ind1 ++;
            ind2 ++;
        }
    }
}

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

pyr_2d_smoothing_picture (Imag, Smooth, Nl, Nc, Ind_Pyramid)
float *Imag, *Smooth;
int Nl, Nc, Ind_Pyramid;
{
    int i,j;
    int im,ii,ip,jm,jp;
    int indi = Ind_Pyramid;
    
    for (i = 0; i < Nl; i++)
    {
        if (i == 0) im = 0;
        else im = (i-1) * Nc;
        if (i >= Nl-1) ip = i * Nc;
        else ip = (i + 1) * Nc;
        ii = i * Nc;
        
        for (j = 0; j < Nc; j++)
        {
            if (j == 0) jm = j;
            else jm = j - 1;
            if (j >= Nc-1) jp = j;
            else jp = j + 1;
            
            Smooth [indi] = 1./16. * (Imag[im + jm]
                                     + Imag[im + jp] 
                                     + Imag[ip + jp]
                                     + Imag[ip + jm])
                          + 1./8. * (Imag[im + j]
                                     + Imag[ii + jp] 
                                     + Imag[ii + jm]
                                     + Imag[ip + j])
                          + 1./4. * Imag[ii+j];
            indi++;
        }
    }    
}            
        
/***************************************************************************/

pyr_2d_position_ind (Tab_Nl, Tab_Col, Tab_Pos, Nl, Nc, Nbr_Etap)
int *Tab_Nl, *Tab_Col, *Tab_Pos, Nbr_Etap, Nl, Nc;
{
    int i;
    
    Tab_Nl [0] = Nl;
    Tab_Col [0] = Nc;
    Tab_Pos [0] = 0;

    for (i = 1; i <= Nbr_Etap; i++)
    {
        Tab_Nl [i] = (Tab_Nl [i-1] - 1) / 2 + 1;
        Tab_Col [i] = (Tab_Col [i-1] - 1) / 2 + 1;
        Tab_Pos [i] = Tab_Pos [i-1] + Tab_Nl [i-1] * Tab_Col [i-1];
    }
}
 
/***************************************************************************/

pyr_2d_pyramid_build (Ima, Tab_Nl, Tab_Col, Tab_Pos, Pyramid, 
                      Nbr_Etap, Type_Transform)
float *Ima, *Pyramid;
int *Tab_Nl, *Tab_Col, *Tab_Pos, Nbr_Etap, Type_Transform;
{
    int i,Size;
    float *Pict;

    Size = Tab_Nl [0]*Tab_Col [0];   
    Pict = f_vector_alloc (Size);
    for (i=0;i<Size;i++) Pict[i] = Ima[i];

    for (i = 0; i < Nbr_Etap; i++)
    {
        /* smooth the image */
        if (Type_Transform == TO_PYR_LINEAR)
        {
            pyr_2d_smoothing_picture (Pict, 
                                      Pyramid, 
                                      Tab_Nl [i], 
                                      Tab_Col [i], 
                                      Tab_Pos [i]);
        }
        else
        {
            pyr_2d_bspline_smoothing_pict (Pict, 
                                              Pyramid, 
                                              Tab_Nl [i], 
                                              Tab_Col [i], 
                                              Tab_Pos [i]);
        }

        /* move the plane inside the pyramid */
        pyr_2d_move_plan_pyr (Pyramid, Tab_Nl, Tab_Col, Tab_Pos, i);
        
        /* difference between an image and itself smoothed */
        pyr_2d_difference_pict (Pict, 
                                Pyramid, 
                                Tab_Nl [i], 
                                Tab_Col [i], 
                                Tab_Pos [i]);

        /* copy the image in the framework */
        pyr_2d_copy_pyr_to_pict (Pyramid, 
                                 Pict,Tab_Nl [i+1],
                                 Tab_Col[i+1],
                                 Tab_Pos[i+1]);

    }
    free ((char *) Pict);
}

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

pyr_2d_add_pict_pyr (Image, Pyramid, Nl, Nc, Pos)
float *Image, *Pyramid;
int Nl, Nc, Pos;
{
    int i,j; 
    int ind_pyr = Pos;
    int ind = 0;
       
    for (i = 0; i < Nl; i++)
    {
        for (j = 0; j < Nc; j++)
        {
            Image [ind] += Pyramid [ind_pyr];
            ind ++;
            ind_pyr ++;
        }
    }
}

/***************************************************************************/            
static calcul_ind10 (x,Nx,Tab_Dep)
int x, Nx, *Tab_Dep;
{
    int i;
    
    for (i = 0; i < 4; i ++)
    {
        Tab_Dep [i] = x + i - 4;
        if (Tab_Dep [i] < 0) Tab_Dep [i] = 0;
    }
    Tab_Dep [4] = x;
    for (i = 5; i < SIZE_TAB_COEF_INTER; i ++)
    {
        Tab_Dep [i] = x + i - 4;
        if (Tab_Dep [i] >= Nx) Tab_Dep [i] = Nx - 1;
    }
}
              
/***************************************************************************/

pyr_2d_interpol_pyr (Pict, Pyramid, Nl, Nc, Nl2, Nc2, Pos)
float *Pict, *Pyramid;
int Nl, Nc, Nl2, Nc2, Pos;
{
    int i,j,k;
    int Pos_I,Tab_Dep[SIZE_TAB_COEF_INTER];
    int ind_dep_pyr,Ind_Pict,Ind_Pyr;
    
    /* Pixels at node interpolation */
    for (i = 0; i < Nl2; i ++)
    {
        for (j = 0; j < Nc2; j ++)
        {
            Ind_Pict = 2 * i * Nc + 2 * j;
            Ind_Pyr = Pos + i * Nc2 + j;
            Pict [Ind_Pict] = Pyramid [Ind_Pyr];
        }
    
        /* pixel node in column */
        Pos_I = i * Nc2 + Pos;
        for (j = 0; j < Nc2; j ++)
        {
            calcul_ind10 (j,Nc2,Tab_Dep);
            Ind_Pict = 2 * i * Nc + 2 * j + 1;
            Pict [Ind_Pict] = 0.;
            for (k = 0; k < SIZE_TAB_COEF_INTER; k ++)
            {
                ind_dep_pyr = Pos_I + Tab_Dep[ k];
                Pict [Ind_Pict] += Tab_Coef_Inter [k] * Pyramid [ind_dep_pyr];
            }
        }
    }
    
    /* pixel node in line */
    for (i = 0; i < Nl2; i ++)
    {
        calcul_ind10 (i,Nl2,Tab_Dep);
        for (j = 0; j < Nc2; j ++)
        {
            Ind_Pict = (2 * i + 1) * Nc + 2 * j;
            Pict [Ind_Pict] = 0.;
            for (k = 0; k < SIZE_TAB_COEF_INTER; k ++)
            {
                ind_dep_pyr = Pos + Tab_Dep[ k] * Nc2 + j;
                Pict [Ind_Pict] += Tab_Coef_Inter [k] * Pyramid [ind_dep_pyr];
            }
        }
        
        /* pixel node in line and column */
        for (j = 0; j < Nc2; j ++)
        {
            Ind_Pict = (2 * i + 1) * Nc + 2 * j + 1;
            Pict [Ind_Pict] = 0.;
            calcul_ind10 (j,Nc2,Tab_Dep);
            for (k = 0; k < SIZE_TAB_COEF_INTER; k ++)
            {
                ind_dep_pyr = Ind_Pict + 2 * (Tab_Dep[ k] - j) - 1;
                Pict [Ind_Pict] += Tab_Coef_Inter [k] * Pict [ind_dep_pyr];
            }
        }
    }
}


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

pyr_2d_build_pict_from_pyr (Pyr, Pict, Tab_Nl, Tab_Col, Tab_Pos, Nbr_Etap)
float *Pict, *Pyr;
int *Tab_Nl, *Tab_Col, *Tab_Pos, Nbr_Etap;
{
    int i,Size;
    float *Pyramid;
   
    Size = wave_io_size_pyr (Tab_Nl[0], Tab_Col[0], Nbr_Etap+1);
    Pyramid = f_vector_alloc(Size);
    for (i = 0; i < Size; i++) Pyramid[i] = Pyr[i];

    for (i = Nbr_Etap - 1; i >= 0; i--)
    {
        /* pyramid interpolation plane */

        pyr_2d_interpol_pyr (Pict, 
                             Pyramid, 
                             Tab_Nl [i], 
                             Tab_Col [i], 
                             Tab_Nl [i+1], 
                             Tab_Col [i+1],
                             Tab_Pos [i+1]);
        
        /* addition of the pyramid plane to the image */
        pyr_2d_add_pict_pyr (Pict, Pyramid, Tab_Nl[i], Tab_Col[i], Tab_Pos[i]);
        
        /* copy the image to the pyramid */
        pyr_2d_copy_pict_to_pyr (Pyramid, 
                                 Pict, Tab_Nl[i], 
                                 Tab_Col [i], 
                                 Tab_Pos [i]);
    }
    free ((char *) Pyramid);
}

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

pyr_2d_dist_pyr (Pyr1, Pyr2, Size, Distance)
float *Pyr1, *Pyr2;
int Size;
float *Distance;
{
    int j;
    float Dist = 0.;

    for (j = 0; j < Size; j++) 
    {
        Dist += fabs((Pyr2 [j] - Pyr1 [j]));
    }
    Dist /= (float) Size;
    *Distance = Dist;
}

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

pyr_2d_build_iter_pict_from_pyr (Pyramid, 
                                 Pict, 
                                 Tab_Nl, 
                                 Tab_Col, 
                                 Tab_Pos, 
                                 Nbr_Etap,
                                 Nbr_Iter,
                                 Type_Transform)
float *Pict, *Pyramid;
int *Tab_Nl, *Tab_Col, *Tab_Pos, Nbr_Etap,Nbr_Iter,Type_Transform;
{
    int i,j;
    int Size;
    float *Pyr1,*Pyr2,*Pyr3;
    float Distance;
    
    Size = Tab_Pos [Nbr_Etap] + Tab_Nl[Nbr_Etap] * Tab_Col[Nbr_Etap];
    Pyr1 = f_vector_alloc(Size);
    Pyr2 = f_vector_alloc(Size);
    Pyr3 = f_vector_alloc(Size);
    
    for (j = 0; j < Size; j++) Pyr1 [j] = Pyramid [j];
    
    for (i = 0; i < Nbr_Iter; i++)
    {
        for (j = 0; j < Size; j++) Pyr3 [j] = Pyr1 [j];
        pyr_2d_build_pict_from_pyr (Pyr3, 
                                    Pict, 
                                    Tab_Nl, 
                                    Tab_Col, 
                                    Tab_Pos, 
                                    Nbr_Etap);
        pyr_2d_pyramid_build (Pict, 
                              Tab_Nl, 
                              Tab_Col, 
                              Tab_Pos, 
                              Pyr2, 
                              Nbr_Etap,
                              Type_Transform);
        
        pyr_2d_dist_pyr (Pyramid, Pyr2, Size, &Distance);
#if VISU_PRINT
        fprintf (stderr, "Iter %d : Distance = %f\n", i, Distance);
#endif
        for (j = 0; j < Size; j++)
        {
            Pyr1 [j] += Pyramid [j] - Pyr2 [j];
        }
    }
    pyr_2d_build_pict_from_pyr (Pyr1, 
                                Pict, 
                                Tab_Nl, 
                                Tab_Col, 
                                Tab_Pos, 
                                Nbr_Etap);
    free ((char *) Pyr1);
    free ((char *) Pyr2);
    free ((char *) Pyr3);
}

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

pyr_2d_info (Pyr, Nl, Nc, Pos, Min, Max, Moy, Sig)
float *Pyr, *Min, *Max, *Moy, *Sig;
int Nl, Nc, Pos;
{
    int i,j,ind;
    float Val;
    float Minimum,Maximum,Moyenne,Sigma;
    
    ind = Pos;
    Minimum = Pyr [ind];
    Maximum = Pyr [ind];
    Moyenne = 0.;
    Sigma = 0.;
    for (i = 0; i < Nl; i++)
    {
        for (j = 0; j < Nc; j++)
        {
            Val = Pyr [ind];
            if (Val < Minimum) Minimum = Val;
            if (Val > Maximum) Maximum = Val;
            Moyenne += Val;
            ind ++;
        }
    }
    Moyenne = Moyenne / ((float) (Nl * Nc));
    for (i = 0; i < Nl*Nc; i++)
      {
       Sigma += (Pyr [Pos + i] - Moyenne) * (Pyr [Pos + i] - Moyenne);
      } 
    Sigma /= (float)(Nl*Nc);
    Sigma = sqrt (Sigma);
    
    *Min = Minimum;
    *Max = Maximum;
    *Moy = Moyenne;
    *Sig = Sigma;
}

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

static test_coef ()
{
    int i,j;
    float Energ = 0.;
    
    for (i = 0; i < 5; i++)
    {
        for (j = 0; j < 5; j++)
        {
            printf ("%f ",Coef_Bspline[i]*Coef_Bspline[j]);
            Energ += Coef_Bspline[i]*Coef_Bspline[j];
        }
        printf ("\n");
    }
    printf ("Energ = %f\n",Energ);
}

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

pyr_2d_bspline_smoothing_pict(Imag, Smooth, Nl, Nc, Ind_Pyramid)
float *Imag, *Smooth;
int Nl, Nc, Ind_Pyramid;
{
    int i,j,k,l;
    int ii,jj;
    int indi = Ind_Pyramid;
    float Val;

    for (i = 0; i < Nl; i++)
    {
        for (j = 0; j < Nc; j++)
        {
            Val = 0.;
            for (k = 0; k < 5; k++)
            {
                ii = i - k + 2;
                if (ii < 0) ii = 0;
                else
                {
                    if (ii >= Nl) ii = Nl - 1;
                }
                for (l = 0; l < 5; l++)
                {
                    jj = j - l + 2;
                    if (jj < 0) jj = 0;
                    else
                    {
                       if (jj >= Nc) jj = Nc - 1;
                    }
                    Val += Coef_Bspline[l] * Coef_Bspline[k] * Imag[ii*Nc + jj];
                }
            }
            Smooth [indi] = Val;
            indi++;
        }
    }
}

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

pyr_2d_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,j;
    float *Plan, *Plan_up;
    int Size = Tab_Nl[Num_Plan_Interp]*Tab_Col[Num_Plan_Interp];

    /* memory allocation */
    Plan = f_vector_alloc (Size);
    Plan_up = f_vector_alloc (Size);
    
    /* copy the pyramid plane to a buffer */
    for (i = 0; i < Tab_Nl[Num_Plan]*Tab_Col[Num_Plan]; i++)
    {
        Plan_up[i] = Plan[i] = Pyr[Tab_Pos[Num_Plan]+i];
    }

    /* scale by scale interpolation */
    for (i = Num_Plan; i > Num_Plan_Interp; i--)
    {
        Nl = Tab_Nl[i];
        Nc = Tab_Col[i];
        Nl_i = Tab_Nl[i-1];
        Nc_i = Tab_Col[i-1];
        for (j = 0; j < Nl*Nc; j++) Plan[j] = Plan_up[j];
        pyr_2d_interpol_pyr (Plan_up, Plan, Nl_i, Nc_i, Nl, Nc, 0);
    }

    /* Copy the result to Imag */
    for (j = 0; j < Size; j++)   Imag[j] = Plan_up[j];

    /* Deallocate the memory */
    free ((char *) Plan);
    free ((char *) Plan_up);
}

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