/* @(#)transform.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: transform.c ** ******************************************************************************* ** ** DESCRIPTION This module contains wavelet transform routines ** ----------- ****************************************************************************** ** ** wavelet_transform_file (File_Name_Imag, File_Name_Transform, ** Type_Transform, Fc, Nbr_Plan) ** char *File_Name_Imag, *File_Name_Transform; ** int Type_Transform; ** float Fc; ** int Nbr_Plan; ** ** Computes the wavelet transform of the image with a file name File_Name_Imag ** and write the result in the file File_Name_Transform ** ** File_Name_Imag = File name of the imput image ** File_Name_Transform = File name of the output wavelet transform ** Type_Transform = wavelet transform algorithm number ** Fc = cut-off frequency if the algorithm use the FFT ** Nbr_Plan = number of scales ** ****************************************************************************** ** ** wavelet_transform_data (Imag, Nl, Nc, Wavelet, Type_Transform, Fc, Nbr_Plan) ** float *Imag; ** int Nl, Nc; ** wave_transf_des *Wavelet; ** int Type_Transform; ** float Fc; ** int Nbr_Plan; ** ** Computes the wavelet transform of the image Imag ** and write the result in the structure Wavelet ** ** Imag = INPUT:image ** Wavelet = OUTPUT:wavelet ** Nl,Nc = INPUT: number of lines and columns ** Type_Transform = INPUT:wavelet transform algorithm number ** Which_Algorithm = 1...7 ** 1 ==> a trous algorithm with a linear scaling function ** the wavelet function is the difference betwwen two resolutions ** 2 ==> a trous algorithm with a B3-spline scaling function ** the wavelet function is the difference betwwen two resolutions ** 3 ==> algorithm using the Fourier transform: ** without any reduction of the samples between two scales ** the Fourier transform of the scaling function is a b3-spline ** the wavelet function is the difference between two resolutions ** 4 ==> pyramidal algorithm in the direct space, with a linear ** scaling function ** 5 ==> pyramidal algorithm in the direct space, with a b3-spline ** scaling function ** 6 ==> algorithm using the Fourier transform: ** with a reduction of the samples between two scales ** the Fourier transform of the scaling function is a b3-spline ** the wavelet function is the difference between two resolutions ** 7 ==> algorithm using the Fourier transform: ** with a reduction of the samples between two scales ** the Fourier transform of the scaling function is a b3-spline ** the wavelet function is the difference between the square of ** two resolutions ** 8 ==> Mallat's Algorithm with biorthogonal filters. ** Fc = INPUT:cut-off frequency if the algorithm use the FFT ** Nbr_Plan = INPUT:number of scales ** ******************************************************************************/ static char sccsid[] = "@(#)transform.c 17.1.1.1 02/01/25 ESO @(#)"; #include #include #include #include "Def_Math.h" #include "Def_Mem.h" #include "Def_Wavelet.h" /*****************************************************************************/ wavelet_transform_file (File_Name_Imag, File_Name_Transform, Type_Transform, Fc, Nbr_Plan) char *File_Name_Imag, *File_Name_Transform; int Type_Transform; float Fc; int Nbr_Plan; { float *Imag; wave_transf_des Wavelet; int Nl, Nc; /* read the input image */ io_read_file_to_pict_f (File_Name_Imag, &Imag, &Nl, &Nc); strcpy (Wavelet.Name_Imag, File_Name_Imag); wavelet_transform_data (Imag, Nl, Nc, &Wavelet, Type_Transform, Fc, Nbr_Plan); wave_io_write (File_Name_Transform, &Wavelet); wave_io_free (&Wavelet); free ((char *) Imag); } /*****************************************************************************/ static test_line_column (Nl, Nc) int Nl, Nc; /* We must check Nl to make sure it is a power of 2 */ { int len_exp,temp; if (Nl != Nc) io_err_message_exit (ERR_IMAGE_SQUARE, " "); else { len_exp = (int)(0.3+log((double)(Nl))/(log(2.0))); INT_POW(2,len_exp,temp); if (Nl != temp) io_err_message_exit (ERR_POWER_OF_2, " "); } } /*****************************************************************************/ wavelet_transform_data (Imag, Nl, Nc, Wavelet, Type_Transform, Fc, Nbr_Plan) float *Imag; int Nl, Nc; wave_transf_des *Wavelet; int Type_Transform; float Fc; int Nbr_Plan; { float *Pyr, *Pave; double Exp; int Size,Min,temp; Wavelet->Nbr_Ligne = Nl; Wavelet->Nbr_Col = Nc; Wavelet->Nbr_Plan = Nbr_Plan; Wavelet->Type_Wave_Transform = Type_Transform; /* test if the number of planes is not too high */ Min = (Nl < Nc) ? Nl : Nc; Exp = (double) Nbr_Plan + 2.; temp = pow(2., Exp) + 0.5; if (Min < temp) io_err_message_exit (ERR_NUMBER_OF_PLANES, " "); switch (Type_Transform) { case TO_PAVE_LINEAR: case TO_PAVE_BSPLINE: Size = Nl * Nc * Nbr_Plan; Wavelet->Pave.Data = f_vector_alloc (Size); Pave = Wavelet->Pave.Data; pave_2d_tfo (Imag, Pave, Nl, Nc, Nbr_Plan, Type_Transform); break; case TO_PAVE_BSPLINE_FFT: /* We must check Nl and Nc to make sure it is a power of 2 */ test_line_column (Nl, Nc); Wavelet->Pave.Freq_Coup = Fc; Size = Nl * Nc * Nbr_Plan; Wavelet->Pave.Data = f_vector_alloc (Size); Pave = Wavelet->Pave.Data; pave_2d_cf_transform (Imag, Pave, Nl, Nc, Nbr_Plan, Type_Transform, Fc); break; case TO_PYR_LINEAR: case TO_PYR_BSPLINE: Size = wave_io_size_pyr (Nl, Nc, Nbr_Plan); Wavelet->Pyramid.Size = Size; Wavelet->Pyramid.Data = f_vector_alloc (Size); Pyr = Wavelet->Pyramid.Data; wave_io_position_ind_pyr (Wavelet->Pyramid.Tab_Nl, Wavelet->Pyramid.Tab_Col, Wavelet->Pyramid.Tab_Pos, Nl, Nc, Nbr_Plan-1); pyr_2d_pyramid_build (Imag, Wavelet->Pyramid.Tab_Nl, Wavelet->Pyramid.Tab_Col, Wavelet->Pyramid.Tab_Pos, Pyr, Nbr_Plan - 1, Type_Transform); break; case TO_PYR_FFT_DIFF_RESOL: /* We must check Nl and Nc to make sure it is a power of 2 */ test_line_column (Nl, Nc); Wavelet->Pyramid.Freq_Coup = Fc; Size = wave_io_size_pyr (Nl, Nc, Nbr_Plan); Wavelet->Pyramid.Size = Size; Wavelet->Pyramid.Data = f_vector_alloc (Size); Pyr = Wavelet->Pyramid.Data; wave_io_position_ind_pyr (Wavelet->Pyramid.Tab_Nl, Wavelet->Pyramid.Tab_Col, Wavelet->Pyramid.Tab_Pos, Nl, Nc, Nbr_Plan); pyr_2d_cf_transform (Imag, Pyr, Wavelet->Pyramid.Tab_Nl, Wavelet->Pyramid.Tab_Col, Wavelet->Pyramid.Tab_Pos, Nbr_Plan, Type_Transform, Fc); break; case TO_PYR_FFT_DIFF_SQUARE_RESOL: /* We must check Nl and Nc to make sure it is a power of 2 */ test_line_column (Nl, Nc); Wavelet->Pyramid.Freq_Coup = Fc; Size = wave_io_size_pyr (Nl, Nc, Nbr_Plan); Wavelet->Pyramid.Size = Size; Wavelet->Pyramid.Data = f_vector_alloc (Size); Pyr = Wavelet->Pyramid.Data; wave_io_position_ind_pyr (Wavelet->Pyramid.Tab_Nl, Wavelet->Pyramid.Tab_Col, Wavelet->Pyramid.Tab_Pos, Nl, Nc, Nbr_Plan); pyr_2d_cf_transform (Imag, Pyr, Wavelet->Pyramid.Tab_Nl, Wavelet->Pyramid.Tab_Col, Wavelet->Pyramid.Tab_Pos, Nbr_Plan, Type_Transform, Fc); break; case TO_MALLAT_BARLAUD: wave_io_mallat_alloc (&(Wavelet->Mallat), Nbr_Plan-1, Nl, Nc); mallat_2d_transform (Imag,&(Wavelet->Mallat), Nl, Nc, Nbr_Plan); break; default: io_err_message_exit (ERR_TRANSF, " "); break; } } /*****************************************************************************/