/* @(#)oper_ima.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: oper_ima.c ** ******************************************************************************* ** ** DESCRIPTION Tools for image processing ** ----------- ** ***************************************************************************** ** ** float lib_mat_correl (Tab_Im1, Tab_Im2, Nl, Nc) ** float *Tab_Im1, *Tab_Im2; ** int Nl, Nc; ** ** gives the correlation between two images ** ***************************************************************************** ** ** float lib_mat_moy (Imag, Nl, Nc) ** float *Imag; ** int Nl, Nc; ** ** gives the average of an image ** ******************************************************************************* ** ** float lib_mat_ecart_type (Imag, Nl, Nc) ** float *Imag; ** int Nl, Nc; ** ** gives the standard deviation ** ******************************************************************************* ** ** lib_mat_moy_ecart_type (Imag, Nl, Nc, Ecart, Moy) ** float *Imag; ** int Nl, Nc; ** float *Ecart, *Moy; ** ** gives the average and the standard deviation of an image ** ******************************************************************************* ** ** lib_mat_detect_snr (Nc,Nl,Image,Option,Nit,Moyenne,Sigma) ** int Nc,Nl; ** float *Image; ** int Option,Nit; ** float *Moyenne,*Sigma; ** ** Estimation of the noise in a image ** INPUT ** Image ** Nl,Nc ** Option ** if Option = 0, the average of the image should be equal to 0 ** Nit = iteration number (3 est une bonne valeur) ** OUTPUT ** *Moyenne = average of the image ** *Sigma = standard deviation of the noise ** ******************************************************************************* ** ** lib_mat_convolve (Imag1, Imag2, Imag3, Nl, Nc) ** float *Imag1, *Imag2, *Imag3; ** int Nl, Nc; ** ** Convolution of two images by using the Fourier transform ** Imag3 = Imag1 x Imag2 ** ******************************************************************************* ** ** lib_mat_convolve_direct (I1, Nl1, Nc1, I2, Nl2, Nc2, I3) ** float *I1, *I2, *I3; ** int Nl1, Nc1, Nl2, Nc2; ** ** Convolution of two images in the direct space ** Nl1, Nc1 = size of I1 ** Nl2, Nc2 = size of I2 ** with Nl2 < NL1 and Nc2 < Nc1 ** ** Results in I3 of size Nl1,Nc1 ** ******************************************************************************/ static char sccsid[] = "@(#)oper_ima.c 17.1.1.1 02/01/25 ESO @(#)"; #include #include #include "Def_Math.h" #include "Def_Mem.h" /***************************************************************************/ float lib_mat_correl (Tab_Im1, Tab_Im2, Nl, Nc) float *Tab_Im1, *Tab_Im2; int Nl, Nc; { float Sum_X2,Sum_Y2,Sum_XY; int i; float Coef; Sum_X2 = Sum_Y2 = Sum_XY = 0.; for (i = 0; i < Nl * Nc; i++) { Sum_X2 += Tab_Im1 [i] * Tab_Im1 [i]; Sum_Y2 += Tab_Im2 [i] * Tab_Im2 [i]; Sum_XY += Tab_Im1 [i] * Tab_Im2 [i]; } Coef = Sum_XY / sqrt (Sum_X2 * Sum_Y2); return (Coef); } /***************************************************************************/ float lib_mat_moy (Imag, Nl, Nc) float *Imag; int Nl, Nc; { int j; float Moy = 0.; for (j = 0; j < Nl * Nc; j++) { Moy += Imag [j]; } Moy /= (float) (Nl * Nc); return (Moy); } /***************************************************************************/ float lib_mat_ecart_type (Imag, Nl, Nc) float *Imag; int Nl, Nc; { int j; float Moy = 0., Ecart = 0.; for (j = 0; j < Nl * Nc; j++) { Moy += Imag [j]; Ecart += Imag [j] * Imag [j]; } Moy /= (float) (Nl * Nc); Ecart /= (float) (Nl * Nc); Ecart = sqrt(Ecart - Moy*Moy); return (Ecart); } /***************************************************************************/ lib_mat_moy_ecart_type (Imag, Nl, Nc, Ecart, Moy) float *Imag; int Nl, Nc; float *Ecart, *Moy; { int j; *Moy = *Ecart = 0.; for (j = 0; j < Nl * Nc; j++) { *Moy += Imag [j]; *Ecart += Imag [j] * Imag [j]; } *Moy /= (float) (Nl * Nc); *Ecart /= (float) (Nl * Nc); *Ecart = sqrt(*Ecart - *Moy * *Moy); } /***************************************************************************/ lib_mat_detect_snr (Nc,Nl,Image,Option,Nit,Moyenne,Sigma) int Nc,Nl; float *Image; int Option,Nit; float *Moyenne,*Sigma; { int It,Il; float S0,S1,S2,Sm=0,x; for (It = 0; It < Nit; It++) { S0 = S1 = S2 = 0.; for (Il = 0; Il < Nl*Nc; Il++) { x = Image[Il]; if ((It == 0) || (fabs(x - *Moyenne) < Sm)) { S0 ++; S1 += x; S2 += x*x; } } if (Option == 1) { *Moyenne = S1 / S0; *Sigma = sqrt(S2/S0- *Moyenne * *Moyenne); } else { *Moyenne = 0.; *Sigma = sqrt(S2/S0); } Sm = 3. * *Sigma; } } /***************************************************************************/ lib_mat_convolve (Imag1, Imag2, Imag3, Nl, Nc) float *Imag1, *Imag2, *Imag3; int Nl, Nc; { int i; complex_float *Fft_Pict1,*Fft_Pict2; complex_float Val; Fft_Pict1 = cf_vector_alloc (Nl * Nc); prepare_fft_real (Imag1, Fft_Pict1, Nl); ft_cf_any_power_of_2 (Fft_Pict1, 1, Nl); Fft_Pict2 = cf_vector_alloc (Nl * Nc); prepare_fft_real (Imag2, Fft_Pict2, Nl); ft_cf_any_power_of_2 (Fft_Pict2, 1, Nl); for (i = 0; i < Nl*Nc; i++) { CF_MLT(Fft_Pict1[i], Fft_Pict2[i], Val); Fft_Pict1[i].re = Val.re; Fft_Pict1[i].im = Val.im; } ft_cf_any_power_of_2 (Fft_Pict1,-1,Nl); for (i = 0; i < Nl*Nc; i++) Imag3[i] = Fft_Pict1[i].re; free ((char *) Fft_Pict1); free ((char *) Fft_Pict2); } /***************************************************************************/ lib_mat_convolve_direct (I1, Nl1, Nc1, I2, Nl2, Nc2, I3) float *I1, *I2, *I3; int Nl1, Nc1, Nl2, Nc2; { int i,j,k,l; int i1,j1; float Val; int Nl_2,Nc_2; float *Ptr_i1, *Ptr_i2; int ii,jj; Nl_2 = Nl2 / 2; Nc_2 = Nc2 / 2; for (i = 0; i < Nl1; i++) { ii = i + Nl_2; for (j = 0; j < Nc1; j++) { Val = 0.; jj = j + Nc_2; for (k = 0; k < Nl2; k++) { i1 = ii - k; if (i1 < 0) i1 = 0; else if (i1 >= Nl1) i1 = Nl1 - 1; Ptr_i1 = I1 + i1*Nc1; Ptr_i2 = I2 + k*Nc2; for (l = 0; l < Nc2; l++) { j1 = jj - l; if (j1 < 0) j1 = 0; else if (j1 >= Nc1) j1 = Nl1 - 1; Val += Ptr_i1[j1] * Ptr_i2[l]; } } I3 [i*Nc1+j] = Val; } } }