/* @(#)io_wave.c 17.1.1.1 (ES0-DMD) 01/25/02 17:21:22 */ /*=========================================================================== 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: io_wave.c ** ******************************************************************************* ** ** DESCRIPTION inpout output routines for the wavelet files ** ----------- ** ******************************************************************************* ** ** wave_io_read (File_Name, Wave_Trans) ** char *File_Name; ** wave_transf_des *Wave_Trans; ** ** read a wavelet transform file and stores the results in a structure ** ** File_Name = INPUT: file name ** Wave_Trans = OUTPUT: wavelet ** ******************************************************************************* ** ** wave_io_write (File_Name, Wave_Trans) ** char *File_Name; ** wave_transf_des *Wave_Trans; ** ** writes a wavelet transform in a file ** ** File_Name = INPUT: file name ** Wave_Trans = INPUT: wavelet ** ******************************************************************************* ** ** wave_io_free (Wave_Trans) ** wave_transf_des *Wave_Trans; ** ** deallocates the memory of a wavelet ** ******************************************************************************* ** ** wave_io_position_ind_pyr (Tab_Nl, Tab_Col, Tab_Pos, Nl, Nc, Nbr_Etap) ** int *Tab_Nl, *Tab_Col, *Tab_Pos; ** int Nl, Nc, Nbr_Etap; ** ** computes positions of index of each plane of a pyramidal transform ** Tab_Nl = OUTPUT: number of lines array ** Tab_Nc = OUTPUT: number of columns array ** Tab_Pos = OUTPUT: position array ** Nl, Nc = original image size ** Nbr_Etap = number of sclales -1 ** ** ex: the plane number 3 has Tab_Nl[2] lines, Tab_Col[2] columns ** and the first pixel is at Tab_Pos[2] in the buffer of the ** wavelet transform ** ******************************************************************************* ** ** int wave_io_size_pyr (Nl, Nc, Nbr_Plan) ** int Nl, Nc, Nbr_Plan; ** ** computes the pyramid sys ** ** Nl, Nc = original image size ** Nbr_Plan = number of sclales ** ** return the size ** ******************************************************************************* ** ** wave_io_mallat_alloc (Des_Wave, Nbr_Etap, Nl, Nc) ** struct mallat_plan_des *Des_Wave; ** int Nbr_Etap, Nl, Nc; ** ** allocates the memory of the structure Des_Wave for a Mallat's transform ** with a number of scales equal to Nbr_Etap + 1 ** ** Nl, Nc = original image size ** Nbr_Etap = number of sclales -1 ** ******************************************************************************* ** ** int wave_io_size_data (Nl, Nc, Nbr_Plan, Type_Wave_Transform) ** int Nl, Nc, Type_Wave_Transform ** ** computes the necessary memory size for a wavelet transform algorithm ** defined by Type_Wave_Transform and for a number of scales Nbr_Plan ** ** Nl, Nc = original image size ** Nbr_Plan = number of scales = number of planes ** Type_Wave_Transform = choosen wavelet transform algorithm ** ******************************************************************************* ** ** wave_io_alloc (Wave_Trans, Type_Transform, Nbr_Plan, Nl, Nc) ** wave_transf_des *Wave_Trans; ** int Type_Transform; ** int Nbr_Plan, Nl, Nc; ** ** allocates the memory for a wavelet transform algorithm ** defined by Type_Wave_Transform and for a number of scales Nbr_Plan ** ** Nl, Nc = original image size ** Nbr_Plan = number of scales = number of planes ** Type_Wave_Transform = choosen wavelet transform algorithm ** ******************************************************************************/ static char sccsid[] = "@(#)io_wave.c 17.1.1.1 02/01/25 ESO @(#)"; #include #include #include #include "Def_Math.h" #include "Def_Mem.h" #include "Def_Wavelet.h" /****************************************************************************/ wave_io_position_ind_pyr (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]; } } /****************************************************************************/ int wave_io_size_pyr (Nl, Nc, Nbr_Plan) int Nl, Nc, Nbr_Plan; { int Size,i; int Nl1,Nc1,Nl2,Nc2,Pos1,Pos2; /* extern char *sys_errlist []; extern int errno; */ Nl1 = Nl; Nc1 = Nc; Pos1 = 0; for (i = 1; i < Nbr_Plan; i++) { Nl2 = (Nl1 - 1) / 2 + 1; Nc2 = (Nc1 - 1) / 2 + 1; Pos2 = Pos1 + Nl1 * Nc1; Nl1 = Nl2; Nc1 = Nc2; Pos1 = Pos2; } Size = Pos2 + Nl2 * Nc2; return (Size); } /****************************************************************************/ int wave_io_size_data (Nl, Nc, Nbr_Plan, Type_Wave_Transform) int Nl, Nc, Type_Wave_Transform; { int Size; switch (Type_Wave_Transform) { case TO_PAVE_LINEAR: case TO_PAVE_BSPLINE: case TO_PAVE_BSPLINE_FFT: Size = Nbr_Plan * Nl * Nc; break; case TO_PYR_LINEAR: case TO_PYR_BSPLINE: case TO_PYR_FFT_DIFF_RESOL: case TO_PYR_FFT_DIFF_SQUARE_RESOL: Size = wave_io_size_pyr (Nl, Nc, Nbr_Plan); break; case TO_MALLAT_BARLAUD: Size = Nl * Nc; break; default : io_err_message_exit (ERR_TRANSF, " "); break; } return (Size); } /***************************************************************************/ static wave_io_name (File_Name_In, File_Name_Out) char *File_Name_In, *File_Name_Out; { int L; strcpy (File_Name_Out, File_Name_In); L = strlen (File_Name_In); if ((L < 5) || (File_Name_In[L-1] != 'e') || (File_Name_In[L-2] != 'v') || (File_Name_In[L-3] != 'a') || (File_Name_In[L-4] != 'w') || (File_Name_In[L-5] != '.')) { strcat (File_Name_Out, ".wave"); } } /****************************************************************************/ static wave_io_mallat_alloc_plan (Des_Wave, Num_Etap, Nl, Nc, Nbr_Etap) struct mallat_plan_des *Des_Wave; int Nbr_Etap, Nl, Nc, Num_Etap; { int Size; unsigned int Size_Struct; Size_Struct = sizeof(struct mallat_plan_des); Size = Nl * Nc; Des_Wave -> Coef_Horiz = f_vector_alloc (Size); Des_Wave -> Coef_Diag = f_vector_alloc (Size); Des_Wave -> Coef_Vert = f_vector_alloc (Size); if (Num_Etap == Nbr_Etap) { Des_Wave -> Low_Resol = f_vector_alloc(Size); Des_Wave -> Smooth_Imag = NULL; } else { Des_Wave -> Smooth_Imag = (struct mallat_plan_des *) calloc (Size_Struct,1); if (Des_Wave -> Smooth_Imag == NULL) io_err_message_exit (ERR_ALLOC_MEMO, " "); Des_Wave -> Low_Resol = NULL; } } /****************************************************************************/ static wave_io_read_mallat (Ptr_Mallat, Nbr_Etap, File_Des) struct mallat_plan_des *Ptr_Mallat; int Nbr_Etap; FILE *File_Des; { struct mallat_plan_des *Ptr_Des; int Size, Nbr; char *Ptr; int i,Nl,Nc; Ptr_Des = Ptr_Mallat; for (i = 1; i <= Nbr_Etap; i++) { Nl = Ptr_Des -> Nl; Nc = Ptr_Des -> Nc; /* memory allocation for one plane */ Size = Ptr_Des -> Nl * Ptr_Des -> Nc; wave_io_mallat_alloc_plan (Ptr_Des, i, Nl, Nc, Nbr_Etap); /* read the horizontal details */ Ptr = (char *) (Ptr_Des -> Coef_Horiz); Nbr = fread (Ptr, sizeof(float), Size, File_Des); if (Nbr <= 0) io_err_message_exit (ERR_READ_DATA, " "); /* read the diagonals details */ Ptr = (char *) (Ptr_Des -> Coef_Diag); Nbr = fread (Ptr, sizeof(float), Size, File_Des); if (Nbr <= 0) io_err_message_exit (ERR_READ_DATA, " "); /* read the verticals details */ Ptr = (char *) (Ptr_Des -> Coef_Vert); Nbr = fread (Ptr, sizeof(float), Size, File_Des); if (Nbr <= 0) io_err_message_exit (ERR_READ_DATA, " "); if (i < Nbr_Etap) { /* read the next structure */ Ptr = (char *) (Ptr_Des -> Smooth_Imag); Nbr = fread (Ptr, sizeof(struct mallat_plan_des), 1, File_Des); if (Nbr <= 0) io_err_message_exit (ERR_READ_DATA, " "); Ptr_Des = Ptr_Des -> Smooth_Imag; } else { /* read the low resolution plane */ Ptr = (char *) (Ptr_Des -> Low_Resol); Nbr = fread (Ptr, sizeof(float), Size, File_Des); if (Nbr <= 0) io_err_message_exit (ERR_READ_DATA, " "); } } } /***************************************************************************/ wave_io_read (File_Name_In, Wave_Trans) char *File_Name_In; wave_transf_des *Wave_Trans; { FILE *File_Des; int Nbr,Nl,Nc,Nbr_Plan,Size; pyramid_f_des Ptr_Pyr; char *Ptr; string File_Name; struct mallat_plan_des *Ptr_Mallat; /* find the exact file name which must finish by .wave */ wave_io_name (File_Name_In, File_Name); /* open the file */ File_Des = fopen (File_Name, "r"); if (File_Des == NULL) io_err_message_exit (ERR_OPEN_FILE, File_Name); /* read the descriptor */ Nbr = fread ((char *) Wave_Trans, sizeof(wave_transf_des), 1, File_Des); if (Nbr <= 0) io_err_message_exit (ERR_READ_DATA, " "); Nl = Wave_Trans -> Nbr_Ligne; Nc = Wave_Trans -> Nbr_Col; Nbr_Plan = Wave_Trans -> Nbr_Plan; switch (Wave_Trans -> Type_Wave_Transform) { case TO_PAVE_LINEAR: case TO_PAVE_BSPLINE: case TO_PAVE_BSPLINE_FFT: Size = Nbr_Plan * Nl * Nc; (Wave_Trans -> Pave).Data = f_vector_alloc(Size); Ptr = (char *) ((Wave_Trans -> Pave).Data); Nbr = fread (Ptr, sizeof(float), Size, File_Des); if (Nbr <= 0) io_err_message_exit (ERR_READ_DATA, " "); break; case TO_PYR_LINEAR: case TO_PYR_BSPLINE: case TO_PYR_FFT_DIFF_RESOL: case TO_PYR_FFT_DIFF_SQUARE_RESOL: Ptr_Pyr = Wave_Trans -> Pyramid; Size = Ptr_Pyr.Size; wave_io_position_ind_pyr (Ptr_Pyr.Tab_Nl, Ptr_Pyr.Tab_Col, Ptr_Pyr.Tab_Pos, Nl, Nc, Nbr_Plan-1); (Wave_Trans -> Pyramid).Data = f_vector_alloc (Size); Nbr = fread ((char *) (Wave_Trans -> Pyramid).Data, sizeof(float), Size, File_Des); if (Nbr <= 0) io_err_message_exit (ERR_READ_DATA, " "); break; case TO_MALLAT_BARLAUD: Ptr_Mallat = &(Wave_Trans -> Mallat); wave_io_read_mallat (Ptr_Mallat, Nbr_Plan - 1, File_Des); break; default : io_err_message_exit (ERR_TRANSF, " "); break; } if (fclose (File_Des) != 0) io_err_message_exit (ERR_CLOSE_FILE, File_Name); } /****************************************************************************/ static wave_io_write_mallat (Ptr_Mallat, Nbr_Etap, File_Des) struct mallat_plan_des *Ptr_Mallat; int Nbr_Etap; FILE *File_Des; { struct mallat_plan_des *Ptr_Des; int Size, Nbr; char *Ptr; Ptr_Des = Ptr_Mallat; /* write the horizontal details */ Size = Ptr_Des -> Nl * Ptr_Des -> Nc; Ptr = (char *) (Ptr_Des -> Coef_Horiz); Nbr = fwrite (Ptr, sizeof(float), Size, File_Des); if (Nbr <= 0) io_err_message_exit (ERR_WRITE_DATA, " "); /* write the diagonal details */ Ptr = (char *) (Ptr_Des -> Coef_Diag); Nbr = fwrite (Ptr, sizeof(float), Size, File_Des); if (Nbr <= 0) io_err_message_exit (ERR_WRITE_DATA, " "); /* write the vertical details */ Ptr = (char *) (Ptr_Des -> Coef_Vert); Nbr = fwrite (Ptr, sizeof(float), Size, File_Des); if (Nbr <= 0) io_err_message_exit (ERR_WRITE_DATA, " "); if (Nbr_Etap > 1) { /* write the next structure */ Ptr = (char *) (Ptr_Des -> Smooth_Imag); Nbr = fwrite (Ptr, sizeof(struct mallat_plan_des), 1, File_Des); if (Nbr <= 0) io_err_message_exit (ERR_WRITE_DATA, " "); /* recursive call for the next structure */ wave_io_write_mallat (Ptr_Des->Smooth_Imag, Nbr_Etap-1, File_Des); } else { /* write the low resolution plane */ Ptr = (char *) (Ptr_Des -> Low_Resol); Nbr = fwrite (Ptr, sizeof(float), Size, File_Des); if (Nbr <= 0) io_err_message_exit (ERR_WRITE_DATA, " "); } } /****************************************************************************/ wave_io_write (File_Name_In, Wave_Trans) char *File_Name_In; wave_transf_des *Wave_Trans; { FILE *File_Des; int Nbr,Nl,Nc,Nbr_Plan,Size; pyramid_f_des Ptr_Pyr; char *Ptr; string File_Name; struct mallat_plan_des *Ptr_Mallat; /* find the exact file name */ wave_io_name (File_Name_In, File_Name); /* open the file */ File_Des = fopen (File_Name, "w"); if (File_Des == NULL) io_err_message_exit (ERR_OPEN_FILE, File_Name); /* write the descriptor */ Nbr = fwrite ((char *) Wave_Trans, sizeof(wave_transf_des), 1, File_Des); if (Nbr <= 0) io_err_message_exit (ERR_WRITE_DATA, " "); Nl = Wave_Trans -> Nbr_Ligne; Nc = Wave_Trans -> Nbr_Col; Nbr_Plan = Wave_Trans -> Nbr_Plan; switch (Wave_Trans -> Type_Wave_Transform) { case TO_PAVE_LINEAR: case TO_PAVE_BSPLINE: case TO_PAVE_BSPLINE_FFT: Ptr = (char *) ((Wave_Trans -> Pave).Data); Size = Nbr_Plan * Nl * Nc; Nbr = fwrite (Ptr, sizeof(float), Size, File_Des); if (Nbr <= 0) io_err_message_exit (ERR_WRITE_DATA, " "); break; case TO_PYR_LINEAR: case TO_PYR_BSPLINE: case TO_PYR_FFT_DIFF_RESOL: case TO_PYR_FFT_DIFF_SQUARE_RESOL: Ptr_Pyr = Wave_Trans -> Pyramid; Size = Ptr_Pyr.Size; Ptr = (char *) Ptr_Pyr.Data; Nbr = fwrite (Ptr, sizeof(float), Size, File_Des); if (Nbr <= 0) io_err_message_exit (ERR_WRITE_DATA, " "); break; case TO_MALLAT_BARLAUD: Ptr_Mallat = &(Wave_Trans -> Mallat); wave_io_write_mallat (Ptr_Mallat, Nbr_Plan - 1, File_Des); break; default : io_err_message_exit (ERR_TRANSF, " "); break; } if (fclose (File_Des) != 0) io_err_message_exit (ERR_CLOSE_FILE, File_Name); } /****************************************************************************/ static wave_io_free_mallat (Des_Wave, Nbr_Etap) struct mallat_plan_des *Des_Wave; int Nbr_Etap; { struct mallat_plan_des *Ptr_Wave; Ptr_Wave = Des_Wave; if (Nbr_Etap > 1) { /* Liberation des buffers de details */ free ((char *) (Ptr_Wave -> Coef_Horiz)); free ((char *) (Ptr_Wave -> Coef_Diag)); free ((char *) (Ptr_Wave -> Coef_Vert)); /* recursive call for the next structure */ wave_io_free_mallat (Ptr_Wave -> Smooth_Imag, Nbr_Etap - 1); free ((char *) (Ptr_Wave -> Smooth_Imag)); } else { free ((char *) Ptr_Wave -> Coef_Horiz); free ((char *) Ptr_Wave -> Coef_Diag); free ((char *) Ptr_Wave -> Coef_Vert); free ((char *) Ptr_Wave -> Low_Resol); } } /****************************************************************************/ wave_io_free (Wave_Trans) wave_transf_des *Wave_Trans; { struct mallat_plan_des *Des_Mallat; int Nbr_Etap; switch (Wave_Trans -> Type_Wave_Transform) { case TO_PAVE_LINEAR: case TO_PAVE_BSPLINE: case TO_PAVE_BSPLINE_FFT: free ((char *) ((Wave_Trans -> Pave).Data)); break; case TO_PYR_LINEAR: case TO_PYR_BSPLINE: case TO_PYR_FFT_DIFF_RESOL: case TO_PYR_FFT_DIFF_SQUARE_RESOL: free ((char *) ((Wave_Trans -> Pyramid).Data)); break; case TO_MALLAT_BARLAUD: Nbr_Etap = Wave_Trans -> Nbr_Plan - 1; Des_Mallat = &(Wave_Trans -> Mallat); wave_io_free_mallat (Des_Mallat, Nbr_Etap); break; default : io_err_message_exit (ERR_TRANSF, " "); break; } } /****************************************************************************/ wave_io_mallat_alloc (Des_Wave, Nbr_Etap, Nl, Nc) struct mallat_plan_des *Des_Wave; int Nbr_Etap, Nl, Nc; { int Nl_Plan, Nc_Plan, i; struct mallat_plan_des *Ptr_Wave; Ptr_Wave = Des_Wave; Nl_Plan = Nl; Nc_Plan = Nc; for (i = 1; i <= Nbr_Etap; i++) { wave_io_mallat_alloc_plan (Ptr_Wave, i, Nl_Plan, Nc_Plan, Nbr_Etap); Ptr_Wave = Ptr_Wave -> Smooth_Imag; } } /****************************************************************************/ wave_io_alloc (Wave_Trans, Type_Transform, Nbr_Plan, Nl, Nc, Fc) wave_transf_des *Wave_Trans; int Type_Transform; int Nbr_Plan, Nl, Nc; float Fc; { int Size; Wave_Trans -> Nbr_Ligne = Nl; Wave_Trans -> Nbr_Col = Nc; Wave_Trans -> Nbr_Plan = Nbr_Plan; Wave_Trans -> Type_Wave_Transform = Type_Transform; Wave_Trans->Pyramid.Freq_Coup = Fc; switch (Wave_Trans -> Type_Wave_Transform) { case TO_PAVE_LINEAR: case TO_PAVE_BSPLINE: case TO_PAVE_BSPLINE_FFT: Size = Nbr_Plan * Nl * Nc; (Wave_Trans -> Pave).Data = f_vector_alloc(Size); break; case TO_PYR_LINEAR: case TO_PYR_BSPLINE: case TO_PYR_FFT_DIFF_RESOL: case TO_PYR_FFT_DIFF_SQUARE_RESOL: Size = wave_io_size_pyr (Nl, Nc, Nbr_Plan); Wave_Trans->Pyramid.Size = Size; wave_io_position_ind_pyr (Wave_Trans->Pyramid.Tab_Nl, Wave_Trans->Pyramid.Tab_Col, Wave_Trans->Pyramid.Tab_Pos, Nl, Nc, Nbr_Plan-1); (Wave_Trans -> Pyramid).Data = f_vector_alloc (Size); break; case TO_MALLAT_BARLAUD: wave_io_mallat_alloc (&(Wave_Trans->Mallat), Nbr_Plan-1, Nl, Nc); break; default : io_err_message_exit (ERR_TRANSF, " "); break; } } /****************************************************************************/