Program Stat_image Implicit none Character(Len=100) option(50) Character(Len=80) Nom,nom2,out,string,mouch Integer*4 stat,iBP,NIter,iPas,Nkeys,I,IIter,Iparam,Jparam,Kparam,Imask,Ios Real*4 FClip,Eps Real*4, dimension(:), allocatable :: StatVal Ibp=-32 Call getline(nkeys,option) if (nkeys.eq.0) then Write(*,*) 'Syntax : -i image -m mask -n iteration -c clip_factor -e Eps -s step -o outputfile -l logfile' Write(*,*) 'Some statistical parameters of an or a part of an image' Write(*,*) '-i image : name of the image to be process (no def.)' Write(*,*) '-m mask name of the optional mask (def. no mask)' Write(*,*) '-n Iteration: number of iterations with rejection (def. 10)' Write(*,*) '-c Clip_Factor: halfsize in sigma of the interval (def. 3)' Write(*,*) '-e Precision: precision thresholding (def. 0.001)' Write(*,*) '-s Step: step between two pixels (def. 1)' Write(*,*) '-o outputfile : file containing the last statistical values (def. stat.lst)' Write(*,*) '-l logfile : logfile (def. mouch.lis)' STOP endif Niter=10 Fclip=3. Eps=0.001 Ipas=1 Nom='' Nom2='' Out='stat.lst' Mouch='mouch.lis' do i=1,nkeys,2 if(option(i).eq.'-i')nom =option(i+1) if(option(i).eq.'-m')nom2=option(i+1) if(option(i).eq.'-n')Read(option(i+1),*) Niter if(option(i).eq.'-c')Read(option(i+1),*) Fclip if(option(i).eq.'-s')Read(option(i+1),*) iPas if(option(i).eq.'-e')Read(option(i+1),*) Eps if(option(i).eq.'-o')out=option(i+1) if(option(i).eq.'-l')mouch=option(i+1) enddo Call Debut('stat_image',mouch) Call WC('Some statistical parameters of an or a part of an image') if ( Nom2 == '' ) then imask=0 else call WC('Mask :'//Trim(nom2)) imask=1 endif call WC('Image File : '//Trim(nom)) call WI('Number of iterations',Niter,1) Call WR('Clip factor ',Fclip,1) Call WR('Precision limit ',Eps,1) Call WI('Step between pixels ',iPas,1) !************************************************************************** allocate(StatVal(15*NIter),stat=ios) if (ios /= 0) stop " STOP allocation error" call calst (Nom,Nom2,iBP,StatVal,NIter,FClip,Eps,iMask,iPas) open(unit=11,file=out) Call WC('Statistical values : '//Trim(out)) Call WC(' ') DO IIter = 1,NIter Write(string,"(' ---- Iter: ',I)")IIter Call WC(string) IParam = 1 + 15*(IIter-1) JParam = 6 + 15*(IIter-1) KParam = 11 + 15*(IIter-1) write(string,*),' NPix / fractions:',StatVal(IParam),& StatVal(JParam) !, & StatVal(KParam) call WC(string) IParam = 12 + 15*(IIter-1) JParam = 13 + 15*(IIter-1) Write(String,*)' Extrema :',StatVal(IParam),StatVal(JParam) Call WC(string) IParam = 2 + 15*(IIter-1) JParam = 7 + 15*(IIter-1) write(string,*),' Mean / err :',StatVal(IParam),StatVal(JParam) call WC(string) if(iiter.eq.niter) then write(11,*) StatVal(IParam) endif IParam = 3 + 15*(IIter-1) JParam = 8 + 15*(IIter-1) write(string,*),' Deviation / err :',StatVal(IParam),StatVal(JParam) call Wc(string) if(iiter.eq.niter) then write(11,*) StatVal(IParam) endif IParam = 4 + 15*(IIter-1) JParam = 9 + 15*(IIter-1) KParam = 14 + 15*(IIter-1) write(string,*),' Skewness / err :',StatVal(IParam), & StatVal(JParam) !, & StatVal(KParam) call WC (string) if(iiter.eq.niter) then write(11,*) StatVal(IParam) endif IParam = 5 + 15*(IIter-1) JParam = 10 + 15*(IIter-1) KParam = 15 + 15*(IIter-1) write(string,*),' Kurtosis / err :',StatVal(IParam), & StatVal(JParam) !, & StatVal(KParam) call WC (string) if(iiter.eq.niter) then write(11,*) StatVal(IParam) endif ENDDO close(11) deallocate(StatVal,STAT=stat) if (stat /= 0) stop 'stat_image, Desallocation Problem StatVal' Call Lafin('Stat_image') end program Stat_image ! =================================================================== SubRoutine calst (Nom,Nom2,iBP,StatVal,NIter,FClip,Eps,iMask,iPas) ! =================================================================== ! Calcul iteratif sur une image des statistiques elementaires qui ! suivent avec soit leur incertitude soit la valeur attendue pour ! une distribution mormale et la probabilite qui en decoule : ! - nombre de pixels concernes (fraction) ; ! - extrema de la population retenu ; ! - esperance mathematique (error) ; ! - deviation standard (rms error) ; ! - skewness de Fisher (vs. Gauss) ; ! - kurtosis de Fisher (vs. Gauss). ! ! iBP : nombre de bit par pixel ! StatVal : valeur des parametres statistiques ! NIter : nombre d'iterations maximal et utilise ! FClip : facteur de "clipping" sur la mediane ! Eps : critere de convergence sur l'ecart-type ! iPas : saut sur les pixels ! iMask : masquage de certains pixels ! =================================================================== parameter (NParam = 5) parameter (NIM = 50) ! 5(x3) parametres memorises ! au plus NIM = 50 iterations character Nom1*100,Comment*80 character (Len=*), Intent (inout) :: Nom,Nom2 DOUBLE PRECISION & dData,dData2,dDataMin,dDataMax, & dPixMin(NIM),dPixMax(NIM), & dNPix,dPixNb(NIM), & dS0,dS1,dS2,dS3,dS4, & dS01,dS02,dS03, & dM1,dM1b,dM1c,dM1d, & dCM2,dCM3,dCM4,& dKStat1,dKStat2,dKStat3,dKStat4, & dExpect(NIM),dStdDev(NIM),dFSkew(NIM),dFKurt(NIM), & dEExpec(NIM),dEStDev(NIM),dEFSkew(NIM),dEFKurt(NIM), & dERFC,dSkew,dPFSkew(NIM),dKurt,dPFKurt(NIM), & dFClip,dDelta,dEps ! real, dimension(:), Intent(inout) :: statval Real statval(1) real, dimension(:), allocatable :: rLigne1,rligne2 real, dimension(:,:), allocatable :: imData,rMask real Dx,Dy,Xd,Yd integer, dimension(:),allocatable :: IndC integer stat ! Nparam=5 ! Nim=50 NIter = max( 1,NIter ) ! garde-fous NIter = min( NIM,NIter ) FClip = max( 0.,FClip ) dFClip= dble( FClip ) dEps = dble( Eps ) if ( iMask.lt.0 .or. iMask.gt.1 ) then iMask = 0 endif if ( iPas.lt.1 ) then iPas = 1 endif ! =================================================================== Lu1 = 1 Lu2 = 2 ! --- ouverture de l'image --- Nom1 = Nom call OFits (Lu1,Nom1,iBP,Comment,Npl,Nl,Dx,Dy,Xd,Yd) allocate(rLigne1(Npl)) allocate(rLigne2(Npl)) allocate(imData(Npl,Nl)) allocate(rMask(Npl,Nl)) allocate(IndC(Npl*Nl)) ! --- saisie du masque --- if ( iMask.eq.1 ) then ! Nom1 doit etre les details de la W.T. call OFits (Lu2,Nom2,iBP,Comment,Npl,Nl,Dx,Dy,Xd,Yd) else DO Ipl = 1,Npl rLigne2(Ipl) = 0. ENDDO endif ! --- fenetrage --- Nl1 = 1 ! Fenetre sur les lignes Nl2 = Nl Nl3 = iPas Npl1 = 1 ! Fenetre sur les points Npl2 = Npl Npl3 = iPas NlS = (Nl2-Nl1)/Nl3 + 1 NplS = (Npl2-Npl1)/Npl3 + 1 iDim = 1 ! rangement mondimensionnel des valeurs ! === premiere estimation en considerant tous les pixels de la fenetre === IIter = 1 dPixMin(IIter) = 1.D9 dPixMax(IIter) = -1.D9 dS0 = 0.D0 ! Initialisation des "power sums" dS1 = 0.D0 dS2 = 0.D0 dS3 = 0.D0 dS4 = 0.D0 DO Il = Nl1,Nl2,Nl3 call RFits (Lu1,Il,Npl,iBP,rLigne1) if ( iMask.eq.1 ) then !call RFits (Lu2,Il,Npl,iBP,rLigne2) endif IlS = (Il-Nl1)/Nl3 + 1 DO Ipl = Npl1,Npl2,Npl3 IplS = (Ipl-Npl1)/Npl3 + 1 imData(IplS,IlS) = rLigne1(Ipl) rMask(IplS,IlS) = rLigne2(Ipl) if ( rMask(IplS,IlS).eq.0. ) then dData = dble( imData(IplS,IlS) ) dData2 = dData*dData dS0 = dS0 + 1.D0 dS1 = dS1 + dData dS2 = dS2 + dData2 dS3 = dS3 + dData2*dData dS4 = dS4 + dData2*dData2 dPixMin(IIter) = min( dData,dPixMin(IIter) ) dPixMax(IIter) = max( dData,dPixMax(IIter) ) endif ENDDO ENDDO call ftClos (Lu1,iStatus) call ftClos (Lu2,iStatus) dNPix = dS0 dPixNb(IIter) = dNPix dExpect(IIter) = 0.D0 dStdDev(IIter) = 0.D0 dFSkew(IIter) = 0.D0 dFKurt(IIter) = 0.D0 if ( dS0.ge.1.D0 ) then ! --- moments centres de l'echantillon --- dM1 = dS1/dS0 dM1b = dM1*dM1 dM1c = dM1*dM1b dM1d = dM1*dM1c dCM2 = dS2/dS0 - dM1b dCM3 = dS3/dS0 - 3.D0*dM1*dS2/dS0 + 2.D0*dM1c dCM4 = dS4/dS0 - 4.D0*dM1*dS3/dS0 + 6.D0*dM1b*dS2/dS0 & - 3.D0*dM1d ! --- estimateurs non biaises des moments centres de la population --- dS01 = dS0 - 1.D0 dS02 = dS0 - 2.D0 dS03 = dS0 - 3.D0 dKStat1 = dM1 if ( dS01.ge.1.D0 ) then dKStat2 = dCM2 * dS0/dS01 if ( dS02.ge.1.D0 ) then dKStat3 = dCM3 * dS0*dS0/( dS01*dS02 ) if ( dS03.ge.1.D0 ) then dKStat4 = (dS0+1.D0)*dCM4 - 3.D0*dS01*dCM2*dCM2 dKStat4 = dKStat4 * dS0*dS0/( dS01*dS02*dS03 ) else dKStat4 = 0.D0 endif else dKStat3 = 0.D0 endif else dKStat2 = 0.D0 endif ! --- descripteurs statistiques avec leur deviation standard --- dStdDev(IIter) = sqrt( dKStat2 ) ! deviation standard dEStDev(IIter) = 2.D0*dS0*dKStat2*dKStat2 + dS01*dKStat4 dEStDev(IIter) = dEStDev(IIter) / ( dS0*(dS0+1.D0) ) dEStDev(IIter) = dEStDev(IIter) / ( 4.D0*dKStat2 ) dEStDev(IIter) = sqrt( dEStDev(IIter) ) dExpect(IIter) = dKStat1 ! esperance mathematiquedEFSkew dEExpec(IIter) = dStdDev(IIter) / dsqrt( dS0 ) if ( dKStat2.ne.0.D0 .and. dS02.ge.1.D0 ) then dFSkew(IIter) = dKStat3 / dKStat2**1.5D0 ! facteur d'asymetrie dEFSkew(IIter) = dS02*(dS0+1.D0)*(dS0+3.D0) dEFSkew(IIter) = 6.D0*dS0*dS01 / dEFSkew(IIter) dEFSkew(IIter) = sqrt( dEFSkew(IIter) ) ! distribution normale dSkew = dFSkew(IIter) / dEFSkew(IIter) dPFSkew(IIter) = 0.5D0*dERFC( dSkew/dsqrt(2.D0) ) if ( dSkew.lt.0.D0 ) then dPFSkew(IIter) = 1.D0 - dPFSkew(IIter) endif else dEFSkew(IIter) = 0.D0 dPFSkew(IIter) = 0.D0 endif if ( dKStat2.ne.0.D0 .and. dS03.ge.1.D0 ) then dFKurt(IIter) = dKStat4 / dKStat2**2.0D0 ! facteur d'aplatissement dEFKurt(IIter) = dS02*dS03*(dS0+3.D0)*(dS0+5.D0) dEFKurt(IIter) = 24.D0*dS0*dS01*dS01 / dEFKurt(IIter) dEFKurt(IIter) = dsqrt( dEFKurt(IIter) ) ! distribution normale dKurt = dFKurt(IIter) / dEFKurt(IIter) dPFKurt(IIter) = 0.5D0*dERFC( dKurt/dsqrt(2.D0) ) if ( dKurt.lt.0.D0 ) then dPFKurt(IIter) = 1.D0 - dPFKurt(IIter) endif else dEFKurt(IIter) = 0.D0 dPFKurt(IIter) = 0.D0 endif else dEStDev(IIter) = 0.D0 dEFSkew(IIter) = 0.D0 dPFSkew(IIter) = 0.D0 dEFKurt(IIter) = 0.D0 dPFKurt(IIter) = 0.D0 endif ! === calcul avec rejet iteratif a concurence de la convergence === IIter = 2 dDelta = dEps + 1.D0 ! au moins une passe si NIter > 1 DO WHILE ( IIter.le.NIter .and. dDelta.gt.dEps .and. dS0.ne.0.D0 ) JIter = IIter - 1 dDataMin = dExpect(JIter) - dFClip*dStdDev(JIter) dDataMax = dExpect(JIter) + dFClip*dStdDev(JIter) dPixMin(IIter) = dDataMax dPixMax(IIter) = dDataMin dS0 = 0.D0 dS1 = 0.D0 dS2 = 0.D0 dS3 = 0.D0 dS4 = 0.D0 DO Il = 1,NlS DO Ipl = 1,NplS if ( rMask(Ipl,Il).eq.0 ) then dData = dble( imData(Ipl,Il) ) if ( dData.gt.dDataMin & .and. dData.le.dDataMax ) then dData2 = dData*dData dS0 = dS0 + 1.D0 dS1 = dS1 + dData dS2 = dS2 + dData2 dS3 = dS3 + dData2*dData dS4 = dS4 + dData2*dData2 dPixMin(IIter) = dmin1( dData,dPixMin(IIter) ) dPixMax(IIter) = dmax1( dData,dPixMax(IIter) ) endif endif ENDDO ENDDO dPixNb(IIter) = dS0 dExpect(IIter) = 0.D0 dStdDev(IIter) = 0.D0 dFSkew(IIter) = 0.D0 dFKurt(IIter) = 0.D0 if ( dS0.ge.1.D0 ) then dM1 = dS1/dS0 dM1b = dM1*dM1 dM1c = dM1*dM1b dM1d = dM1*dM1c dCM2 = dS2/dS0 - dM1b dCM3 = dS3/dS0 - 3.D0*dM1*dS2/dS0 + 2.D0*dM1c dCM4 = dS4/dS0 - 4.D0*dM1*dS3/dS0 + 6.D0*dM1b*dS2/dS0 & - 3.D0*dM1d dS01 = dS0 - 1.D0 dS02 = dS0 - 2.D0 dS03 = dS0 - 3.D0 dKStat1 = dM1 if ( dS01.ge.1.D0 ) then dKStat2 = dCM2 * dS0/dS01 if ( dS02.ge.1.D0 ) then dKStat3 = dCM3 * dS0*dS0/( dS01*dS02 ) if ( dS03.ge.1.D0 ) then dKStat4 = (dS0+1.D0)*dCM4 - 3.D0*dS01*dCM2*dCM2 dKStat4 = dKStat4 * dS0*dS0/( dS01*dS02*dS03 ) else dKStat4 = 0.D0 endif else dKStat3 = 0.D0 endif else dKStat2 = 0.D0 endif dStdDev(IIter) = dsqrt( dKStat2 ) dEStDev(IIter) = 2.D0*dS0*dKStat2*dKStat2 + dS01*dKStat4 dEStDev(IIter) = dEStDev(IIter) / ( dS0*(dS0+1.D0) ) dEStDev(IIter) = dEStDev(IIter) / ( 4.D0*dKStat2 ) dEStDev(IIter) = dsqrt( dEStDev(IIter) ) dExpect(IIter) = dKStat1 dEExpec(IIter) = dStdDev(IIter) / dsqrt( dS0 ) if ( dKStat2.ne.0.D0 .and. dS02.ge.1.D0 ) then dFSkew(IIter) = dKStat3 / dKStat2**1.5D0 dEFSkew(IIter) = dS02*(dS0+1.D0)*(dS0+3.D0) dEFSkew(IIter) = 6.D0*dS0*dS01 / dEFSkew(IIter) dEFSkew(IIter) = dsqrt( dEFSkew(IIter) ) dSkew = dFSkew(IIter) / dEFSkew(IIter) dPFSkew(IIter) = 0.5D0*dERFC( dSkew/dsqrt(2.D0) ) if ( dSkew.lt.0.D0 ) then dPFSkew(IIter) = 1.D0 - dPFSkew(IIter) endif else dEFSkew(IIter) = 0.D0 dPFSkew(IIter) = 0.D0 endif if ( dKStat2.ne.0.D0 .and. dS03.ge.1.D0 ) then dFKurt(IIter) = dKStat4 / dKStat2**2.0D0 dEFKurt(IIter) = dS02*dS03*(dS0+3.D0)*(dS0+5.D0) dEFKurt(IIter) = 24.D0*dS0*dS01*dS01 / dEFKurt(IIter) dEFKurt(IIter) = dsqrt( dEFKurt(IIter) ) dKurt = dFKurt(IIter) / dEFKurt(IIter) dPFKurt(IIter) = 0.5D0*dERFC( dKurt/dsqrt(2.D0) ) if ( dKurt.lt.0.D0 ) then dPFKurt(IIter) = 1.D0 - dPFKurt(IIter) endif else dEFKurt(IIter) = 0.D0 dPFKurt(IIter) = 0.D0 endif else dEStDev(IIter) = 0.D0 dEFSkew(IIter) = 0.D0 dPFSkew(IIter) = 0.D0 dEFKurt(IIter) = 0.D0 dPFKurt(IIter) = 0.D0 endif dDelta = (dStdDev(JIter)-dStdDev(IIter))/dStdDev(JIter) IIter = IIter + 1 ENDDO ! --- memorisation des evaluations successives --- NIter = IIter - 1 DO IIter = 1,NIter IParam = 1 + (IIter-1)*(3*NParam) StatVal(IParam) = sngl( dPixNb(IIter) ) IParam = IParam + 1 StatVal(IParam) = sngl( dExpect(IIter) ) IParam = IParam + 1 StatVal(IParam) = sngl( dStdDev(IIter) ) IParam = IParam + 1 StatVal(IParam) = sngl( dFSkew(IIter) ) IParam = IParam + 1 StatVal(IParam) = sngl( dFKurt(IIter) ) IParam = IParam + 1 StatVal(IParam) = sngl( 100.D0* dPixNb(IIter)/dNPix ) IParam = IParam + 1 StatVal(IParam) = sngl( dEExpec(IIter) ) IParam = IParam + 1 StatVal(IParam) = sngl( dEStDev(IIter) ) IParam = IParam + 1 StatVal(IParam) = sngl( dEFSkew(IIter) ) IParam = IParam + 1 StatVal(IParam) = sngl( dEFKurt(IIter) ) IParam = IParam + 1 dS0 = 1.D0- dPixNb(IIter)/dPixNb(max0(1,IIter-1)) StatVal(IParam) = sngl( 100.D0* dS0 ) IParam = IParam + 1 StatVal(IParam) = sngl( dPixMin(IIter) ) IParam = IParam + 1 StatVal(IParam) = sngl( dPixMax(IIter) ) IParam = IParam + 1 StatVal(IParam) = sngl( dPFSkew(IIter) ) IParam = IParam + 1 StatVal(IParam) = sngl( dPFKurt(IIter) ) ENDDO deallocate(rLigne1,rligne2,STAT=stat) if(stat /= 0) stop 'dans Stats, Pb desallocation rligne' deallocate(imData,rMask,STAT=stat) if(stat /= 0) stop 'dans Stats, Pb desallocation imData,rMask' deallocate(IndC,STAT=stat) if(stat /= 0) stop 'dans Stats, Pb desallocation IndC' ! =================================================================== end Subroutine calst DOUBLE PRECISION Function dERFC (dX) ! Calcul de la fonction d'erreur complementaire erfc(x) definie comme le ! produit de 2/sqrt(Pi) par l'integrale entre x et l'infini de exp(-t*t). ! Source : Numerical Recipes ! =================================================================== DOUBLE PRECISION & dX, & dT,dZ, & dA0,dA1,dA2,dA3,dA4,dA5,dA6,dA7,dA8,dA9 dA0 = -1.26551223D0 dA1 = 1.00002368D0 dA2 = 0.37409196D0 dA3 = 0.09678418D0 dA4 = -0.18628806D0 dA5 = 0.27886807D0 dA6 = -1.13520398D0 dA7 = 1.48851587D0 dA8 = -0.82215223D0 dA9 = 0.17087277D0 dZ = abs( dX ) dT = 1.D0 / ( 1.D0 + 0.5D0*dZ ) dERFC = dT * exp( -dZ*dZ + & dA0 + dT*( dA1 + dT*( dA2 + dT*( dA3 + dT*( & dA4 + dT*( dA5 + dT*( dA6 + dT*( dA7 + dT*( & dA8 + dT* dA9 )))))))) ) if ( dX < 0.D0 ) dERFC = 2.D0 - dERFC ! =================================================================== end function dERFC