C++********************************************************************* C C RCTFSS.F ADDED 2D OPERATION NOV 2000 HAIXIAO C INPIC BUG APR 02 ARDEAN LEITH C OPFILEC FEB 03 ARDEAN LEITH C ILIST FIXED MAY 09 ARDEAN LEITH C ********************************************************************** C=* FROM: SPIDER - MODULAR IMAGE PROCESSING SYSTEM. AUTHOR: J.FRANK * C=* Copyright (C) 1985-2009 Health Research Inc. * C=* * C=* HEALTH RESEARCH INCORPORATED (HRI), * C=* ONE UNIVERSITY PLACE, RENSSELAER, NY 12144-3455. * C=* * C=* Email: spider@wadsworth.org * C=* * C=* This program 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=* This program 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=* * C=* You should have received a copy of the GNU General Public License * C=* along with this program; if not, write to the * C=* Free Software Foundation, Inc., * C=* 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. * C=* * C ********************************************************************** C C RCTFSS C C PURPOSE: 2D and 3D CTF CORRECTION USING WIENER FILTERING APPROACH C C F(K)=SUM(H I(K)FI(K))/(SUM(HI(K)^2)+1/SNR), WHERE HI(K) IS C CTF FUNCTION FOR THE ITH IMAGE; C FI(K) IS THE ITH IMAGE; SNR IS SIGNAL-TO-NOISE RATION. C C23456789012345678901234567890123456789012345678901234567890123456789012 C--********************************************************************* SUBROUTINE RCTFSS INCLUDE 'CMBLOCK.INC' INCLUDE 'CMLIMIT.INC' CHARACTER(LEN=MAXNAM) :: FINPAT,FILNAM,FINPIC INTEGER, ALLOCATABLE, DIMENSION(:) :: ILIST COMPLEX, ALLOCATABLE, DIMENSION(:,:,:) :: VOLIN,VOLSUM,CTF2 CHARACTER(LEN=1) :: NULL LOGICAL :: FIRST DATA LUN1,INPIC/10,55/ CALL SET_MPI(ICOMM,MYPID,MPIERR) ! SET MYPID NILMAX = NIMAX ALLOCATE(ILIST(NILMAX),STAT=IRTFLG) IF (IRTFLG .NE. 0) THEN CALL ERRT(46,'TF CTS, ILIST',NILMAX) RETURN ENDIF C ASK FOR DATA FILE CALL FILELIST(.TRUE.,INPIC,FINPAT,NLET,ILIST,NILMAX,NIMA, & 'ENTER TEMPLATE FOR IMAGE FILE',IRTFLG) IF (IRTFLG .NE. 0) GOTO 9999 NULL = CHAR(0) CALL FILERD(FILNAM,NLEP,NULL,'ENTER TEMPLATE FOR CTF',IRTFLG) IF (IRTFLG .NE. 0) GOTO 9999 CALL RDPRM(SNR,NOT_USED,'SNR') IF (SNR .LE. 0.0) THEN CALL ERRT(31,'TF CTF',IER) GOTO 9999 ENDIF S = 1.0/SNR DO K=1,NIMA CALL FILGET(FINPAT,FINPIC,NLET,ILIST(K),IRTFLG) IF (IRTFLG .NE. 0) GOTO 9999 MAXIM = 0 CALL OPFILEC(0,.FALSE.,FINPIC,LUN1,'O',IFORM, & NSAM1,NROW1,NSLICE1, & MAXIM,'DUMMY',.TRUE.,IRTFLG) IF (IRTFLG.NE.0) GOTO 9999 C PROCESS THE VOLUME INUMDEM = 0 CC USE INUMDEM TO RECORD 3D OR 2D CC INUMDEM = 3 3D CC 1 2D IF (IFORM.EQ.-11) THEN LRCL = NSAM1 LS = NSAM1 NSAM1 = NSAM1-1 INUMDEM = 1 ELSEIF (IFORM.EQ.-12) THEN LRCL = NSAM1 LS = NSAM1 NSAM1 = NSAM1-2 INUMDEM = 1 ELSEIF (IFORM .EQ. 1) THEN LS =NSAM1+2-MOD(NSAM1,2) LRCL =NSAM1 INUMDEM = 1 ELSEIF (IFORM .EQ. -21) THEN LRCL = NSAM1 LS = NSAM1 NSAM1 = NSAM1-1 INUMDEM = 3 ELSEIF (IFORM .EQ. -22) THEN LRCL = NSAM1 LS = NSAM1 NSAM1 = NSAM1-2 INUMDEM = 3 ELSEIF (IFORM .EQ. 3) THEN LS = NSAM1+2-MOD(NSAM1,2) LRCL = NSAM1 INUMDEM = 3 ELSE IF (MYPID .LE. 0) CLOSE(LUN1) CALL ERRT(2,'TF CTF',NE) GOTO 9999 ENDIF C GET THE DIMENSION INFORMATION FOR THE FIRST TIME IF (K .EQ. 1) THEN NSAM = NSAM1 NROW = NROW1 NSLICE = NSLICE1 C VOLUMES NEEDED ALLOCATE(VOLIN(LS/2,NROW,NSLICE),STAT=IRTFLG) IF (IRTFLG.NE.0) CALL ERRT(46,'TF CTF, VOLIN',IER) ALLOCATE(VOLSUM(LS/2,NROW,NSLICE),STAT=IRTFLG) IF (IRTFLG.NE.0) CALL ERRT(46,'TF CTF, VOLSUM',IER) ALLOCATE(CTF2(LS/2,NROW,NSLICE),STAT=IRTFLG) IF (IRTFLG.NE.0) CALL ERRT(46,'TF CTF, CTF2',IER) C SET OUTPUT VOLUME AND CTF^2 TO ZERO VOLSUM = 0.0 CTF2 = 0.0 ELSE C C VERIFY THE DIMENSIONS IF(NSAM.NE.NSAM1.OR.NROW.NE.NROW1.OR.NSLICE.NE.NSLICE1) THEN IF (MYPID .LE. 0) CLOSE(LUN1) CALL ERRT(1,'TF CTF',NE) GOTO 9999 ENDIF ENDIF C READ INPUT FILE DO L=1,NSLICE DO J=1,NROW CALL REDLIN(LUN1,VOLIN(1,J,L),LRCL,J+(L-1)*NROW) ENDDO ENDDO IF (MYPID .LE. 0) CLOSE(LUN1) IF ((IFORM.EQ.1) .OR. (IFORM.EQ.3)) THEN C FOURIER TRANSFORM INV = +1 CALL FMRS_3(VOLIN,NSAM,NROW,NSLICE,INV) IF (INV .EQ. 0) THEN CALL ERRT(38,'TF CTF',NE) GOTO 9999 ENDIF ENDIF C GET CTF CALL FILGET(FILNAM,FINPIC,NLEP,ILIST(K),IRTFLG) IF (IRTFLG .NE. 0) GOTO 9999 MAXIM=0 CALL OPFILEC(0,.FALSE.,FINPIC,LUN1,'O',IFORM, & NSAM1,NROW1,NSLICE1, & MAXIM,'DUMMY',.TRUE.,IRTFLG) IF (IRTFLG.NE.0) GOTO 9999 LRCL1 = NSAM1 C ONLY 2D FOURIER INUMCTF = 0 CC USE INUMCTF TO RECORD 3D OR 2D CC INUMCTF = 3 3D C 1 2D IF (IFORM .EQ. -11) THEN NSAM1 = NSAM1-1 INUMCTF = 1 ELSEIF (IFORM .EQ. -12) THEN NSAM1 = NSAM1-2 INUMCTF = 1 ELSEIF (IFORM .EQ. -21) THEN NSAM1 = NSAM1-1 INUMCTF = 3 ELSEIF (IFORM .EQ. -22) THEN NSAM1 = NSAM1-2 INUMCTF = 3 ELSE IF (MYPID .LE. 0) CLOSE(LUN1) CALL ERRT(2,'TF CTF',NE) GOTO 9999 ENDIF IF (INUMCTF .NE. INUMDEM) THEN IF (MYPID .LE. 0) CLOSE(LUN1) CALL ERRT(2,'TF CTF',NE) GOTO 9999 ENDIF c VERIFY DIMENSIONS OF THE CTF IF (NSAM.NE.NSAM1.OR.NROW.NE.NROW1.OR.NSLICE.NE.NSLICE1) THEN IF (MYPID .LE. 0) CLOSE(LUN1) CALL ERRT(1,'TF CTF',NE) GOTO 9999 ENDIF C MUTLIPLY INPUT VOLUME STORED IN THE MEMORY volin, PLACE OUTPUT C IN VOLSUM AND SQUARED CTF IN CTF2. CALL MULFC3(LUN1,VOLIN,VOLSUM,CTF2,LRCL1/2,NROW,NSLICE) IF (MYPID .LE. 0) CLOSE(LUN1) C END OF DO-LOOP OVER ALL VOLUMES ENDDO C DIVIDE BY THE SUM OF H^2+1/SNR VOLIN = VOLSUM/(REAL(CTF2)+S) c CALL DIVCR3(B(KAD2),B(KAD2),B(KAD3),S,LRCL1/2,NROW,NSLICE) INV=-1 CALL FMRS_3(VOLIN,NSAM,NROW,NSLICE,INV) IFORM = INUMCTF MAXIM = 0 CALL OPFILEC(0,.TRUE.,FILNAM,LUN1,'N',IFORM,NSAM,NROW,NSLICE, & MAXIM,'OUTPUT',.FALSE.,IRTFLG) IF (IRTFLG.NE.0) GOTO 9999 C WRITE THE OUTPUT DO L=1,NSLICE DO J=1,NROW CALL WRTLIN(LUN1,VOLIN(1,J,L),NSAM,J+(L-1)*NROW) ENDDO ENDDO IF (MYPID .LE. 0) CLOSE(LUN1) 9999 IF(ALLOCATED(CTF2)) DEALLOCATE(CTF2) IF(ALLOCATED(VOLSUM)) DEALLOCATE(VOLSUM) IF(ALLOCATED(VOLIN)) DEALLOCATE(VOLIN) IF(ALLOCATED(ILIST)) DEALLOCATE(ILIST) END