C ++******************************************************************** C * C SSNRB.F * C OPFILEC FEB 03 ARDEAN LEITH C OPFILES REWRITE DEC 10 ARDEAN LEITH C * C ********************************************************************** C=* * C=* This file is part of: SPIDER - Modular Image Processing System. * C=* SPIDER System Authors: Joachim Frank & ArDean Leith * C=* Copyright 1985-2010 Health Research Inc., * C=* Riverview Center, 150 Broadway, Suite 560, Menands, NY 12204. * C=* Email: spider@wadsworth.org * C=* * C=* SPIDER is free software; you can redistribute it and/or * C=* modify it under the terms of the GNU General Public License as * C=* published by the Free Software Foundation; either version 2 of the * C=* License, or (at your option) any later version. * C=* * C=* SPIDER is distributed in the hope that it will be useful, * C=* but WITHOUT ANY WARRANTY; without even the implied warranty of * C=* merchantability or fitness for a particular purpose. See the GNU * C=* General Public License for more details. * C=* You should have received a copy of the GNU General Public License * C=* along with this program. If not, see * C=* * C ********************************************************************** C SSNRB() C C PARAMETERS: C C OPERATION: 'RF SN' C C PURPOSE: COMPUTE THE SPECTRAL SIGNAL-TO-NOISE RATIO (SSNR) OF A C SERIES OF IMAGES. C C23456789012345678901234567890123456789012345678901234567890123456789012 C*********************************************************************** SUBROUTINE SSNRB() INCLUDE 'CMBLOCK.INC' INCLUDE 'CMLIMIT.INC' CHARACTER(LEN=MAXNAM) :: FILPAT1,FILMASK DOUBLE PRECISION, ALLOCATABLE :: SUMCMP(:,:) DOUBLE COMPLEX, ALLOCATABLE :: FA(:,:) REAL, ALLOCATABLE :: B(:,:), RMSK(:,:) CHARACTER(LEN=MAXNAM) :: DOCNAM REAL :: DLIST(5) INTEGER, ALLOCATABLE :: LR(:) DOUBLE PRECISION, ALLOCATABLE :: SIGNAL(:), FR(:) CHARACTER (LEN=1) :: NULL = CHAR(0) CHARACTER (LEN=60) :: COMMENT LOGICAL :: NEWFILE LOGICAL :: FOUROK = .FALSE. LOGICAL :: ASKNAM = .TRUE. INTEGER, PARAMETER :: LUN1 = 20 INTEGER, PARAMETER :: LUNM = 21 INTEGER, PARAMETER :: LUNDOC = 80 INTEGER, PARAMETER :: LUNXM1 = 81 C OPEN FIRST INPUT FILE, NOT FOURIER MAXIM1 = 0 NILMAX = NIMAX ! INUMBR FROM CMLIMIT CALL OPFILES(0,LUN1,LUNDOC,LUNXM1, & ASKNAM,FILPAT1,NLET1, 'O', & IFORM1,NSAM1,NROW1,NSLICE1,MAXIM1, & NULL, & FOUROK, INUMBR,NILMAX, & NDUM,NGOT1,IMG1, IRTFLG) IF (IRTFLG .NE. 0) RETURN C OPEN MASK FILE (NOT FOURIER) MAXIMM = 0.0 CALL OPFILEC(0,ASKNAM,FILMASK,LUNM,'O',IFORMM, & NSAMM,NROWM,NSLICEM, & MAXIMM,'MASK',FOUROK,IRTFLG) IF (IRTFLG .NE. 0) RETURN CALL SIZCHK(UNUSED,NSAM1,NROW1,NSLICE1,0, & NSAMM,NROWM,NSLICEM,0,IRTFLG) IF (IRTFLG .NE. 0) GOTO 999 C TOTAL NUMBER OF IMAGES WRITE(NOUT,2001) NGOT1 2001 FORMAT(' NUMBER OF IMAGES: ',I5) WI = 1.0 CALL RDPRM1S(WI,NOT_USED,'RING WIDTH',IRTFLG) IF (IRTFLG .NE. 0) GOTO 999 Y1 = FLOAT(MAX(NSAM1,NROW1)) INC = INT(Y1/WI) / 2 + 1 NNNN = NSAM1 + 2 - MOD(NSAM1,2) ALLOCATE (RMSK(NSAM1,NROW1), & SUMCMP(NNNN/2,NROW1), & FA(NNNN/2,NROW1), & B(NNNN,NROW1), & LR(INC), & SIGNAL(INC), & FR(INC), STAT=IRTFLG) IF (IRTFLG .NE. 0) THEN MWANT = NSAM1*NROW1 + NNNN/2*NROW1 + NNNN/2*NROW1 + & NNNN*NROW1 + 3* INC CALL ERRT(46,'RF SN, RMSK',MWANT) GOTO 999 ENDIF C READ MASK CALL REDVOL(LUNM,NSAM1,NROW1,1,1,RMSK,IRTFLG) C INITIALIZE THE SUMS SUMCMP = 0.0D0 FA = CMPLX(0.0,0.0) NINDX1 = 1 DO ! LOOP OVER ALL IMAGES C READ CURRENT IMAGE CALL REDVOL(LUN1,NSAM1,NROW1,1,1,B,IRTFLG) C MASK MULTIPLY c$omp parallel do private(i,j) DO J=1,NROW1 DO I=1,NSAM1 B(I,J) = B(I,J) * RMSK(I,J) ENDDO ENDDO C FORWARD FFT INV = +1 CALL FMRS_2(B,NSAM1,NROW1,INV) IF (INV .NE. 1) GOTO 999 c$omp parallel do private(i,j,ix) DO J=1,NROW1 DO I=1,NNNN,2 IX = (I+1)/2 FA(IX,J) = FA(IX,J) + CMPLX(B(I,J),B(I+1,J)) SUMCMP(IX,J) = SUMCMP(IX,J) + B(I,J) * & DBLE(B(I,J))+B(I+1,J) * DBLE(B(I+1,J)) ENDDO ENDDO C OPEN NEXT INPUT FILE, UPDATES NINDX1 CALL NEXTFILE(NINDX1, INUMBR, & FOUROK, LUNXM1, & NGOT1, MAXIM1, & LUN1,0, FILPAT1, 'O', & IMG1, IRTFLG) IF (IRTFLG .EQ. -1) EXIT ! END OF INPUT STACK IF (IRTFLG .NE. 0) GOTO 999 ! ERROR CALL LUNGETSIZE(LUN1,NSAM,NROW,NSLICE,IRTFLG) CALL SIZCHK(UNUSED,NSAM1,NROW1, 1,0, & NSAM, NROW, NSLICE,0,IRTFLG) IF (IRTFLG .NE. 0) GOTO 999 ! BAD SIZE ENDDO SIGNAL = 0.0D0 FR = 0.0D0 LR = 0 DO 149 J=1,NROW1 JJ = J-1 IF (JJ .GT. NROW1/2) JJ = JJ-NROW1 DO 149 I=1,NNNN/2 IF (I.EQ.1 .AND. JJ.LT.0) GOTO 149 PII = 0.5*SQRT((FLOAT(JJ)/FLOAT(NROW1/2))**2 + & (FLOAT(I-1)/FLOAT(NSAM1/2))**2) IF (PII .LE. 0.5) THEN L = MIN(MAX(NINT(PII*Y1/WI)+1,1),INC) LR(L) = LR(L) + 1 SIGNAL(L) = SIGNAL(L) + ABS(FA(I,J))**2 FR(L) = FR(L) + SUMCMP(I,J) - & ABS(FA(I,J))**2 / NGOT1 ENDIF 149 CONTINUE ! DOES GOTO CYCLE INNER OR OUTER LOOP? al C SAVE RESULTS CALL OPENDOC(DOCNAM,.TRUE.,NLET,LUNDOC,NDOC,ASKNAM, & 'OUTPUT DOCUMENT',.FALSE.,.FALSE.,.TRUE., & NEWFILE,IRTFLG) IF (IRTFLG .NE. 0) GOTO 999 C 123456789 123456789 123456789 123456789 123456789 123456789 COMMENT = 'RADIUS NORM.-RAD. SSNR NRING VAR.' CALL LUNDOCPUTCOM(NDOC,COMMENT(1:51),IRTFLG) DO L=1,INC IF (LR(L).NE.0) THEN DLIST(1) = L DLIST(2) = FLOAT(L-1) / FLOAT(INC-1) * 0.5 DLIST(3) = AMAX1(0.0,SNGL((FLOAT(NGOT1-1)/NGTO1)* & SIGNAL(L) / FR(L))-1.0) DLIST(4) = LR(L) N1 = LR(L) N2 = (NGOT1-1) * LR(L) DLIST(5) = SQRT(2.*((N2+N1-2)*(1+2*DLIST(3)) + & N1 * DLIST(3)**2) / FLOAT(N1*(N2-4))) CALL LUNDOCWRTDAT(NDOC,L,DLIST(2),4,IRTFLG) !!CALL SAVD(NDOC,DLIST,NDLI,IRTFLG) ENDIF ENDDO 999 CLOSE(LUNDOC) CLOSE(LUNM) CLOSE(LUN1) IF (ALLOCATED(SUMCMP)) DEALLOCATE (SUMCMP) IF (ALLOCATED(FA)) DEALLOCATE (FA) IF (ALLOCATED(RMSK)) DEALLOCATE (RMSK) IF (ALLOCATED(B)) DEALLOCATE (B) END