/* @(#)extr_plan.c	17.1.1.1 (ES0-DMD) 01/25/02 17:21:21 */
/*===========================================================================
  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:  extr_plan.c
**
*******************************************************************************
**
**    DESCRIPTION  Extraction of a plane from a wavelet transform
**    ----------- 
**                 
******************************************************************************
**
** wavelet_extract_plan_file (File_Name_Imag, File_Name_Transform, Num_Plan)
** char *File_Name_Imag, *File_Name_Transform;
** int Num_Plan;
** 
** Extracts a plane from a wavelet transform which is in the file of name
** File_Name_Transform. The plane is saved in a image of name File_Name_Imag
** 
** File_Name_Imag = output file name of the image
** File_Name_Transform = input file name of the wavelet transform
** Num_Plan = plane number to extract (Num_Plan = 1 .. Number_of_Planes)
**
******************************************************************************
**
** wavelet_extract_plan (Wavelet, Imag, Nl, Nc, Num_Plan)
** float *Imag;
** wave_transf_des *Wavelet;
** int *Nl, *Nc, Num_Plan;
**
** Extracts a plane from a wavelet transform 
** 
** Wavelet = INPUT: wavelet transform
** Imag = OUTPUT: extracted plane
** Nl = OUTPUT: ptr to the number of lines
** Nc = OUTPUT: ptr to the number of columns
** Num_Plan = plane number to extract (Num_Plan = 1 .. Number_of_Planes)
** 
** The Mallat's wavelet transform is not treated by this routine 
** 
******************************************************************************
**
** wavelet_extract_plan_mallat (Wavelet, Imag, Imag_H, Imag_D, Imag_V,
**                             Nl_Plan, Nc_Plan, Num_Plan)
** float **Imag, **Imag_H, **Imag_D, **Imag_V;
** wave_transf_des *Wavelet;
** int *Nl_Plan, *Nc_Plan, Num_Plan;
**
** Extracts a  plane from a Mallat's wavelet transform 
**    We need:  Wavelet -> Type_Wave_Transform == TO_MALLAT_BARLAUD
**
** Wavelet = INPUT: wavelet transform
** Imag = OUTPUT: extracted plane
** Nl_Plan = OUTPUT: ptr to the number of lines of Imag
** Nc_Plan = OUTPUT: ptr to the number of columns of Imag
** Num_Plan = plane number to extract (Num_Plan = 1 .. Number_of_Planes)
** Imag_H = image which contains horizontals details of the plane number
** Imag_D = image which contains diagonals details of the plane number
** Imag_V = image which contains verticals details of the plane number
**
** The size of these tree images is Nl_Plan/2 x Nc_Plan/2
**
******************************************************************************
**
** wavelet_pointer_plan (Wavelet, Ptr, Nl, Nc, Num_Plan, Choose)
** float **Ptr;
** wave_transf_des *Wavelet;
** int *Nl, *Nc, Num_Plan;
** int Choose;
**
** return a pointer on the plane number Num_Plan of the wavelet transform
** 
** Wavelet = INPUT: wavelet transform
** Ptr = OUTPUT: pointer on a plane
** Nl = OUTPUT: ptr to the number of lines
** Nc = OUTPUT: ptr to the number of columns
** Num_Plan = plane number to point (Num_Plan = 1 .. Number_of_Planes)
** Choose = only use for the mallat transform
**	             Choose =  1 ==>  *Ptr = Ptr_Mallat -> Low_Resol
**	             Choose =  2 ==>  *Ptr = Ptr_Mallat -> Coef_Horiz
**	             Choose =  3 ==>  *Ptr = Ptr_Mallat -> Coef_Vert
**	             Choose =  4 ==>  *Ptr = Ptr_Mallat -> Coef_Diag
**
******************************************************************************/

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

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

