C++******************************************************************** C C DSGR_PM.F C ADDED APSHIFT NOV 03 ARDEAN LEITH C AP_END FEB 04 ARDEAN LEITH C CROSRNG_E SPEEDS UP AUG 04 ARDEAN LEITH C AP_END CALL HAS PARLIST OCT 04 ARDEAN LEITH C CCROT STATISTICS FEB 05 ARDEAN LEITH C APRINGS_ONE MAR 05 ARDEAN LEITH C CROSRNG_E REWRITE MAR 08 ARDEAN LEITH C IREFLIST VAR. NAME MAR 08 ARDEAN LEITH C FMRS_PLAN MAR 08 ARDEAN LEITH C APRINGS_ONE_NEW MAR 08 ARDEAN LEITH C AP_END REDO NOV 08 ARDEAN LEITH C AP_STAT_ADD NOV 08 ARDEAN LEITH C ANGDIF EXPDIR BUG MAR 09 ARDEAN LEITH C CROSRNG_2, TT REMOVED JUN 10 ARDEAN LEITH C AP_STAT NBORDER OCT 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 C DSGR_PM C C PURPOSE: FIND ROTATIONAL AND SHIFT PARAMETERS TO ALIGN A SERIES OF C REFERENCE IMAGES WITH SAMPLE IMAGES 'AP REF' C C ORDER OF PROCESSING: C 1 - CONVERT EACH REFERENCE IMAGE TO RINGS, DO THE FFTS C FOR ALL THE RINGS, APPLY WEIGHTS TO THE RINGS, STORE IT IN CIRCREF. C IN ADDITION, HIGHEST FREQUENCY FOR ALL THE RINGS EXCEPT C MAXRING ARE DIVIDED BY 2. C 2 - CONVERT EACH INPUT IMAGE TO RINGS, DO THE FFTS, C COMPARE EACH INPUT IMAGE WITH ALL THE REFERENCE IMAGES C USING CROSRNG_E. C SINCE Y WERE PRE-WEIGHTED THE RESULT IS ALREADY CORRECT. C C PARAMETERS: C IREFLIST LIST OF REF. IMAGE FILE NUMBERS (INPUT) C NUMREF NO. OF IMAGES (INPUT) C IEXPLIST LIST OF EXP. IMAGE FILE NUMBERS (INPUT) C NUMEXP NO. OF IMAGES (INPUT) C NSAM,NROW ACTUAL (NOT-WINDOWED) IMAGE SIZE (INPUT) C REFANGDOC REF. ANGLES FILE NAME (INPUT) C EXPANGDOC EXP. ANGLES FILE NAME (INPUT) C REFPAT REF. IMAGE SERIES FILE TEMPLATE (INPUT) C EXPPAT EXP. IMAGE SERIES FILE TEMPLATE (INPUT) C C OPERATIONS: 'AP MD', 'AP REF', 'AP RD', 'AP RN' C USED WHEN: (CIRCREF_IN_CORE && NUMTH.GT.1 && NUMEXP.GT.NUMTH) .OR. C CTYPE .EQ. 'T' ) C C23456789 123456789 123456789 123456789 123456789 123456789 123456789 12 C--********************************************************************* SUBROUTINE DSGR_PM(IREFLIST,NUMREF,IEXPLIST,NUMEXP, & NSAM,NROW,ANGDIFTHR,RANGE, & NRING,LCIRC,NUMR,CIRCREF, & MODE, REFANGDOC,EXPANGDOC,SCRFILE,FFTW_PLANS, & REFPAT,EXPPAT,CKMIRROR,CTYPE,ISHRANGE,LUNDOC) INCLUDE 'CMLIMIT.INC' INCLUDE 'CMBLOCK.INC' INTEGER, DIMENSION(NUMREF) :: IREFLIST INTEGER, DIMENSION(NUMEXP) :: IEXPLIST INTEGER, DIMENSION(3,NRING) :: NUMR REAL, DIMENSION(LCIRC,NUMREF) :: CIRCREF CHARACTER(LEN=1) :: MODE,NULL CHARACTER (LEN=*) :: REFPAT,EXPPAT,SCRFILE CHARACTER (LEN=*) :: REFANGDOC,EXPANGDOC LOGICAL :: CKMIRROR CHARACTER (LEN=*) :: CTYPE C ALLOCATABLE ARRAYS REAL, ALLOCATABLE, DIMENSION(:,:,:) :: EXPBUF REAL, ALLOCATABLE, DIMENSION(:) :: AVI,SIGI REAL, ALLOCATABLE, DIMENSION(:) :: CCLIST REAL, ALLOCATABLE, DIMENSION(:) :: RANGNLIST INTEGER, ALLOCATABLE, DIMENSION(:) :: NPROJ INTEGER, ALLOCATABLE, DIMENSION(:) :: ILIST LOGICAL, ALLOCATABLE, DIMENSION(:) :: MIRLIST !integer, allocatable, dimension(:) :: iptr INTEGER *8 :: IPLAN !STRUCTURE POIN LOGICAL :: CIRCREF_IN_CORE LOGICAL :: MIRRORNEW LOGICAL :: GOTREFANG C AUTOMATIC ARRAYS REAL, DIMENSION(6) :: DLIST REAL, DIMENSION(8,NUMEXP) :: ANGEXP REAL, DIMENSION(3,NUMEXP) :: EXPDIR REAL, DIMENSION(3,NUMREF) :: REFDIR,ANGREF INTEGER, PARAMETER :: NLISTMAX = 15 REAL, DIMENSION(NLISTMAX) :: PARLIST REAL, PARAMETER :: QUADPI = 3.14159265358979323846 REAL, PARAMETER :: DGR_TO_RAD = (QUADPI/180) INTEGER :: NBORDER = 0 ! UNUSED INTEGER :: NSUBPIX = 0 ! UNUSED DATA INPIC/77/,INANG/78/,LUNRING/50/ NULL = CHAR(0) C FIND NUMBER OF OMP THREADS CALL GETTHREADS(NUMTH) #ifdef DEBUGD NR = NUMR(1,NRING) ! OUTER RING NUMBER MAXRIN = NUMR(3,NRING) write(88,*) ' DEBUG OUTPUT FROM DSGR_PM' write(88,902) NR,MAXRIN 902 format(' OUTER RING NUMBER: ',I10,' RING LENGTH:',I10) do i = 1, nring write(88,901) NUMR(1,I),NUMR(2,I),NUMR(3,I) 901 format(3i10) enddo #endif C THIS ALTERS RANGE! RANGE = COS(RANGE*DGR_TO_RAD) C MAKE EXPBUF BIG ENOUGH FOR ALL THREADS MWANTX = NSAM * NROW * NUMTH ALLOCATE(EXPBUF(NSAM,NROW,NUMTH), & AVI(NUMTH), & SIGI(NUMTH), & ILIST(NUMTH), & CCLIST(NUMTH), & RANGNLIST(NUMTH), & NPROJ(NUMTH), & MIRLIST(NUMTH), & STAT=IRTFLG) IF (IRTFLG .NE. 0) THEN MWANT = MWANTX + 7*NUMTH CALL ERRT(46,'EXPBUF...',MWANT) GOTO 9999 ENDIF C DUMMY CALL TO INITIALIZE REV. FFTW3 PLAN FOR USE WITHIN OMP || CALL FMRS_PLAN(.FALSE.,DUM,NSAM,NROW,1, 1, -1, IPLAN, IRTFLG) CALL FMRS_PLAN(.FALSE.,DUM,NSAM,NROW,1, 1, +1, IPLAN, IRTFLG) C ZERO DLIST ARRAY DLIST = 0.0 GOTREFANG = .FALSE. IF (CTYPE(1:2) .NE. 'MD') THEN C GET REF. PROJ. ANGLES (ANGREF) FROM DOC. FILE (REFANGDOC) OR C IMAGE FILE (REFPAT) HEAD CALL AP_GETANGA(IREFLIST,NUMREF,0,REFANGDOC,REFPAT, & INPIC,INANG,3,ANGREF,NGOTREF,IRTFLG) IF (IRTFLG .NE. 0) GOTO 9999 GOTREFANG = .TRUE. C CONVERT REF. ANGLES TO UNITARY DIRECTIONAL VECTORS (REFDIR). CALL AP_GETSATA(ANGREF,REFDIR,3,NUMREF,IRTFLG) IF (IRTFLG .NE. 0) GOTO 9999 ENDIF NGOTPAR = 0 IF (EXPANGDOC .NE. NULL) THEN C LOAD EXP. PROJ. ANGLES & ALIGNMENT PARAMETERS (ANGEXP) C FROM DOC. FILE (EXPANGDOC) CALL AP_GETANGA(IEXPLIST,NUMEXP,0,EXPANGDOC,EXPPAT, & INPIC,INANG,8,ANGEXP,NGOTPAR,IRTFLG) IF (IRTFLG .NE. 0) GOTO 9999 C CONVERT EXP. ANGLES TO UNITARY DIRECTIONAL VECTORS (EXPDIR). CALL AP_GETSATA(ANGEXP,EXPDIR,8,NUMEXP,IRTFLG) IF (IRTFLG .NE. 0) GOTO 9999 ELSE ANGEXP = 0.0 ENDIF C READ REFERENCE IMAGES INTO REFERENCE RINGS (CIRCREF) ARRAY CIRCREF_IN_CORE = .TRUE. CALL APRINGS_NEW(IREFLIST,NUMREF, NSAM,NROW, & NRING,LCIRC,NUMR,MODE,FFTW_PLANS, & REFPAT,INPIC,CIRCREF,CIRCREF_IN_CORE, & LUNRING,SCRFILE,IRTFLG) IF (IRTFLG .NE. 0) GOTO 9999 C INITIALIZE CCROT STATISTICS CALL AP_STAT_ADD(-1,CCROT,ANGDIF,ANGDIFTHR,CCROTLAS, & CCROTAVG,IBIGANGDIF,ANGDIFAVG,IMPROVCCROT, & CCROTIMPROV,IWORSECCROT,CCROTWORSE) DO IEXPT=1,NUMEXP,NUMTH C LOOP OVER ALL EXPERIMENTAL (SAMPLE) IMAGES C LOAD NUMTH EXPERIMENTAL IMAGES INTO ARRAY EXPBUF IEND = MIN(NUMEXP,IEXPT+NUMTH-1) CALL AP_GETDATS(IEXPLIST,NUMEXP,NSAM,NROW, & NUMTH,EXPPAT,INPIC, IEXPT,IEND, & EXPBUF,AVI,SIGI, IRTFLG) IF (IRTFLG .NE. 0) GOTO 9999 C FIND CLOSEST MATCHING REFERENCE IMAGE FOR EACH EXP. IMAGE c$omp parallel do private(IEXP,IT) DO IEXP=IEXPT,IEND IT = IEXP-IEXPT+1 CALL DSGR_2(EXPBUF(1,1,IT),NSAM,NROW,NUMR,NRING, & MODE,CIRCREF,LCIRC,NUMREF, & REFDIR,EXPDIR(1,IEXP),RANGE, & ILIST(IT),CCLIST(IT),RANGNLIST(IT),MIRLIST(IT), & NPROJ(IT),CKMIRROR,FFTW_PLANS) ENDDO C LOOP OVER CURRENT SET OF EXPERIMENTAL (SAMPLE) IMAGES DO IEXP=IEXPT,MIN(NUMEXP,IEXPT+NUMTH-1) C IT POINTS TO CURRENT EXP. IMAGE IN THE SUBLIST IT = IEXP - IEXPT + 1 C ILIST(IT) IS LIST NUMBER OF MOST SIMILAR REF. IMAGE IREF = ILIST(IT) IF (IREF .LE. 0) THEN C NO NEARBY REFERENCE IMAGE IMGREF = 0 C IREFT IS FOR REFDIR INDEX IREFT = 1 ELSE IMGREF = IREFLIST(IREF) C IREFT IS FOR REFDIR INDEX IREFT = IREF ENDIF IMGEXP = IEXPLIST(IEXP) CCROT = CCLIST(IT) RANGNEW = RANGNLIST(IT) MIRRORNEW = MIRLIST(IT) PEAKV = 0.0 XSHNEW = 0.0 YSHNEW = 0.0 IF (IMGREF .GT. 0 .AND. ISHRANGE .GT. 0) THEN C DETERMINE SHIFT PARAMETERS FOR EXP. IMAGE IN EXPBUF NSAMP = 2*NSAM+2 NROWP = 2*NROW CALL APSHIFT(INPIC, REFPAT,IMGREF, & NSAM,NROW, NSAMP,NROWP, & EXPBUF(1,1,IT),AVI(IT),SIGI(IT), ISHRANGE, & RANGNEW,XSHNEW,YSHNEW,MIRRORNEW,PEAKV,IRTFLG) IF (IRTFLG .NE. 0) RETURN ENDIF C AP_END WRITES ALIGNMENT PARAMETERS TO DOC FILE CALL AP_END(IEXP,IMGEXP,IMGREF, & ANGREF(1,IREFT),REFDIR(1,IREFT), & ANGEXP(1,IEXP),EXPDIR(1,IEXP),ISHRANGE, & GOTREFANG, NGOTPAR, CCROT,PEAKV, & RANGNEW,XSHNEW,YSHNEW,MIRRORNEW,REFPAT, & NPROJ(IT), CTYPE, LUNDOC,PARLIST) CALL AP_STAT_ADD(NGOTPAR,CCROT,PARLIST(10), & ANGDIFTHR,ANGEXP(8,IEXP), & CCROTAVG,IBIGANGDIF,ANGDIFAVG,IMPROVCCROT, & CCROTIMPROV,IWORSECCROT,CCROTWORSE) ENDDO ENDDO IF (LUNDOC .GT. 0) THEN C SAVE CCROT & ANGULAR DISPLACEMENT STATISTICS CALL AP_STAT(NUMEXP,ANGDIFTHR,IBIGANGDIF, & ANGDIFAVG, CCROTAVG, & IMPROVCCROT,CCROTIMPROV, & IWORSECCROT,CCROTWORSE, & NBORDER,NSUBPIX,LUNDOC) ENDIF 9999 IF (ALLOCATED(EXPBUF)) DEALLOCATE(EXPBUF) IF (ALLOCATED(AVI)) DEALLOCATE(AVI) IF (ALLOCATED(SIGI)) DEALLOCATE(SIGI) IF (ALLOCATED(ILIST)) DEALLOCATE(ILIST) IF (ALLOCATED(CCLIST)) DEALLOCATE(CCLIST) IF (ALLOCATED(RANGNLIST)) DEALLOCATE(RANGNLIST) IF (ALLOCATED(MIRLIST)) DEALLOCATE(MIRLIST) IF (ALLOCATED(NPROJ)) DEALLOCATE(NPROJ) END C++************************** DSGR_2 ********************************** SUBROUTINE DSGR_2(XIM,NSAM,NROW, NUMR,NRING,MODE, & CIRCREF,LCIRC,NUMREF, & SA,TA,RANGE, & IMGREFL,CCROT,RANGNEW,MIRNEW,NPROJ, & CKMIRROR,FFTW_PLANS) C NOTE: RUNS WITHIN OMP PARALLEL SECTION OF CODE! REAL, DIMENSION(NSAM,NROW) :: XIM INTEGER, DIMENSION(3,NRING) :: NUMR CHARACTER(LEN=1) :: MODE REAL, DIMENSION(LCIRC,NUMREF) :: CIRCREF REAL, DIMENSION(3,NUMREF) :: SA REAL, DIMENSION(3) :: TA REAL :: RANGE LOGICAL :: CKMIRROR,ONLYMIRROR LOGICAL :: MIRNEW,LIMITRANGE LOGICAL :: ISMIRRORED LOGICAL :: USE_UN,USE_MIR INTEGER *8 :: FFTW_PLANS(*) DOUBLE PRECISION :: CCROTD,TOTMIN C AUTOMATIC ARRAYS INTEGER, DIMENSION(NUMREF) :: LCG REAL, DIMENSION(LCIRC) :: CIRCEXP MAXRIN = NUMR(3,NRING) #ifdef SP_LIBFFTW3 MAXRIN = NUMR(3,NRING) - 2 #endif WR = 0.0 NPROJ = NUMREF LIMITRANGE = (RANGE .LT. 1.0) IF (LIMITRANGE) THEN C DETERMINE WHICH ONES ARE TO BE COMPARED NPROJ = 0 DO IMI=1,NUMREF C ABS - DIRECTIONS AT 180 DEGREES ARE THE SAME (MIRRORED) C DT NEAR 1.0 = NOT-MIRRORED, DT NEAR -1.0 = MIRRORED DT = (TA(1) * SA(1,IMI) + & TA(2) * SA(2,IMI) + & TA(3) * SA(3,IMI)) DTABS = ABS(DT) IF (DTABS .GE. RANGE) THEN C EITHER MIRRORED OR NON-MIRRORED IS WITHIN RANGE NPROJ = NPROJ + 1 LCG(NPROJ) = IMI IF (DT .LT. 0) LCG(NPROJ) = -IMI ENDIF ENDDO IF (NPROJ .LE. 0) THEN C NO REF. IMAGE WITHIN COMPARISON ANGLE IMGREFL = 0 CCROT = -1.0 RANGNEW = 0.0 MIRNEW = .FALSE. RETURN ENDIF ENDIF C CALCULATE DIMENSIONS FOR NORMALIZING IN APRINGS_ONE CNS2 = NSAM/2+1 CNR2 = NROW/2+1 C EXTRACT EXP. IMAGE RINGS, NORMALIZE & FFT THEM CALL APRINGS_ONE_NEW(NSAM,NROW, CNS2,CNR2, XIM, .FALSE., & MODE,NUMR,NRING,LCIRC, WR,FFTW_PLANS, & CIRCEXP,IRTFLG) CCROTD = -1.0D20 IMGREFL = 0 DO IMIL=1,NPROJ IMI = IMIL IF (LIMITRANGE) IMI = ABS(LCG(IMIL)) IF (CKMIRROR .AND. LIMITRANGE) THEN C ONLY SEARCH EITHER MIRRORED OR NON-MIRRORED USE_UN = (LCG(IMIL) .GE. 0) USE_MIR = (LCG(IMIL) .LT. 0) ELSE C SEARCH BOTH MIRRORED & NON-MIRRORED IF CHKMIR USE_UN = .TRUE. USE_MIR = CKMIRROR ENDIF C CHECK MIRRORED/ NON-MIRRORED POSITION CALL CROSRNG_2(CIRCREF(1,IMI),CIRCEXP, & LCIRC,NRING, MAXRIN,NUMR, & .FALSE.,FFTW_PLANS(1), & USE_UN,USE_MIR, & ISMIRRORED, TOTMIN,TOT) IF (TOTMIN .GE. CCROTD) THEN C GOOD MATCH WITH TOTA (MIRRORED OR NOT) POSITION CCROTD = TOTMIN RANGNEW = TOT MIRNEW = ISMIRRORED IMGREFL = IMI ENDIF ENDDO !END OF: DO IMIL=1,NPROJ IF (MODE .EQ. 'F') THEN RANGNEW = (RANGNEW-1.0) / MAXRIN*360.0 ELSE RANGNEW = (RANGNEW-1.0) / MAXRIN*180.0 ENDIF CCROT = CCROTD END