C++********************************************************************* C C MRQLI_PS.F USED CMLIMIT AUG 00 ARDEAN LEITH C ADDED REF_CIRC FILE APR 01 ARDEAN LEITH C NORMASS -> NORMAS OCT 01 ARDEAN LEITH C PROMPTS JAN 02 ARDEAN LEITH C OPFILEC FEB 03 ARDEAN LEITH C APMASTER REWRITE AUG 03 ARDEAN LEITH C USED AP_OUT FEB 04 ARDEAN LEITH C NORMASS USED FOR ALTIX JUN 04 ARDEAN LEITH C INTERPOLATION ON EDGE BUG AUG 04 ARDEAN LEITH C CROSRNG_E SPEEDS UP AUG 04 ARDEAN LEITH C AP_END CALL HAS PARLIST OCT 04 ARDEAN LEITH C PEAKV = 1 JAN 05 ARDEAN LEITH C RANGNEWT OPT 64 BUG FIXED MAY 05 ARDEAN LEITH C DISCARD MIRROR ... JUN 06 ARDEAN LEITH C AP_STAT CALL JAN 07 ARDEAN LEITH C REF-RINGS FILE ADDED MAR 08 ARDEAN LEITH C CROSRNG_E REWRITE MAR 08 ARDEAN LEITH C IREFLIST VAR. NAME MAR 08 ARDEAN LEITH C APRINGS_ONE_NEW MAR 08 ARDEAN LEITH C AP_END CALL CHANGED NOV 08 ARDEAN LEITH C AP_STAT_ADD NOV 08 ARDEAN LEITH C NPROJ BUG AUG 09 ARDEAN LEITH C ********************************************************************** C=* This file is part of: * C=* SPIDER - Modular Image Processing System. Author: J. FRANK * C=* Copyright 1985-2009 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=* * 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 MRQLI_PS C C PURPOSE: FIND ROTATIONAL AND SHIFT PARAMETERS TO ALIGN A SERIES OF C REFERENCE IMAGES WITH SAMPLE IMAGES C USED IF: CIRCREF_IN_CORE .AND. NUMEXP .GE. NUMTH 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 LSAM,LROW 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 NOTE: IF USING MRQLI_PS, MOST MEMORY DEMAND APPEARS TO BE DEPENDENT C ON LCIRC & NUMREF. LCIRC IS THE TOTAL LENGTH OF THE ARRAY C THAT HOLDS THE CIRCULAR RINGS, SO IT IS DEPENDENT ON C NUMBER OF RINGS AND THEIR RADIUS. NUMREF IS NUMBER OF REFERENCE C IMAGES. BIGGEST ARRAY ALLOCATED IS: CIRCREF(LCIRC,NUMREF) C ANOTHER SMALL ALLOCATED ARRAY IS: A(NSAM,NROW,NUMTH) C FOR 83 IMAGES of 125x125 WITH RINGS AT 5...47 SIZE=45MB C ARRAYS ONLY APPEAR TO TAKE: 3.6MB? C C23456789012345678901234567890123456789012345678901234567890123456789012 C--********************************************************************* C --------------------- NON-MPI CODE -------------------------------- #ifndef USE_MPI SUBROUTINE MRQLI_PS(IREFLIST,NUMREF,IEXPLIST,NUMEXP, & LSAM,LROW,NR,LENTT,ISHRANGE,ISTEP, & NRING,LCIRC,NUMR,CIRCREF,CIRCREF_IN_CORE, & MODE, REFANGDOC,EXPANGDOC,SCRFILE,FFTW_PLANS, & REFPAT,EXPPAT,RANGE,CKMIRROR,CTYPE,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 LOGICAL :: CIRCREF_IN_CORE LOGICAL :: CKMIRROR LOGICAL :: MIRRORNEW LOGICAL :: GOTREFANG,LIMITRANGE CHARACTER(LEN=1) :: MODE CHARACTER (LEN=*) :: REFANGDOC,EXPANGDOC CHARACTER (LEN=*) :: SCRFILE CHARACTER (LEN=*) :: REFPAT,EXPPAT CHARACTER (LEN=*) :: CTYPE INTEGER *8 :: FFTW_PLANS(*) C AUTOMATIC ARRAYS DOUBLE PRECISION, DIMENSION(LENTT) :: TT REAL, DIMENSION(3) :: ANGOUT C ALLOCATED ARRAYS REAL, ALLOCATABLE, DIMENSION(:,:,:) :: A INTEGER, ALLOCATABLE, DIMENSION(:) :: NPROJA REAL, ALLOCATABLE, DIMENSION(:,:) :: DLIST REAL, ALLOCATABLE,DIMENSION(:,:) :: REFDIR,EXPDIR REAL, ALLOCATABLE,DIMENSION(:,:) :: ANGREF,ANGEXP INTEGER, PARAMETER :: NLISTMAX = 15 REAL, DIMENSION(NLISTMAX) :: PARLIST PARAMETER (QUADPI = 3.14159265358979323846) PARAMETER (DGR_TO_RAD = (QUADPI/180)) DATA INPIC/77/,INANG/78/,LUNRING/50/ MYPID = -1 C SET TYPE OF OUTPUT DOC FILES WANTED NWANTOUT = 7 IF (CTYPE(1:2) .EQ. 'SH') NWANTOUT = 15 C INITIALIZE CCROT STATISTICS COUNTERS ANGDIFTHR = 0.0 CALL AP_STAT_ADD(-1,CCROT,ANGDIF,ANGDIFTHR,CCROTLAS, & CCROTAVG,IBIGANGDIF,ANGDIFAVG,IMPROVCCROT, & CCROTIMPROV,IWORSECCROT,CCROTWORSE) LIMITRANGE = (RANGE .GT. 0.0) RANGECOS = COS(RANGE * DGR_TO_RAD) !write(6,*) ' range,rangecos:',range,rangecos C FIND DIVAS, NUMTH, NSAM, & NROW CALL APMASTER_1(MODE,DIVAS,NR,NUMTH,LSAM,LROW,NSAM,NROW, & TT,LENTT) !write(6,*) ' aprings_new called ' C READ REFERENCE IMAGES INTO REFERENCE RINGS (CIRCREF) ARRAY CALL APRINGS_NEW(IREFLIST,NUMREF, LSAM,LROW, & NRING,LCIRC,NUMR,MODE,FFTW_PLANS, & REFPAT,INPIC,CIRCREF,CIRCREF_IN_CORE, & LUNRING,SCRFILE,IRTFLG) NSAM = LSAM NROW = LROW ALLOCATE(A(LSAM,LROW,NUMTH), & DLIST(5,NUMTH), & NPROJA(NUMTH), & STAT=IRTFLG) IF (IRTFLG .NE. 0) THEN MWANT = LSAM*LROW*NUMTH + 6*NUMTH CALL ERRT(46,'A & DLIST',MWANT) RETURN ENDIF DLIST = 0.0 ! ZEROS WHOLE DLIST GOTREFANG = .FALSE. IF (LIMITRANGE .OR. CTYPE(1:2) .EQ. 'SH') THEN C REFANGLES FILE FOR RESTRICTED ANGULAR SEARCH OR 'SH' ALLOCATE(REFDIR(3,NUMREF), & ANGREF(3,NUMREF), STAT=IRTFLG) IF (IRTFLG .NE. 0) THEN MWANT = 6*NUMREF CALL ERRT(46,'REFDIR, ANGREF',MWANT) RETURN ENDIF GOTREFANG = .TRUE. C READ REF. ANGLES INTO ANGREF CALL AP_GETANGA(IREFLIST,NUMREF,0,REFANGDOC,REFPAT, & INPIC,INANG,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) ENDIF NGOTPAR = 0 IF (CTYPE(1:2) .EQ. 'RQ' .OR. & (CTYPE(1:2) .EQ. 'SH' .AND. & EXPANGDOC .NE. CHAR(0))) THEN C READ EXP. ANGLES INTO ANGEXP ALLOCATE(ANGEXP(8,NUMEXP), EXPDIR(3,NUMEXP), STAT=IRTFLG) IF (IRTFLG .NE. 0) THEN MWANT = 11*NUMEXP CALL ERRT(46,'ANGEXP....',MWANT) RETURN ENDIF 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) ELSE C DUMMY ALLOCATE TO AVOID BUS ERROR ON SOME SYSTEMS ALLOCATE(ANGEXP(8,1), EXPDIR(3,1), STAT=IRTFLG) IF (IRTFLG .NE. 0) THEN MWANT = 11 CALL ERRT(46,'ANGEXP....',MWANT) RETURN ENDIF ENDIF C LOOP OVER ALL SETS OF EXPERIMENTAL (SAMPLE) IMAGES ----------- DO IEXPT=1,NUMEXP,NUMTH C LOAD EXP. IMAGE DATA FOR THIS SET OF IMAGES IEND = MIN(NUMEXP,IEXPT+NUMTH-1) CALL AP_GETDAT(IEXPLIST,NUMEXP,LSAM,LROW,LSAM,LROW, & NUMTH,EXPPAT,INPIC, IEXPT,IEND, & 1,LROW,1,LSAM, A, & IRTFLG) IF (IRTFLG .NE. 0) GOTO 9999 C NUMTH INPUT IMAGES READY TO BE ALIGNED c$omp parallel do private(iexp,it) DO IEXP=IEXPT,MIN(NUMEXP,IEXPT+NUMTH-1) IT = IEXP - IEXPT + 1 CALL APRQ2D(A(1,1,IT),CIRCREF,TT,NUMR, & NSAM,NROW,ISHRANGE,ISTEP, & LCIRC,NRING,NUMREF,MODE, & REFDIR,EXPDIR(1,IEXP),RANGECOS, & DLIST(1,IT),DLIST(2,IT), & DLIST(3,IT),DLIST(4,IT), & DLIST(5,IT),NPROJA(IT), & CKMIRROR,LIMITRANGE,FFTW_PLANS) ENDDO c$omp end parallel do C OUTPUT FROM DLIST POSITION C 1 - NUMBER OF MOST SIMILAR REFERENCE PROJECTION. C 2 - NOT-NORMALIZED CORRELATION COEFFICIENT. C 3 - PSI ANGLE. (IN=PLANE ROTATION) RANGNEW C 4 - SX SHIFT C 5 - SY SHIFT DO IEXP=IEXPT,MIN(NUMEXP,IEXPT+NUMTH-1) IT = IEXP - IEXPT + 1 IMGEXP = IEXPLIST(IEXP) C DLIST(1,IT) IS LIST NUMBER OF MOST SIMILAR REF. IMAGE C (<0 IF MIRRORED, 0 IF NO SIMILAR IMAGE ) IREF = INT(DLIST(1,IT)) IF (IREF .LT. 0) THEN C MIRRORED REFERENCE IMAGE IMGREF = IREFLIST(-IREF) C IREFT IS FOR REFDIR INDEX IREFT = -IREF MIRRORNEW = .TRUE. ELSEIF (IREF .EQ. 0) THEN C NO NEARBY REFERENCE IMAGE IMGREF = 0 C IREFT IS FOR REFDIR INDEX IREFT = 1 MIRRORNEW = .FALSE. ELSE IMGREF = IREFLIST(IREF) C IREFT IS FOR REFDIR INDEX IREFT = IREF MIRRORNEW = .FALSE. ENDIF CCROT = DLIST(2,IT) RANGNEW = DLIST(3,IT) XSHNEW = DLIST(4,IT) YSHNEW = DLIST(5,IT) NPROJ = NPROJA(IT) PEAKV = 1.0 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, 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,LUNDOC) ENDIF 9999 CONTINUE C DEALLOCATE ARRAYS IF (ALLOCATED(DLIST)) DEALLOCATE(DLIST) IF (ALLOCATED(NPROJA)) DEALLOCATE(NPROJA) IF (ALLOCATED(A)) DEALLOCATE(A) IF (ALLOCATED(REFDIR)) DEALLOCATE(REFDIR) IF (ALLOCATED(ANGEXP)) DEALLOCATE(ANGEXP) IF (ALLOCATED(ANGREF)) DEALLOCATE(ANGREF) END #else C------------------------ MPI SPECIFIC SUBROUTINE -------------------- SUBROUTINE MRQLI_PS(IREFLIST,NUMREF,IEXPLIST,NUMEXP, & LSAM,LROW,NR,LENTT,ISHRANGE,ISTEP, & NRING,LCIRC,NUMR,CIRCREF,CIRCREF_IN_CORE, & MODE, REFANGDOC,EXPANGDOC,SCRFILE,FFTW_PLANS, & REFPAT,EXPPAT,RANGE,CKMIRROR,CTYPE,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 LOGICAL :: CIRCREF_IN_CORE LOGICAL :: CKMIRROR,MIRRORNEW LOGICAL :: GOTREFANG,LIMITRANGE CHARACTER (LEN=1) :: MODE CHARACTER (LEN=*) :: REFANGDOC,EXPANGDOC CHARACTER (LEN=*) :: SCRFILE CHARACTER (LEN=*) :: REFPAT,EXPPAT CHARACTER (LEN=*) :: CTYPE INTEGER *8 :: FFTW_PLANS(*) C AUTOMATIC ARRAYS DOUBLE PRECISION, DIMENSION(LENTT) :: TT REAL, DIMENSION(3) :: ANGOUT C ALLOCATED ARRAYS REAL, ALLOCATABLE, DIMENSION(:,:,:) :: A REAL, ALLOCATABLE, DIMENSION(:,:) :: DLIST REAL, ALLOCATABLE,DIMENSION(:,:) :: REFDIR,EXPDIR REAL, ALLOCATABLE,DIMENSION(:,:) :: ANGREF,ANGEXP INTEGER, PARAMETER :: NLISTMAX = 15 REAL, DIMENSION(NLISTMAX) :: PARLIST PARAMETER (QUADPI = 3.14159265358979323846) PARAMETER (DGR_TO_RAD = (QUADPI/180)) DATA INPIC/77/,INANG/78/,NSCF/50/ include 'mpif.h' integer mypid, comm, mpierr, nprocs, nexploc, nrem integer istat(mpi_status_size) real , allocatable, dimension(:,:,:) :: aloc real , allocatable, dimension(:,:,:) :: abuf real , allocatable, dimension(:,:) :: DLISTLOC integer, allocatable, dimension(:) :: NPROJA real , allocatable, dimension(:,:) :: partab,partabloc integer, allocatable, dimension(:) :: nbase, psize integer iproc, isam, jrow, jloc, iloc, tag, master, & iti, imit, iglb, ibeg, iend #ifdef MPI_DEBUG double precision tcom1, tcom2 #endif c comm = mpi_comm_world call mpi_comm_rank(comm, mypid , mpierr) call mpi_comm_size(comm, nprocs, mpierr) master = 0 C SET TYPE OF OUTPUT DOC FILES WANTED NWANTOUT = 7 IF (CTYPE(1:2) .EQ. 'SH') NWANTOUT = 15 C INITIALIZE CCROT STATISTICS COUNTERS ANGDIFTHR = 0.0 C INITIALIZE CCROT STATISTICS CALL AP_STAT_ADD(-1,CCROT,ANGDIF,ANGDIFTHR,CCROTLAS, & CCROTAVG,IBIGANGDIF,ANGDIFAVG,IMPROVCCROT, & CCROTIMPROV,IWORSECCROT,CCROTWORSE) LIMITRANGE = (RANGE .GT. 0.0) MAXRIN = NUMR(3,NRING) RANGECOS = COS(RANGE*DGR_TO_RAD) C FIND DIVAS, NUMTH, NSAM, & NROW CALL APMASTER_1(MODE,DIVAS,NR,NUMTH,LSAM,LROW,NSAM,NROW, & TT,LENTT) C READ REFERENCE IMAGES INTO REFERENCE RINGS (CIRCREF) ARRAY CALL APRINGS_NEW(IREFLIST,NUMREF, LSAM,LROW, & NRING,LCIRC,NUMR,MODE,FFTW_PLANS, & REFPAT,INPIC,CIRCREF,CIRCREF_IN_CORE, & NSCF,SCRFILE,IRTFLG) NSAM = LSAM NROW = LROW GOTREFANG = .FALSE. IF (LIMITRANGE .OR. CTYPE(1:2) .EQ. 'SH') THEN C REFANGLES FILE FOR RESTRICTED ANGULAR SEARCH OR 'SH' ALLOCATE(REFDIR(3,NUMREF), & ANGREF(3,NUMREF), STAT=IRTFLG) IF (IRTFLG .NE. 0) THEN MWANT = 6*NUMREF CALL ERRT(46,'REFDIR, ANGREF',MWANT) RETURN ENDIF GOTREFANG = .TRUE. C READ REF. ANGLES INTO ANGREF CALL AP_GETANGA(IREFLIST,NUMREF,0,REFANGDOC,REFPAT, & INPIC,INANG,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) ENDIF NGOTPAR = 0 IF (CTYPE(1:2) .EQ. 'RQ' .OR. & (CTYPE(1:2) .EQ. 'SH' .AND. & EXPANGDOC .NE. CHAR(0))) THEN C READ EXP. ANGLES INTO ANGEXP ALLOCATE(ANGEXP(8,NUMEXP), EXPDIR(3,NUMEXP), STAT=IRTFLG) IF (IRTFLG .NE. 0) THEN MWANT = 11*NUMEXP CALL ERRT(46,'ANGEXP....',MWANT) RETURN ENDIF 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,7,NUMEXP,IRTFLG) ELSE C DUMMY ALLOCATE TO AVOID BUS ERROR ON SOME SYSTEMS ALLOCATE(ANGEXP(8,1), EXPDIR(3,1), STAT=IRTFLG) IF (IRTFLG .NE. 0) THEN MWANT = 10*1 CALL ERRT(46,'ANGEXP....',MWANT) RETURN ENDIF ENDIF C === PARTITION AND DISTRIBUTE EXPERIMENTAL IMAGES === ALLOCATE(PSIZE(NPROCS), & NBASE(NPROCS), STAT=IRTFLG) IF (IRTFLG .NE. 0) THEN WRITE(6,*) 'MRQLI_PS: FAILED TO ALLOCATE PSIZE & NBASE' RETURN ENDIF CALL SETPART(NUMEXP, PSIZE, NBASE) NEXPLOC = PSIZE(MYPID+1) ALLOCATE(ALOC(NSAM,NROW,NEXPLOC), & ABUF(NSAM,NROW,PSIZE(1)),STAT=IRTFLG) IF (IRTFLG .NE. 0) THEN WRITE(6,*) ' MRQLI_PS: FAILED TO ALLOCATE ABUF & ALOC' RETURN ENDIF ALOC = 0.0 #ifdef MPI_DEBUG WRITE(6,111) NBASE(MYPID+1), MYPID CALL FLUSHFILE(6) 111 FORMAT(' MRQLI_PS: NBASE = ', I5, ' MYPID = ', I5) TCOM0 = MPI_WTIME() #endif C PROCESS 0 READS IMAGES & DISTRIBUTES TO DIFFERENT PROCESSORS. === DO IPROC = 1, NPROCS NLOC = PSIZE(IPROC) C === CALCULATE THE GLOBAL INDEX === IBEG = NBASE(IPROC) + 1 IEND = NBASE(IPROC) + NLOC C ALTHOUGH THE FOLLOWING IS CALLED BY ALL PROCESSORS, C ONLY ONE PROCESSOR (MYPID=0) READS IMAGES INTO ABUF CALL AP_GETDAT1P(IEXPLIST,NUMEXP,LSAM,LROW,LSAM,LROW, & NUMTH,EXPPAT,INPIC, IBEG,IEND, & 1,LROW,1,LSAM, ABUF,IRTFLG) IF (IRTFLG .NE. 0) STOP IF (IPROC .GT. 1) THEN IF (MYPID .EQ. 0) THEN #ifdef MPI_DEBUG WRITE(6,222) IPROC-1 CALL FLUSHFILE(6) 222 FORMAT(' MRQLI_PS: SENDING TO PID = ', I3) #endif CALL MPI_SEND(ABUF , NSAM*NROW*NLOC, MPI_REAL, & IPROC-1, IPROC-1 , COMM , & MPIERR) ELSE IF (MYPID .EQ. IPROC-1) THEN C === SLAVES RECEIVE LOCAL PIECES === CALL MPI_RECV(ALOC , NSAM*NROW*NLOC, MPI_REAL, & 0 , MPI_ANY_TAG , COMM , & ISTAT, MPIERR) IF (MPIERR .NE. 0) THEN WRITE(6,*) ' RECV FAILED' STOP ENDIF #ifdef MPI_DEBUG WRITE(6,223) MYPID CALL FLUSHFILE(6) 223 FORMAT(' MRQLI_PS: RECEIVED BY MYPID = ', I3) #endif ENDIF ELSE IF (MYPID .EQ. 0) THEN C === SIMPLY COPY FROM ABUF TO ALOC === DO JLOC = 1, NLOC DO ISAM = 1, NSAM DO JROW = 1, NROW ALOC(ISAM,JROW,JLOC) = ABUF(ISAM,JROW,JLOC) ENDDO ENDDO ENDDO ENDIF ENDDO IF (ALLOCATED(ABUF)) DEALLOCATE(ABUF) ALLOCATE(PARTAB(15,NUMEXP), & PARTABLOC(15,NEXPLOC), & DLISTLOC(5,NEXPLOC), & NPROJA(NEXPLOC), & STAT=IRTFLG) IF (IRTFLG .NE. 0) THEN MWANT = 15*NUMEXP + 15*NEXPLOC + 6*NEXPLOC CALL ERRT(46,'MRQLI_PS; PARTAB..',MWANT) RETURN ENDIF #ifdef MPI_DEBUG TCOM1 = MPI_WTIME() IF (MYPID .EQ. 0) WRITE(6,440) TCOM1-TCOM0 440 FORMAT(' MRQLI_PS: DATA DIST TIME = ', 1PE11.3) WRITE(6,444) MYPID 444 FORMAT(' MRQLI_PS: CALLING APRQ2D.., MYPID = ', I3) #endif PARTAB = 0.0 PARTABLOC = 0.0 DO IT=1,NEXPLOC IGLB = NBASE(MYPID+1) + IT CALL APRQ2D(ALOC(1,1,IT),CIRCREF,TT,NUMR, & NSAM,NROW,ISHRANGE,ISTEP, & LCIRC,NRING,NUMREF,MODE, & REFDIR,EXPDIR(1,IGLB),RANGECOS, & DLISTLOC(1,IT),DLISTLOC(2,IT), & DLISTLOC(3,IT),DLISTLOC(4,IT), & DLISTLOC(5,IT),NPROJA(IT), & CKMIRROR,LIMITRANGE, FFTW_PLANS) C OUTPUT IN DLISTLOC C 1 - NUMBER OF THE MOST SIMILAR REFERENCE PROJECTION. C 2 - NOT-NORMALIZED CORRELATION COEFFICIENT. C 3 - PSI ANGLE. (IN=PLANE ROTATION) C 4 - SX SHIFT C 5 - SY SHIFT IMGEXP = IEXPLIST(IGLB) C DLISTLOC(1,IT) IS LIST NUMBER OF MOST SIMILAR REF. IMAGE C (<0 IF MIRRORED, 0 IF NONE ) IREF = INT(DLISTLOC(1,IT)) IF (IREF .LT. 0) THEN C MIRRORED REFERENCE IMAGE IMGREF = IREFLIST(-IREF) C IREFT IS FOR REFDIR INDEX IREFT = -IREF MIRRORNEW = .TRUE. ELSEIF (IREF .EQ. 0) THEN C NO NEARBY REFERENCE IMAGE IMGREF = 0 C IREFT IS FOR REFDIR INDEX IREFT = 1 MIRRORNEW = .FALSE. ELSE IMGREF = IREFLIST(IREF) C IREFT IS FOR REFDIR INDEX IREFT = IREF MIRRORNEW = .FALSE. ENDIF CCROT = DLISTLOC(2,IT) RANGNEW = DLISTLOC(3,IT) XSHNEW = DLISTLOC(4,IT) YSHNEW = DLISTLOC(5,IT) NPROJ = NPROJA(IT) PEAKV = 1.0 CALL AP_END(IGLB,IMGEXP,IMGREF, & ANGREF(1,IREFT),REFDIR(1,IREFT), & ANGEXP(1,IGLB), EXPDIR(1,IGLB),ISHRANGE, & GOTREFANG, NGOTPAR, CCROT,PEAKV, & RANGNEW,XSHNEW,YSHNEW,MIRRORNEW, REFPAT, & NPROJ, CTYPE, LUNDOC, PARTABLOC(1,IT)) ENDDO C 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, COMM , MPIERR) C WRITE HAS BEEN SYNCHRONIZED WITHIN LUNDOCWRTDAT IF (NWANTOUT .GT. 7) THEN IF (LUNDOC .GT. 0) THEN C SAVE IN ALIGNMENT DOC FILE C <,<,<, MIR-REF#,IMG#,INPLANE<, SX,SY,NPROJ,DIF,CCROT,INPLANE<,SX,SY DO IT = 1, NUMEXP CALL LUNDOCWRTDAT(LUNDOC,IT,PARTAB(1,IT), & NWANTOUT,IRTFLG) ENDDO ENDIF ELSE IF (LUNDOC .GT. 0) THEN C SAVE IN ALIGNMENT DOC FILE C MIR-REF#, CCROT, INPLANE<, SX,SY, IMG#, < DIFF DO IT = 1, NUMEXP CALL LUNDOCWRTDAT(LUNDOC,IT,PARTAB(1,IT), & NWANTOUT,IRTFLG) ENDDO ENDIF ENDIF 9999 CONTINUE C DEALLOCATE ARRAYS IF (ALLOCATED(ALOC)) DEALLOCATE(ALOC) IF (ALLOCATED(PSIZE)) DEALLOCATE(PSIZE) IF (ALLOCATED(NBASE)) DEALLOCATE(NBASE) IF (ALLOCATED(NPROJA)) DEALLOCATE(NPROJA) IF (ALLOCATED(DLISTLOC)) DEALLOCATE(DLISTLOC) IF (ALLOCATED(PARTAB)) DEALLOCATE(PARTAB) IF (ALLOCATED(PARTABLOC)) DEALLOCATE(PARTABLOC) IF (ALLOCATED(REFDIR)) DEALLOCATE(REFDIR) IF (ALLOCATED(ANGEXP)) DEALLOCATE(ANGEXP) IF (ALLOCATED(ANGREF)) DEALLOCATE(ANGREF) END C ------------------------- END OF MPI CODE ----------------------- #endif C+********************************************************************** C C APRQ2D.F C C PARAMETERS: C DIREF NUMBER OF MOST SIMILAR REF. PROJ. (OUTPUT) C (NEGATIVE IF MIRRORED) C CCROT CORR COEFF. (OUTPUT) C RANGNEW INPLANE ANGLE (OUTPUT) C XSHSUM SHIFT (OUTPUT) C YSHSUM SHIFT (OUTPUT) C NPROJ (OUTPUT) C C-********************************************************************** SUBROUTINE APRQ2D(EXPBUF,CIRCREF,TT,NUMR, & NSAM,NROW,ISHRANGE,ISTEP, & LCIRC,NRING,NUMREF,MODE, & REFDIR,EXPDIR,RANGECOS, & DIREF,CCROT,RANGNEW,XSHSUM,YSHSUM,NPROJ, & CKMIRROR,LIMITRANGE,FFTW_PLANS) C NOTE: RUNS WITHIN OMP PARALLEL SECTION OF CODE IF NOT UNDER MPI! REAL :: EXPBUF(NSAM,NROW) REAL :: CIRCREF(LCIRC,NUMREF) INTEGER :: NUMR(3,NRING) DOUBLE PRECISION :: FITP(-1:1,-1:1) DOUBLE PRECISION :: TT(*) CHARACTER (LEN=1) :: MODE REAL :: REFDIR(3,NUMREF), EXPDIR(3) INTEGER *8 :: FFTW_PLANS(*) C AUTOMATIC ARRAYS DOUBLE PRECISION :: FIT(-ISTEP:ISTEP,-ISTEP:ISTEP) REAL :: ROTMP(-ISTEP:ISTEP,-ISTEP:ISTEP) REAL :: CIRCEXP(LCIRC) INTEGER, ALLOCATABLE :: LCG(:) DOUBLE PRECISION :: CCROTD,PEAK,CCROTD_INTERP DOUBLE PRECISION :: CCOA,CCMA LOGICAL :: CKMIRROR,LIMITRANGE LOGICAL :: MIRRORED PARAMETER (QUADPI = 3.14159265358979323846) PARAMETER (DGR_TO_RAD = (QUADPI/180)) PEAK = 0.0 WR = 0.0 ! DUMMY VALUE FLAG FOR APRINGS CALL MAXRIN = NUMR(3,NRING) ! ACTUAL LENGTH OF LONGEST RING #ifdef SP_LIBFFTW3 MAXRIN = NUMR(3,NRING) - 2 #endif IEND = NUMREF IF (.NOT. LIMITRANGE) THEN C DUMMY ALLOCATE TO AVOID BUS ERROR ON SOME SYSTEMS ALLOCATE(LCG(1),STAT=IRTFLG) NPROJ = NUMREF ELSE C RESTRICTED RANGE SEARCH ALLOCATE(LCG(NUMREF),STAT=IRTFLG) IF (IRTFLG .NE. 0) THEN CALL ERRT(46,'LCG',NUMREF) RETURN ENDIF C MAKE LIST OF NEARBY REFERENCE IMAGES NPROJ = 0 C write(6,*) 'numref,range: ',numref,rangecos,'expdir: ',expdir DO IREF=1,NUMREF C DT NEAR 1.0 = NOT-MIRRORED, DT NEAR -1.0 = MIRRORED DT = (EXPDIR(1) * REFDIR(1,IREF) + & EXPDIR(2) * REFDIR(2,IREF) + & EXPDIR(3) * REFDIR(3,IREF)) DTABS = ABS(DT) IF (DTABS .GE. RANGECOS) THEN C MIRRORED OR NON-MIRRORED IS WITHIN RANGE IF (CKMIRROR .OR. DT .GT. 0) THEN C DO NOT DISCARD IF NOT MIRRORED OR WANT MIRRORED NPROJ = NPROJ + 1 LCG(NPROJ) = IREF IF (DT .LT. 0) LCG(NPROJ) = -IREF c if (mod(iref,50) .le. 0) then c write(110,*)'in iref: ',iref,' dt: ',dt,lcg(nproj) c endif ENDIF ENDIF ENDDO IF (NPROJ .LE. 0) THEN C THERE IS NO REFERENCE WITHIN SEARCH RANGE XSHSUM = 0 YSHSUM = 0 DIREF = 0 RANGNEW = 0 CCROT = -1.0 RETURN ENDIF IEND = NPROJ ENDIF ! END OF RESTRICTED RANGE SEARCH CCROTD = -HUGE(CCROTD) NEWMAX = .FALSE. c GO THROUGH CENTERS FOR SHIFT ALIGNMENT ------------------------ DO JT=-ISHRANGE,ISHRANGE,ISTEP CNR2 = NROW / 2 + 1 + JT DO IT=-ISHRANGE,ISHRANGE,ISTEP CNS2 = NSAM / 2 + 1 + IT C NORMALIZE IMAGE VALUES UNDER THE MASK OVER VARIANCE RANGE C INTERPOLATE TO POLAR COORD., CREATE FOURIER OF: A_CIRC C NO WEIGHTING OF RINGS !write(6,*) ' aprings one:,',jt,it,cns2,cnr2 CALL APRINGS_ONE_NEW(NSAM,NROW, CNS2,CNR2, EXPBUF,.FALSE., & MODE,NUMR,NRING,LCIRC,WR,FFTW_PLANS, & CIRCEXP,IRTFLG) C COMPARE EXP. IMAGE WITH ALL REFERENCE IMAGES ------------ DO IRR=1,IEND IR = IRR IF (LIMITRANGE) IR = ABS(LCG(IRR)) IF ((CKMIRROR .AND. LIMITRANGE) .OR. & (.NOT. CKMIRROR)) THEN IF (.NOT. CKMIRROR) MIRRORED = .FALSE. IF (LIMITRANGE) MIRRORED = (LCG(IRR) .LT. 0) CALL CROSRNG_EP_NEW(CIRCREF(1,IR),CIRCEXP, & LCIRC,NRING, MAXRIN,NUMR, & CCOA,RANGO, & TT,MIRRORED, FFTW_PLANS(1)) ELSE C CHECK BOTH NON-MIRRORED & MIRRORED POSITIONS CALL CROSRNG_MSP_NEW(CIRCREF(1,IR),CIRCEXP, & LCIRC,NRING,MAXRIN,NUMR, & CCOA,RANGO, CCMA,RANGM, & TT, FFTW_PLANS(1)) #ifdef NEVER if(it .eq. -10 .and. jt .eq. -1 .and. ir .eq. 47)then rt1 = ang_n(tot,mode,maxrin) write(6,821) it,jt, tota,rt1 821 format(' tota(',i3,',',i3,'): ',f14.4,' ',2f8.2) rt1 = ang_n(tmt,mode,maxrin) write(6,826) it,jt, tmta,rt1 826 format(' tmta(',i3,',',i3,'): ',f14.4,' ',2f8.2) endif #endif ENDIF IF (CCOA .GE. CCROTD) THEN C GOOD MATCH WITH CCOA (MIRRORED OR NOT) POSITION CCROTD = CCOA IBE = IR ISX = IT ISY = JT RANGNEW = ANG_N(RANGO,MODE,MAXRIN) IDIS = IR IF (LIMITRANGE .AND. (LCG(IRR) .LT. 0)) IDIS = -IR ENDIF IF (CKMIRROR .AND. .NOT. LIMITRANGE) THEN C HAVE TO COMPARE WITH MIRRORED POSITION IF (CCMA .GE. CCROTD) THEN C GOOD MATCH WITH MIRRORED POSITION CCROTD = CCMA IBE = IR ISX = IT ISY = JT RANGNEW = ANG_N(RANGM,MODE,MAXRIN) IDIS = -IR ENDIF ENDIF ENDDO ! END OF: DO IRR=1,IEND ---------------------- REFS. ENDDO ! END OF: DO IT=-ISHRANGE,ISHRANGE,ISTEP ENDDO ! END OF: DO JT=-ISHRANGE,ISHRANGE,ISTEP SX = ISX SY = ISY CCROT = CCROTD #ifdef NEVER write(6,921) idis,isx,isy,ccrotd,rangnew 921 format(' 1 ',i5,' (',i3,',',i3,'): ',f12.4,' ',2f8.2,f6.1) #endif C CHECK LOCATIONS WITHIN ISHRANGE AROUND MAX ------------------ DIREF = IDIS C WHEN (IDIS .LT. 0) CHECK MIRRORED ONLY MIRRORED = (IDIS .LT. 0) C DO NOT INTERPOLATE FOR POINT ON THE EDGE IF (IABS(ISX).NE.ISHRANGE .AND. IABS(ISY).NE.ISHRANGE) THEN C NOT ON BOUNDARY, HAVE TO FIND NEIGHBOURING VALUES FIT(0,0) = CCROTD ROTMP(0,0) = RANGNEW DO JT=-ISTEP,ISTEP CNR2 = NROW / 2 + 1 + JT + ISY DO IT=-ISTEP,ISTEP CNS2 = NSAM / 2 + 1 + IT + ISX IF (IT.NE.0 .OR. JT.NE.0) THEN C NORMALIZE IMAGE VALUES UNDER THE MASK OVER VARIANCE RANGE C INTERPOLATE INTO POLAR COORDINATES C CREATES FOURIER OF: CIRCEXP !!write(6,*) ' aprings max:,',jt,it,cns2,cnr2 CALL APRINGS_ONE_NEW(NSAM,NROW, CNS2,CNR2, & EXPBUF,.FALSE., & MODE,NUMR,NRING,LCIRC,WR,FFTW_PLANS, & CIRCEXP,IRTFLG) CALL CROSRNG_EP_NEW(CIRCREF(1,IBE),CIRCEXP, & LCIRC,NRING,MAXRIN,NUMR, & FIT(IT,JT),ROTMP(IT,JT), & TT,MIRRORED, FFTW_PLANS(1)) C RECORD BEST ANGLE IN ROTMP ROTMP(IT,JT) = ANG_N(ROTMP(IT,JT),MODE,MAXRIN) ENDIF ! END OF: IF (IT.NE.0 .OR. JT.NE.0) ENDDO ! END OF: DO IT=-ISTEP,ISTEP ENDDO ! END OF: DO JT=-ISTEP,ISTEP #ifdef NEVER fit(it,jt) = fit(it,jt) write(6,961) it,jt, fit(it,jt),cns2,cnr2, rotmp(it,jt) 961 format(' fit(',i2,',',i2,') : ',f12.4,' ',2f8.2,f6.1) #endif C FIND THE MAXIMUM CC ANGLE WITHIN +/-ISTEP C MAXIMUM CANNOT BE ON THE EDGE, I.E., IT,JT/=ISTEP AFIT = FIT(0,0) JTMA = 0 ITMA = 0 RANGNEWT = ROTMP(0,0) IF (ISTEP .GT. 1) THEN DO JT=-ISTEP+1,ISTEP-1 DO IT=-ISTEP+1,ISTEP-1 IF (FIT(IT,JT) .GT. AFIT) THEN AFIT = FIT(IT,JT) RANGNEWT = ROTMP(IT,JT) !compiler bug on OPT64 ITMA = IT JTMA = JT ENDIF ENDDO ! END OF: DO IT=-ISTEP,ISTEP ENDDO ! END OF: DO JT=-ISTEP,ISTEP ENDIF C TEMP VARIABLE OVERCOMES COMPILER BUG ON OPT 64 PGI 6.0 RANGNEW = RANGNEWT CCROTD = AFIT C COPY VALUES AROUND THE PEAK. DO JT=-1,1 DO IT=-1,1 FITP(IT,JT) = FIT(ITMA+IT,JTMA+JT) c write(6,910) it,jt,fitp(it,jt) 910 format(' fitp(',i5,',',i5,') : ',f14.4) ENDDO ENDDO C UPDATE INTEGRAL LOCATION OF PEAK ISX = ISX + ITMA ISY = ISY + JTMA SX = ISX SY = ISY #ifdef NEVER cnr2 = nrow / 2 + 1 + isy cns2 = nsam / 2 + 1 + isx write(6,905) idis,isx,isy, ccrotd,rangnew, cns2,cnr2,sx,sy 905 format(' 2 ',i5,' (',i3,',',i3,'): ',f12.4,' ',5f8.2) #endif C SUB-PIXEL INTERPOLATION ----------------------------------- CALL PARABLD(FITP,SSX,SSY,PEAK) IF (ABS(SSX) .LT. 1.0 .AND. ABS(SSY) .LT. 1.0) THEN C NOT ON EDGE OF 3x3 INTERPOLATION AREA CNS2 = NSAM/2+1 + SX + SSX CNR2 = NROW/2+1 + SY + SSY C NORMALIZE IMAGE VALUES UNDER THE MASK OVER VARIANCE RANGE C INTERPOLATE INTO POLAR COORDINATES, C CREATE FOURIER OF: CIRCEXP !write(6,*) ' aprings sub:,',sx,sy,cns2,cnr2 CALL APRINGS_ONE_NEW(NSAM,NROW, CNS2,CNR2, EXPBUF,.FALSE., & MODE,NUMR,NRING,LCIRC,WR,FFTW_PLANS, & CIRCEXP,IRTFLG) CALL CROSRNG_EP_NEW(CIRCREF(1,IBE),CIRCEXP,LCIRC,NRING, & MAXRIN,NUMR, CCROTD_INTERP,RANGNEW_INTERP, & TT,MIRRORED, & FFTW_PLANS(1)) #ifdef NEVER rt1 = ang_n(rangnew_interp,mode,maxrin) write(6,904) idis,isx,isy, ccrotd_interp,rt1, & cns2,cnr2, sx+ssx,sy+ssy 904 format(' 3 ',i5,' (',i3,',',i3,'): ',f12.4,' ',5f8.2) #endif IF (CCROTD_INTERP .GT. CCROTD) THEN C USE SUB-PIXEL LOCATION CCROTD = CCROTD_INTERP RANGNEW = ANG_N(RANGNEW_INTERP,MODE,MAXRIN) SX = SX + SSX SY = SY + SSY ENDIF ENDIF ENDIF CCROT = CCROTD SX = -SX SY = -SY C HAVE TO CHANGE ORDER OF SHIFT & ROTATION. C IN THIS PROGRAM IMAGE IS SHIFTED FIRST, ROTATED SECOND. C IN 'RT SQ' IT IS ROTATION FIRST, SHIFT SECOND. C THIS CODE CORRESPONDS TO 'SA P'. CO = COS(RANGNEW * DGR_TO_RAD) SO = -SIN(RANGNEW * DGR_TO_RAD) XSHSUM = SX*CO - SY*SO YSHSUM = SX*SO + SY*CO C ALMOST ZERO IS LIKELY TO BE ZERO IF (ABS(XSHSUM) .LT. 0.04) XSHSUM = 0.0 IF (ABS(YSHSUM) .LT. 0.04) YSHSUM = 0.0 IF (ABS(RANGNEW) .LT. 0.04) RANGNEW = 0.0 #ifdef NEVER cns2 = nsam / 2 + 1 - sx cnr2 = nrow / 2 + 1 - sy write(6,906) idis,isx,isy, ccrotd,rangnew, & cns2,cnr2, xshsum,yshsum 906 format(' 4 ',i5,' (',i3,',',i3,'): ',f12.4,' ',5f8.2) write(6,*) ' ------------------------------------- ' write(6,*) ' ' #endif 9999 IF (ALLOCATED(LCG)) DEALLOCATE(LCG) END