wavelet_extract_plan_file (File_Name_Imag, File_Name_Transform, Num_Plan)
char *File_Name_Imag, *File_Name_Transform;
int Num_Plan;
{
    float *Imag, *Imag_H, *Imag_D, *Imag_V;
    wave_transf_des Wavelet;
    int Nl, Nc;

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

    /* plane extraction */
    if (Wavelet.Type_Wave_Transform == TO_MALLAT_BARLAUD)
    {
        wavelet_extract_plan_mallat (&Wavelet, &Imag, &Imag_H, &Imag_D, &Imag_V,
                             &Nl, &Nc, Num_Plan);

        /* write the details images
        if (Num_Plan < Nbr_Plan)
        {
            io_write_pict_f_to_file ("m_horiz", Imag_H, Nl/2, Nc/2);
            io_write_pict_f_to_file ("m_diag", Imag_D, Nl/2, Nc/2);
            io_write_pict_f_to_file ("m_vert", Imag_V, Nl/2, Nc/2);
        } */
    }
    else wavelet_extract_plan (&Wavelet, &Imag, &Nl, &Nc, Num_Plan);

    /* save the result */
    io_write_pict_f_to_file (File_Name_Imag, Imag, Nl, Nc);

    wave_io_free (&Wavelet);
    free ((char *) Imag);
}

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

wavelet_extract_plan (Wavelet, Imag, Nl, Nc, Num_Plan)
float **Imag;
wave_transf_des *Wavelet;
int *Nl, *Nc, Num_Plan;
{
    float *Ptr;
    int i, Nbr_Plan;
    string Mes_Memory;

    /* number of planes */
    Nbr_Plan = Wavelet -> Nbr_Plan;

    /* Test the plane number Num_Plan */
    if ((Num_Plan <= 0) || (Num_Plan > Nbr_Plan))
    {
        /* The plane number Num_Plan is false */
        sprintf (Mes_Memory,", Number of scales = %d\n", Nbr_Plan);
        io_err_message_exit (ERR_PLANE_NUMBER, Mes_Memory);
    }

    switch (Wavelet->Type_Wave_Transform)
    {
        case TO_PAVE_LINEAR:
        case TO_PAVE_BSPLINE:
        case TO_PAVE_BSPLINE_FFT:
                   *Nl = Wavelet->Nbr_Ligne;
                   *Nc = Wavelet->Nbr_Col;
                   Ptr = Wavelet->Pave.Data + (Num_Plan - 1) * *Nl * *Nc;
                   *Imag = f_vector_alloc (*Nl * *Nc);
                   for (i = 0; i < *Nl * *Nc; i++) (*Imag)[i] = Ptr[i];
                   break;
        case TO_PYR_LINEAR:
        case TO_PYR_BSPLINE:
        case TO_PYR_FFT_DIFF_RESOL:
        case TO_PYR_FFT_DIFF_SQUARE_RESOL:
                   *Nl = (Wavelet->Pyramid.Tab_Nl) [Num_Plan - 1];
                   *Nc = (Wavelet->Pyramid.Tab_Col)[Num_Plan - 1];
                    Ptr = Wavelet->Pyramid.Data +
                                    (Wavelet->Pyramid.Tab_Pos)[Num_Plan - 1];
                   *Imag = f_vector_alloc (*Nl * *Nc);
                   for (i = 0; i < *Nl * *Nc; i++) (*Imag)[i] = Ptr[i];
                   break;
        case TO_MALLAT_BARLAUD:
                  fprintf (stderr,"Error: See wavelet_extract_plan_mallat\n");
                   break;
        default:
                    io_err_message_exit (ERR_TRANSF, " ");
                    break;
    }   
}

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

wavelet_extract_plan_mallat (Wavelet, Imag, Imag_H, Imag_D, Imag_V,
                             Nl_Plan, Nc_Plan, Num_Plan)
