C++********************************************************************* C C DSGR_P.F C CROSRNG_E SPEEDS UP AUG 04 ARDEAN LEITH C AP_END CALL HAS PARLIST OCT 04 ARDEAN LEITH C SPIDER REF_RINGS FILE FEB 05 ARDEAN LEITH C CCROT STATISTICS FEB 05 ARDEAN LEITH C NON-INCORE PARALLEL FEB 05 ARDEAN LEITH C APRINGS_ONE MAR 05 ARDEAN LEITH C FFTW_PLANS MAR 08 ARDEAN LEITH C APRINGS_ONE_NEW MAR 08 ARDEAN LEITH C MPI PARTAB BUG FIXED OCT 08 ARDEAN LEITH C MPI LOCAL READ NOW NOV 08 ARDEAN LEITH C AP_STAT_ADD NOV 08 ARDEAN LEITH C C ********************************************************************** C=* FROM: SPIDER - MODULAR IMAGE PROCESSING SYSTEM. AUTHOR: J.FRANK * C=* Copyright (C) 1985-2008 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 DSGR_P(IREFLIST,NUMREF,IEXPLIST,NUMEXP, C NSAM,NROW,NR,LENTT,RANGE, C NRING,LCIRC,NUMR,CIRCREF,CIRCREF_IN_CORE, C MODE,REFANGDOC,EXPANGDOC,SCRFILE,FFTW_PLANS, C REFPAT,EXPPAT,CKMIRROR,CTYPE,ISHRANGE,LUNDOC) C C PURPOSE: FIND ROTATIONAL AND SHIFT PARAMETERS TO ALIGN A SERIES OF C REFERENCE IMAGES WITH SAMPLE IMAGES 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_**. C SINCE Y WERE PRE-WEIGHTED THE RESULT ARE 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 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 REF', 'AP RD', 'AP RN' C C--********************************************************************* #ifndef USE_MPI C --------------------- NON-MPI CODE -------------------------------- SUBROUTINE DSGR_P(IREFLIST,NUMREF, IEXPLIST,NUMEXP, & NSAM,NROW,LENTT,RANGE,ANGDIFTHR, & NRING,LCIRC,NUMR,CIRCREF,CIRCREF_IN_CORE, & 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 CHARACTER (LEN=*) :: REFANGDOC CHARACTER (LEN=*) :: EXPANGDOC CHARACTER (LEN=*) :: SCRFILE CHARACTER (LEN=*) :: REFPAT,EXPPAT CHARACTER (LEN=*) :: CTYPE CHARACTER (LEN=74) :: COMMENT DOUBLE PRECISION :: CCROTD CHARACTER(LEN=1) :: NULL LOGICAL :: CIRCREF_IN_CORE LOGICAL :: CKMIRROR LOGICAL :: MIRRORNEW LOGICAL :: GOTREFANG LOGICAL :: ONLYMIRROR LOGICAL :: LIMITRANGE LOGICAL :: MIRRORED INTEGER *8 :: FFTW_PLANS(*) C ALLOCATABLE ARRAYS REAL, ALLOCATABLE, DIMENSION(:) :: EXPBUF C AUTOMATIC ARRAYS DOUBLE PRECISION, DIMENSION(NUMREF) :: TOTMIN,TOTMIR DOUBLE PRECISION, DIMENSION(LENTT) :: TT REAL, DIMENSION(NUMREF) :: ROTANGT,ROTANGM REAL, DIMENSION(LCIRC) :: CIRCEXP REAL, DIMENSION(6) :: DLIST REAL, DIMENSION(3) :: EXPDIR REAL, DIMENSION(8,NUMEXP) :: ANGEXP REAL, DIMENSION(3,NUMREF) :: ANGREF,REFDIR INTEGER, DIMENSION(NUMREF) :: LCG INTEGER :: NSAID = 0 INTEGER, PARAMETER :: NLISTMAX = 15 REAL, DIMENSION(NLISTMAX) :: PARLIST PARAMETER (QUADPI = 3.1415926535897932384626) PARAMETER (DGR_TO_RAD = (QUADPI/180)) DATA INPIC,LUNANG,LUNRING/77,78,50/ MYPID = -1 ! NOT USING MPI NULL = CHAR(0) MAXRIN = NUMR(3,NRING) ! ACTUAL LENGTH OF LONGEST RING #ifdef SP_LIBFFTW3 MAXRIN = NUMR(3,NRING) - 2 #endif C THIS ALTERS RANGE! RANGE = COS(RANGE*DGR_TO_RAD) C FIND NUMBER OF OMP THREADS CALL GETTHREADS(NUMTH) IF (MODE .EQ. 'H') THEN DIVAS = 180.0 ELSE DIVAS = 360.0 ENDIF C MAKE EXPBUF AT LEAST NSAM*NROW FOR USE BY APSHIFT MWANTX = NSAM * NROW ALLOCATE(EXPBUF(MWANTX),STAT=IRTFLG) IF (IRTFLG .NE. 0) THEN CALL ERRT(46,'DSGR_P, EXPBUF',MWANTX) GOTO 9999 ENDIF C LOAD REF. PROJ. ANGLES (ANGREF) FROM DOC. FILE (REFANGDOC) OR C REF. IMAGE FILE (REFPAT) HEAD CALL AP_GETANGA(IREFLIST,NUMREF,0,REFANGDOC,REFPAT, & INPIC,LUNANG,3,ANGREF,NGOTREF,IRTFLG) IF (IRTFLG .NE. 0) GOTO 9999 C CONVERT REF. ANGLES TO UNITARY DIRECTIONAL VECTORS (REFDIR). CALL AP_GETSATA(ANGREF,REFDIR,3,NUMREF,IRTFLG) IF (IRTFLG .NE. 0) GOTO 9999 GOTREFANG = .TRUE. C READ REFERENCE IMAGES INTO REFERENCE RINGS (CIRCREF) ARRAY OR C CREATE REFERENCE RINGS FILE FOR LATER READING 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 IF (NSAID .LE. 0) THEN IF (CIRCREF_IN_CORE) THEN WRITE(NOUT,91)NUMTH 91 FORMAT(' Ref. rings in core, Threads: ',I4) ELSE WRITE(NOUT,92)NUMTH 92 FORMAT(' Ref. rings not in core, Threads: ',I4) ENDIF NSAID = NSAID + 1 ENDIF C LOAD EXP. PROJ. ANGLES & ALIGNMENT PARAMETERS (ANGEXP) C FROM DOC. FILE (EXPANGDOC) 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,LUNANG,8,ANGEXP,NGOTPAR,IRTFLG) IF (IRTFLG .NE. 0) GOTO 9999 ELSE ANGEXP = 0.0 ENDIF C INITIALIZE CCROT STATISTICS CALL AP_STAT_ADD(-1,CCROT,ANGDIF,ANGDIFTHR,CCROTLAS, & CCROTAVG,IBIGANGDIF,ANGDIFAVG,IMPROVCCROT, & CCROTIMPROV,IWORSECCROT,CCROTWORSE) C CALCULATE DIMENSIONS FOR APRINGS CNS2 = NSAM / 2 + 1 CNR2 = NROW / 2 + 1 DO IEXP=1,NUMEXP C LOOP OVER ALL EXPERIMENTAL (SAMPLE) IMAGES ------------------ C CONVERT EXP. ANGLE TO UNITARY DIRECTIONAL VECTORS (EXPDIR). CALL AP_GETSATA(ANGEXP(1,IEXP),EXPDIR,8,1,IRTFLG) C LOAD CURRENT (SINGLE) EXPERIMENTAL IMAGE INTO ARRAY EXPBUF CALL AP_GETDATS(IEXPLIST,NUMEXP,NSAM,NROW, & 1,EXPPAT,INPIC, IEXP,IEXP, & EXPBUF,AVI,SIGI,IRTFLG) IF (IRTFLG .NE. 0) GOTO 9999 C EXTRACT EXP. IMAGE POLAR COORD. RINGS, NORMALIZE & FFT THEM CALL APRINGS_ONE_NEW(NSAM,NROW, CNS2,CNR2, EXPBUF,.FALSE., & MODE,NUMR,NRING,LCIRC, 0.0,FFTW_PLANS, & CIRCEXP,IRTFLG) IF (IRTFLG .NE. 0) GOTO 9999 C DETERMINE WHICH REF IMAGES ARE TO BE COMPARED NUMREFLCG = 0 IEND = NUMREF CCROTD = -1.0D20 LIMITRANGE = RANGE .LT. 1.0 IF (LIMITRANGE) THEN C DETERMINE WHICH REF IMAGES ARE TO BE COMPARED DO IMI=1,NUMREF C LOOP OVER ALL REF. IMAGES C ABS - DIRECTIONS AT 180 DEGREES ARE DIFFERENT C (- DO NOT CHECK MIRRORED) C DT NEAR 1.0 = NOT-MIRRORED, DT NEAR -1.0 = MIRRORED DT = (EXPDIR(1) * REFDIR(1,IMI) + & EXPDIR(2) * REFDIR(2,IMI) + & EXPDIR(3) * REFDIR(3,IMI)) DTABS = ABS(DT) IF (DTABS .GE. RANGE) THEN C NON-MIRRORED OR MIRRORED REF. IS WITHIN RANGE NUMREFLCG = NUMREFLCG + 1 LCG(NUMREFLCG) = IMI IF (DT .LT. 0) LCG(NUMREFLCG) = -IMI ENDIF ENDDO IF (NUMREFLCG .LE. 0) THEN C NO REF. IMAGE WITHIN COMPARISON ANGLE IREF = 0 GOTO 1000 ENDIF IEND = NUMREFLCG ENDIF C REF. IMAGES FOUND WITHIN COMPARISION RANGE IF (CIRCREF_IN_CORE) THEN C USE CIRCREF FOR REFERENCE RINGS c$omp parallel do private(imil,imi,mirrored) DO IMIL=1,IEND IMI = IMIL IF (LIMITRANGE) IMI = ABS(LCG(IMIL)) IF ((CKMIRROR .AND. LIMITRANGE) .OR. & (.NOT. CKMIRROR)) THEN IF (.NOT. CKMIRROR) MIRRORED = .FALSE. IF (LIMITRANGE) MIRRORED = (LCG(IMIL) .LT. 0) C CHECK EITHER MIRRORED OR NON-MIRRORED POSITION CALL CROSRNG_EP_NEW(CIRCREF(1,IMI),CIRCEXP, & LCIRC,NRING, MAXRIN,NUMR, & TOTMIN(IMI),ROTANGT(IMI), & TT,MIRRORED,FFTW_PLANS(1)) ELSE C CHECK BOTH NON-MIRRORED & MIRRORED POSITIONS CALL CROSRNG_MSP_NEW(CIRCREF(1,IMI),CIRCEXP, & LCIRC,NRING, MAXRIN,NUMR, & TOTMIN(IMI),ROTANGT(IMI), & TOTMIR(IMI),ROTANGM(IMI), & TT,FFTW_PLANS(1)) ENDIF ENDDO c$omp end parallel do ELSE C USE REFERENCE RINGS FILE (MIGHT BE AN INCORE FILE) c$omp parallel do private(imil,imi,mirrored,ithread) c$omp& schedule(static,1) DO IMIL=1,IEND !LOOP OVER ALL REFERENCE IMAGES IMI = IMIL IF (LIMITRANGE) IMI = ABS(LCG(IMIL)) C FIND THREAD NUMBER ITHREAD = MOD((IMIL-1),NUMTH) + 1 C FILL CIRCREF FROM REFERENCE RINGS FILE c$omp critical CALL REDLIN(LUNRING,CIRCREF(1,ITHREAD),LCIRC,IMI) c$omp end critical IF ((CKMIRROR .AND. LIMITRANGE) .OR. & (.NOT. CKMIRROR)) THEN IF (.NOT. CKMIRROR) MIRRORED = .FALSE. IF (LIMITRANGE) MIRRORED = (LCG(IMIL) .LT. 0) C CHECK EITHER MIRRORED OR NON-MIRRORED POSITION CALL CROSRNG_EP_NEW(CIRCREF(1,ITHREAD),CIRCEXP, & LCIRC,NRING, MAXRIN,NUMR, & TOTMIN(IMI),ROTANGT(IMI), & TT,MIRRORED,FFTW_PLANS(1)) ELSE C CHECK BOTH NON-MIRRORED & MIRRORED POSITIONS CALL CROSRNG_MSP_NEW(CIRCREF(1,ITHREAD),CIRCEXP, & LCIRC,NRING, MAXRIN,NUMR, & TOTMIN(IMI),ROTANGT(IMI), & TOTMIR(IMI),ROTANGM(IMI), & TT,FFTW_PLANS(1)) ENDIF ENDDO comp end parallel do ENDIF C LOOP OVER ALL RELEVANT REF. IMAGES ------------------------- IREF = 0 CCROTD = -1.0D20 DO IMIL=1,IEND IMI = IMIL IF (NUMREFLCG .GT. 0) IMI = ABS(LCG(IMIL)) IF (TOTMIN(IMI) .GE. CCROTD) THEN C GOOD MATCH WITH TOTMIN (MIRRORED OR NOT) POSITION CCROTD = TOTMIN(IMI) RANGNEW = ROTANGT(IMI) MIRRORNEW = (LIMITRANGE .AND. (LCG(IMIL) .LT. 0)) IREF = IMI ENDIF IF (CKMIRROR .AND. .NOT. LIMITRANGE) THEN C HAVE TO COMPARE WITH MIRRORED POSITION IF (TOTMIR(IMI) .GE. CCROTD) THEN C GOOD MATCH, MIRRORED POSITION IS BETTER CCROTD = TOTMIR(IMI) RANGNEW = ROTANGM(IMI) MIRRORNEW = .TRUE. IREF = IMI ENDIF ENDIF ! END OF: IF (CKMIRROR) ENDDO ! END OF: DO IMIL=1,IEND ------------------------- 1000 RANGNEW = (RANGNEW-1) / MAXRIN * DIVAS CCROT = CCROTD IMGEXP = IEXPLIST(IEXP) PEAKV = 0.0 XSHNEW = 0.0 YSHNEW = 0.0 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 IF (IMGREF .GT. 0 .AND. ISHRANGE .GT. 0) THEN C DETERMINE SHIFT PARAMETERS, NO NEED TO RELOAD EXP. IMAGE NSAMP = 2*NSAM+2 NROWP = 2*NROW CALL APSHIFT(INPIC, REFPAT,IMGREF, & NSAM,NROW, NSAMP,NROWP, & EXPBUF,AVI,SIGI, ISHRANGE, & RANGNEW,XSHNEW,YSHNEW,MIRRORNEW,PEAKV,IRTFLG) IF (IRTFLG .NE. 0) RETURN ENDIF C AP_END WRITES ALIGNMENT PARAMETERS TO DOC FILE NPROJ = NUMREF IF (LIMITRANGE) NPROJ = NUMREFLCG CALL AP_END(IEXP,IMGEXP,IMGREF, & ANGREF(1,IREFT),REFDIR(1,IREFT), & ANGEXP(1,IEXP),EXPDIR,ISHRANGE, & GOTREFANG, NGOTPAR, CCROT,PEAKV, & RANGNEW,XSHNEW,YSHNEW,MIRRORNEW,REFPAT, & NPROJ,CTYPE, LUNDOC,PARLIST) CALL AP_STAT_ADD(NGOTPAR,CCROT,PARLIST(10), & ANGDIFTHR,ANGEXP(8,IEXP), & CCROTAVG,IBIGANGDIF,ANGDIFAVG,IMPROVCCROT, & CCROTIMPROV,IWORSECCROT,CCROTWORSE) ENDDO ! END OF: DO IEXP=1,NUMEXP LOOP OVER ALL EXP. IMAGES ---- IF (NUMEXP .GT. 1) THEN C SAVE CCROT & ANGULAR DISPLACEMENT STATISTICS CALL AP_STAT(NUMEXP,ANGDIFTHR,IBIGANGDIF, & ANGDIFAVG, CCROTAVG, & IMPROVCCROT,CCROTIMPROV, & IWORSECCROT,CCROTWORSE,LUNDOC) ENDIF 9999 CLOSE(LUNANG) IF (.NOT. CIRCREF_IN_CORE) CLOSE(LUNRING) IF (ALLOCATED(EXPBUF)) DEALLOCATE(EXPBUF) END #else C------------------------ MPI SPECIFIC SUBROUTINE -------------------- SUBROUTINE DSGR_P(IREFLIST,NUMREF,IEXPLIST,NUMEXP, & NSAM,NROW,LENTT,RANGE,ANGDIFTHR, & NRING,LCIRC,NUMR,CIRCREF,CIRCREF_IN_CORE, & 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 CHARACTER (LEN=*) :: REFANGDOC CHARACTER (LEN=*) :: EXPANGDOC CHARACTER (LEN=*) :: SCRFILE CHARACTER (LEN=*) :: REFPAT,EXPPAT CHARACTER (LEN=*) :: CTYPE CHARACTER (LEN=74) :: COMMENT DOUBLE PRECISION :: CCROTD CHARACTER(LEN=1) :: NULL LOGICAL :: CIRCREF_IN_CORE LOGICAL :: CKMIRROR LOGICAL :: MIRRORNEW LOGICAL :: GOTREFANG LOGICAL :: ONLYMIRROR LOGICAL :: LIMITRANGE LOGICAL :: MIRRORED INTEGER *8 :: FFTW_PLANS(*) C ALLOCATABLE ARRAYS REAL, ALLOCATABLE, DIMENSION(:) :: EXPBUF C AUTOMATIC ARRAYS DOUBLE PRECISION, DIMENSION(NUMREF) :: TOTMIN,TOTMIR DOUBLE PRECISION, DIMENSION(LENTT) :: TT REAL, DIMENSION(NUMREF) :: ROTANGT,ROTANGM REAL, DIMENSION(LCIRC) :: CIRCEXP REAL, DIMENSION(6) :: DLIST REAL, DIMENSION(3) :: EXPDIR REAL, DIMENSION(8,NUMEXP) :: ANGEXP REAL, DIMENSION(3,NUMREF) :: ANGREF,REFDIR INTEGER, DIMENSION(NUMREF) :: LCG INTEGER :: NSAID = 0 INTEGER, PARAMETER :: NLISTMAX = 15 REAL, DIMENSION(NLISTMAX) :: PARLIST PARAMETER (QUADPI = 3.1415926535897932384626) PARAMETER (DGR_TO_RAD = (QUADPI/180)) DATA INPIC,LUNANG,LUNRING/77,78,50/ INCLUDE 'mpif.h' INTEGER, ALLOCATABLE, DIMENSION(:) :: PSIZE INTEGER, ALLOCATABLE, DIMENSION(:) :: NBASE REAL, ALLOCATABLE, DIMENSION(:,:) :: PARTAB,PARTABLOC #ifdef MPI_DEBUG DOUBLE PRECISION TCOM0, TCOM1 #endif LOGICAL :: ONLYONE_RED,ONLYONE_WRT COMMON /COMM_MPI/ONLYONE_RED,ONLYONE_WRT ICOMM = MPI_COMM_WORLD CALL MPI_COMM_RANK(ICOMM, MYPID, MPIERR) CALL MPI_COMM_SIZE(ICOMM, NPROCS, MPIERR) NULL = CHAR(0) NR = NUMR(1,NRING) ! OUTER RING NUMBER MAXRIN = NUMR(3,NRING) ! ACTUAL LENGTH OF LONGEST RING #ifdef SP_LIBFFTW3 MAXRIN = NUMR(3,NRING) - 2 #endif C THIS ALTERS RANGE! RANGE = COS(RANGE*DGR_TO_RAD) C FIND NUMBER OF OMP THREADS CALL GETTHREADS(NUMTH) IF (MODE .EQ. 'H') THEN DIVAS = 180.0 ELSE DIVAS = 360.0 ENDIF C LOAD REF. PROJ. ANGLES (ANGREF) FROM DOC. FILE (REFANGDOC) OR C REF. IMAGE FILE (REFPAT) HEAD CALL AP_GETANGA(IREFLIST,NUMREF,0,REFANGDOC,REFPAT, & INPIC,LUNANG,3,ANGREF,NGOTREF,IRTFLG) IF (IRTFLG .NE. 0) GOTO 9999 C CONVERT REF. ANGLES TO UNITARY DIRECTIONAL VECTORS (REFDIR). CALL AP_GETSATA(ANGREF,REFDIR,3,NUMREF,IRTFLG) IF (IRTFLG .NE. 0) GOTO 9999 GOTREFANG = .TRUE. C READ REFERENCE IMAGES INTO REFERENCE RINGS (CIRCREF) ARRAY OR C CREATE REFERENCE RINGS FILE FOR LATER READING 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 IF (NSAID .LE. 0) THEN IF (CIRCREF_IN_CORE) THEN IF (MYPID .LE. 0) WRITE(NOUT,91)NUMTH 91 FORMAT(' Ref. rings in core, Threads: ',I4) ELSE IF (MYPID .LE. 0) WRITE(NOUT,92)NUMTH 92 FORMAT(' Ref. rings not in core, Threads: ',I4) ENDIF NSAID = NSAID + 1 ENDIF C LOAD EXP. PROJ. ANGLES & ALIGNMENT PARAMETERS (ANGEXP) C FROM DOC. FILE (EXPANGDOC) 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,LUNANG,8,ANGEXP,NGOTPAR,IRTFLG) IF (IRTFLG .NE. 0) GOTO 9999 ELSE ANGEXP = 0.0 ENDIF C INITIALIZE CCROT STATISTICS CALL AP_STAT_ADD(-1,CCROT,ANGDIF,ANGDIFTHR,CCROTLAS, & CCROTAVG,IBIGANGDIF,ANGDIFAVG,IMPROVCCROT, & CCROTIMPROV,IWORSECCROT,CCROTWORSE) C DISTRIBUTE EXPERIMENTAL IMAGES ALLOCATE(PSIZE(NPROCS), NBASE(NPROCS), STAT=IRTFLG) IF (IRTFLG .NE. 0) THEN CALL ERRT(46,'DSGR_P: NBASE',2*NPROCS) GOTO 9999 ENDIF C FILL PSIZE WITH PARTITION LIMITS, PSIZE(1) = NUMEXP / NPROCS CALL SETPART(NUMEXP, PSIZE, NBASE) #ifdef MPI_DEBUG WRITE(6,111) NBASE(MYPID+1), MYPID 111 FORMAT(' DSGR_P: NBASE(MYPID+1): ', I5, ' MYPID: ', I5) CALL FLUSHFILE(6) CALL MPI_BARRIER(ICOMM,MPIERR) #endif NLOC = PSIZE(MYPID+1) ALLOCATE(EXPBUF(NSAM*NROW), & PARTAB(15,NUMEXP), & PARTABLOC(15,NLOC), & STAT=IRTFLG) IF (IRTFLG .NE. 0) THEN MWANT = NSAM*NROW + 15*(NUMEXP+NLOC) CALL ERRT(46,'DSGR_P: PARTAB...',MWANT) GOTO 9999 ENDIF C ZERO THE ARRAYS PARTAB = 0.0 PARTABLOC = 0.0 C CALCULATE DIMENSIONS FOR APRINGS CNS2 = NSAM / 2 + 1 CNR2 = NROW / 2 + 1 C FOLLOWING LOOP PERFORMED BY ALL PROCESSORS ON LOCAL DATA c write(0,*)' dsgr_p; barrier 0,nloc,nprocs: ',nloc,nprocs,mypid CALL MPI_BARRIER(ICOMM,MPIERR) ONLYONE_RED = .FALSE. ! DO NOT BROADCAST LOADED DATA ONLYONE_WRT = .FALSE. ! NEEDED FOR NORM3 IN AP_GETDATS DO JLOC = 1, NLOC ! --------------------------------------- INUM = NBASE(MYPID+1) + JLOC IMGEXP = IEXPLIST(INUM) c write(0,*)' dsgr_p; jloc,inum: ',jloc,inum,imgexp,mypid C LOAD CURRENT EXPERIMENTAL IMAGE INTO ARRAY EXPBUF CALL AP_GETDATS(IEXPLIST,NUMEXP,NSAM,NROW, & 1,EXPPAT,INPIC, INUM,INUM, & EXPBUF,AVI,SIGI, IRTFLG) IF (IRTFLG .NE. 0) GOTO 9999 C CONVERT EXP. ANGLE TO UNITARY DIRECTIONAL VECTORS (EXPDIR). CALL AP_GETSATA(ANGEXP(1,INUM),EXPDIR,8,1,IRTFLG) C EXTRACT EXP. IMAGE POLAR COORD. RINGS, NORMALIZE & FFT THEM CALL APRINGS_ONE_NEW(NSAM,NROW, CNS2,CNR2, EXPBUF, & .FALSE.,MODE,NUMR,NRING,LCIRC, 0.0,FFTW_PLANS, & CIRCEXP,IRTFLG) IF (IRTFLG .NE. 0) GOTO 9999 C DETERMINE WHICH REF IMAGES ARE TO BE COMPARED NUMREFLCG = 0 IEND = NUMREF CCROTD = -1.0D20 LIMITRANGE = RANGE .LT. 1.0 IF (LIMITRANGE) THEN C DETERMINE WHICH REF IMAGES ARE TO BE COMPARED DO IMI=1,NUMREF C LOOP OVER ALL REF. IMAGES C ABS - DIRECTIONS AT 180 DEGREES ARE DIFFERENT C (- DO NOT CHECK MIRRORED) C DT NEAR 1.0 = NOT-MIRRORED, DT NEAR -1.0 = MIRRORED DT = (EXPDIR(1) * REFDIR(1,IMI) + & EXPDIR(2) * REFDIR(2,IMI) + & EXPDIR(3) * REFDIR(3,IMI)) DTABS = ABS(DT) IF (DTABS .GE. RANGE) THEN C NON-MIRRORED OR MIRRORED REF. IS WITHIN RANGE NUMREFLCG = NUMREFLCG + 1 LCG(NUMREFLCG) = IMI IF (DT .LT. 0) LCG(NUMREFLCG) = -IMI ENDIF ENDDO IF (NUMREFLCG .LE. 0) THEN C NO REF. IMAGE WITHIN COMPARISON ANGLE IREF = 0 GOTO 1000 ENDIF IEND = NUMREFLCG ENDIF C REF. IMAGES FOUND WITHIN COMPARISION RANGE IF (CIRCREF_IN_CORE) THEN C USE CIRCREF FOR REFERENCE RINGS DO IMIL=1,IEND IMI = IMIL IF (LIMITRANGE) IMI = ABS(LCG(IMIL)) IF ((CKMIRROR .AND. LIMITRANGE) .OR. & (.NOT. CKMIRROR)) THEN IF (.NOT. CKMIRROR) MIRRORED = .FALSE. IF (LIMITRANGE) MIRRORED = (LCG(IMIL) .LT. 0) C CHECK EITHER MIRRORED OR NON-MIRRORED POSITION CALL CROSRNG_EP_NEW(CIRCREF(1,IMI),CIRCEXP, & LCIRC,NRING, MAXRIN,NUMR, & TOTMIN(IMI),ROTANGT(IMI), & TT,MIRRORED,FFTW_PLANS(1)) ELSE C CHECK BOTH NON-MIRRORED & MIRRORED POSITIONS CALL CROSRNG_MSP_NEW(CIRCREF(1,IMI),CIRCEXP, & LCIRC,NRING, MAXRIN,NUMR, & TOTMIN(IMI),ROTANGT(IMI), & TOTMIR(IMI),ROTANGM(IMI), TT,FFTW_PLANS(1)) ENDIF ENDDO ! END OF: DO IMIL=1,IEND ELSE C USE REFERENCE RINGS FILE (MIGHT BE AN INCORE FILE) DO IMIL=1,IEND !LOOP OVER ALL REFERENCE IMAGES IMI = IMIL IF (LIMITRANGE) IMI = ABS(LCG(IMIL)) C FIND THREAD NUMBER ITHREAD = MOD((IMIL-1),NUMTH) + 1 C FILL CIRCREF FROM REFERENCE RINGS FILE CALL REDLIN(LUNRING,CIRCREF(1,ITHREAD),LCIRC,IMI) IF ((CKMIRROR .AND. LIMITRANGE) .OR. & (.NOT. CKMIRROR)) THEN IF (.NOT. CKMIRROR) MIRRORED = .FALSE. IF (LIMITRANGE) MIRRORED = (LCG(IMIL) .LT. 0) C CHECK EITHER MIRRORED OR NON-MIRRORED POSITION CALL CROSRNG_EP_NEW(CIRCREF(1,ITHREAD),CIRCEXP, & LCIRC,NRING, MAXRIN,NUMR, & TOTMIN(IMI),ROTANGT(IMI), & TT,MIRRORED,FFTW_PLANS(1)) ELSE C CHECK BOTH NON-MIRRORED & MIRRORED POSITIONS CALL CROSRNG_MSP_NEW(CIRCREF(1,ITHREAD),CIRCEXP, & LCIRC,NRING, MAXRIN,NUMR, & TOTMIN(IMI),ROTANGT(IMI), & TOTMIR(IMI),ROTANGM(IMI), & TT,FFTW_PLANS(1)) ENDIF ENDDO ! END OF: DO IMIL=1,IEND ENDIF C LOOP OVER ALL RELEVANT REF. IMAGES IREF = 0 CCROTD = -1.0D20 DO IMIL=1,IEND IMI = IMIL IF (NUMREFLCG .GT. 0) IMI = ABS(LCG(IMIL)) IF (TOTMIN(IMI) .GE. CCROTD) THEN C GOOD MATCH WITH TOTMIN (MIRRORED OR NOT) POSITION CCROTD = TOTMIN(IMI) RANGNEW = ROTANGT(IMI) MIRRORNEW = (LIMITRANGE .AND. (LCG(IMIL) .LT. 0)) IREF = IMI ENDIF IF (CKMIRROR .AND. .NOT. LIMITRANGE) THEN C HAVE TO COMPARE WITH MIRRORED POSITION IF (TOTMIR(IMI) .GE. CCROTD) THEN C GOOD MATCH, MIRRORED POSITION IS BETTER CCROTD = TOTMIR(IMI) RANGNEW = ROTANGM(IMI) MIRRORNEW = .TRUE. IREF = IMI ENDIF ENDIF ! END OF: IF (CKMIRROR) ENDDO ! END OF: DO IMIL=1,IEND 1000 RANGNEW = (RANGNEW-1) / MAXRIN * DIVAS CCROT = CCROTD PEAKV = 0.0 XSHNEW = 0.0 YSHNEW = 0.0 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 c write(0,*)' dsgr_p; above apshiftn,imgexp: ',imgexp,mypid 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,AVI,SIGI, ISHRANGE, & RANGNEW,XSHNEW,YSHNEW,MIRRORNEW,PEAKV,IRTFLG) IF (IRTFLG .NE. 0) RETURN ENDIF NPROJ = NUMREF IF (LIMITRANGE) NPROJ = NUMREFLCG C AP_END PUTS ALIGNMENT PARAMETERS FOR THIS IMAGE IN PARTABLOC CALL AP_END(IMGEXP,IMGEXP,IMGREF, & ANGREF(1,IREFT),REFDIR(1,IREFT), & ANGEXP(1,INUM),EXPDIR,ISHRANGE, & GOTREFANG, NGOTPAR, CCROT,PEAKV, & RANGNEW,XSHNEW,YSHNEW,MIRRORNEW,REFPAT, & NPROJ, CTYPE, LUNDOC,PARTABLOC(1,JLOC)) CALL AP_STAT_ADD(NGOTPAR,CCROT,PARTABLOC(10,JLOC), & ANGDIFTHR,ANGEXP(8,INUM), & CCROTAVG,IBIGANGDIF,ANGDIFAVG,IMPROVCCROT, & CCROTIMPROV,IWORSECCROT,CCROTWORSE) ENDDO ! END OF: DO JLOC = 1, NLOC ---------------------- c write(0,*) ' dsgr_p, barrier 2, mypid: ' ,mypid CALL MPI_BARRIER(ICOMM,MPIERR) ONLYONE_RED = .TRUE. ! SET BCAST FLAG BACK TO DEFAULT ONLYONE_WRT = .TRUE. ! SET BCAST FLAG BACK TO DEFAULT C COLLECT ALL ALIGNMENT PARAMETERS INTO PARTAB DO IPROC = 1, NPROCS PSIZE(IPROC) = 15 * PSIZE(IPROC) NBASE(IPROC) = 15 * NBASE(IPROC) ENDDO CALL MPI_ALLGATHERV(PARTABLOC, PSIZE(MYPID+1), & MPI_REAL, PARTAB, PSIZE, NBASE, & MPI_REAL, ICOMM, MPIERR) C WRITE IS SYNCHRONIZED WITHIN LUNDOCWRTDAT IF (LUNDOC .GT. 0) THEN C SAVE IN ALIGNMENT DOC FILE C <,<,<, MIR-REF#,IMG#,INPLANE<, SX,SY,NPROJ, C