float **Imag, **Imag_H, **Imag_D, **Imag_V;
wave_transf_des *Wavelet;
int *Nl_Plan, *Nc_Plan, Num_Plan;
{
    struct mallat_plan_des *Ptr_Mallat;
    int Num_Etap, i, Size, Nbr_Plan;
    string Mes_Memory;


    /* number of planes */
    Nbr_Plan = Wavelet -> Nbr_Plan;

    /* Test the plane number Num_Plan */
    if ((Num_Plan <= 0) || (Num_Plan > Nbr_Plan))
    {
        /* The plane number Num_Plan is false */
        sprintf (Mes_Memory,", Number of scales = %d\n", Nbr_Plan);
        io_err_message_exit (ERR_PLANE_NUMBER, Mes_Memory);
    }

    Ptr_Mallat = &(Wavelet -> Mallat);

    /* if Num_Plan = Nbr_Plan, there is only one plane to extract,
       the last resolution plane */

    if (Num_Plan == Nbr_Plan)
    {
        /* we go to the last plane */
        for (Num_Etap = 1; Num_Etap < Nbr_Plan - 1; Num_Etap++)
        {
            Ptr_Mallat = Ptr_Mallat -> Smooth_Imag;
        }
        *Nl_Plan = Ptr_Mallat -> Nl;
        *Nc_Plan = Ptr_Mallat -> Nc;

        /* dynamic allocation for the plane to extract */
        *Imag = f_vector_alloc (*Nl_Plan * *Nc_Plan);

        /* copy of the last resolution plane to Imag */
        for (i = 0; i < *Nl_Plan * *Nc_Plan; i++)
                  (*Imag)[i] =  (Ptr_Mallat -> Low_Resol) [i];
    }
    else
    {
        /* we go the plane to extract */
        for (Num_Etap = 1; Num_Etap < Num_Plan; Num_Etap++)
        {
            Ptr_Mallat = Ptr_Mallat -> Smooth_Imag;
        }

        /* size of the image which constains this plane plus the next */
        *Nl_Plan = Ptr_Mallat -> Nl * 2;
        *Nc_Plan = Ptr_Mallat -> Nc * 2;
        *Imag = f_vector_alloc (*Nl_Plan * *Nc_Plan);

        /* detail image size */
        Size =  Ptr_Mallat -> Nl *  Ptr_Mallat -> Nc;
        *Imag_D = f_vector_alloc (Size);
        *Imag_H = f_vector_alloc (Size);
        *Imag_V = f_vector_alloc (Size);

        /* extraction of the image + the details */
        mallat_2d_extract_plan (*Imag, *Nl_Plan, *Nc_Plan, 
                                *Imag_H, *Imag_D, *Imag_V, 
                                Ptr_Mallat, Nbr_Plan - Num_Plan + 1);
    }
}

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

wavelet_pointer_plan (Wavelet, Ptr, Nl, Nc, Num_Plan, Choose)
float **Ptr;
wave_transf_des *Wavelet;
int *Nl, *Nc, Num_Plan;
int Choose;
{
    int i;
    struct mallat_plan_des *Ptr_Mallat;
    int Num_Etap;

    switch (Wavelet->Type_Wave_Transform)
    {
        case TO_PAVE_LINEAR:
        case TO_PAVE_BSPLINE:
        case TO_PAVE_BSPLINE_FFT:
                   *Nl = Wavelet->Nbr_Ligne;
                   *Nc = Wavelet->Nbr_Col;
                   *Ptr = Wavelet->Pave.Data + (Num_Plan - 1) * *Nl * *Nc;
                   break;
        case TO_PYR_LINEAR:
        case TO_PYR_BSPLINE:
        case TO_PYR_FFT_DIFF_RESOL:
        case TO_PYR_FFT_DIFF_SQUARE_RESOL:
                   *Nl = (Wavelet->Pyramid.Tab_Nl) [Num_Plan - 1];
                   *Nc = (Wavelet->Pyramid.Tab_Col)[Num_Plan - 1];
                    *Ptr = Wavelet->Pyramid.Data +
                                    (Wavelet->Pyramid.Tab_Pos)[Num_Plan - 1];
                   break;
        case TO_MALLAT_BARLAUD:
                   Ptr_Mallat = &(Wavelet -> Mallat);

                   /* Positionnement sur le plan a extraire */
                   for (Num_Etap = 1; Num_Etap < Num_Plan; Num_Etap++)
                     Ptr_Mallat = Ptr_Mallat -> Smooth_Imag;

                   /* Taille de l'image contenant ce plan + les suivant */
                   *Nl = Ptr_Mallat -> Nl * 2;
                   *Nc = Ptr_Mallat -> Nc * 2;
      
                  switch (Choose)
	          {
	             case 1: *Ptr = Ptr_Mallat -> Low_Resol;break;
	             case 2: *Ptr = Ptr_Mallat -> Coef_Horiz;break;
	             case 3: *Ptr = Ptr_Mallat -> Coef_Vert;break;
	             case 4: *Ptr = Ptr_Mallat -> Coef_Diag;break;
	          }
                  break;
        default:
                   printf ("Error: See wavelet_extract_plan_mallat\n");
                   break;
    }   
}

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