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 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 MRQLI_PS C C PURPOSE: FIND ROTATIONAL AND SHIFT PARAMETERS TO ALIGN A SERIES OF C REFERENCE IMAGES WITH SAMPLE IMAGES 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 REAL, ALLOCATABLE, DIMENSION(:,:) :: DLIST REAL, ALLOCATABLE,DIMENSION(:,:) :: REFDIR,EXPDIR REAL, ALLOCATABLE,DIMENSION(:,:) :: ANGREF,ANGEXP INTEGER, PARAMETER :: NLISTMAX = 15 REAL, DIMENSION(NLISTMAX) :: PARLIST PARAMETER (NLIST=7) 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 IBIGANGDIF = 0 ANGDIFAVG = 0.0 CCROTAVG = 0.0 IMPROVCCROT = 0 CCROTIMPROV = 0.0 IWORSECCROT = 0 CCROTWORSE = 0.0 LIMITRANGE = (RANGE .GT. 0.0 .AND. RANGE .LT. 360.0) RANGE = COS(RANGE*DGR_TO_RAD) NIMALCG = NUMREF 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, & LUNRING,SCRFILE,IRTFLG) NSAM = LSAM NROW = LROW ALLOCATE(A(LSAM,LROW,NUMTH), & DLIST(NLIST,NUMTH), & STAT=IRTFLG) IF (IRTFLG .NE. 0) THEN MWANT = LSAM*LROW*NUMTH + NLIST*NUMTH CALL ERRT(46,'A & DLIST',MWANT) RETURN ENDIF DLIST = 0.0 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(7,NUMEXP), EXPDIR(3,NUMEXP), STAT=IRTFLG) IF (IRTFLG .NE. 0) THEN MWANT = 10*NUMEXP CALL ERRT(46,'ANGEXP....',MWANT) RETURN ENDIF CALL AP_GETANGA(IEXPLIST,NUMEXP,0,EXPANGDOC,EXPPAT, & INPIC,INANG,7,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(7,1), EXPDIR(3,1), STAT=IRTFLG) IF (IRTFLG .NE. 0) THEN MWANT = 10 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),RANGE, & DLIST(2,IT),DLIST(3,IT), & DLIST(4,IT),DLIST(5,IT),DLIST(6,IT), & NIMALCG,CKMIRROR,LIMITRANGE,FFTW_PLANS) ENDDO C OUTPUT (IN DLIST POSITION IS INCREASED BY 1, NO.1 IS THE KEY). 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 C 6 - INPUT IMAGE NUMBER. DO IEXP=IEXPT,MIN(NUMEXP,IEXPT+NUMTH-1) IT = IEXP-IEXPT+1 IMGEXP = IEXPLIST(IEXP) C DLIST(2,IT IS LIST NUMBER OF MOST SIMILAR REF. IMAGE C (<0 IF MIRRORED, 0 IF NONE ) IREF = INT(DLIST(2,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(3,IT) RANGNEW = DLIST(4,IT) XSHNEW = DLIST(5,IT) YSHNEW = DLIST(6,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,LSAM,LROW,CCROT,PEAKV, & RANGNEW,XSHNEW,YSHNEW,MIRRORNEW,EXPPAT,REFPAT, & NIMALCG, CTYPE, A,LUNDOC,PARLIST) CCROTAVG = CCROTAVG + CCROT IF (NGOTPAR .GE. 8) THEN C COMPILE CCROT CHANGE STATISTICS ANGDIF = PARLIST(10) IF (ANGDIF .GT. ANGDIFTHR)IBIGANGDIF = IBIGANGDIF + 1 CCROTLAS = ANGEXP(8,IEXP) ANGDIFAVG = ANGDIFAVG + PARLIST(10) IF (CCROT .GE. CCROTLAS) THEN IMPROVCCROT = IMPROVCCROT + 1 CCROTIMPROV = CCROTIMPROV + CCROT ELSE IWORSECCROT = IWORSECCROT + 1 CCROTWORSE = CCROTWORSE + CCROT ENDIF ENDIF ! END OF: IF (NGOTPAR .GE. 8) 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(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 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 REAL, ALLOCATABLE, DIMENSION(:,:) :: DLIST REAL, ALLOCATABLE,DIMENSION(:,:) :: REFDIR,EXPDIR REAL, ALLOCATABLE,DIMENSION(:,:) :: ANGREF,ANGEXP INTEGER, PARAMETER :: NLISTMAX = 15 REAL, DIMENSION(NLISTMAX) :: PARLIST PARAMETER (NLIST=7) 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 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 IBIGANGDIF = 0 ANGDIFAVG = 0.0 CCROTAVG = 0.0 IMPROVCCROT = 0 CCROTIMPROV = 0.0 IWORSECCROT = 0 CCROTWORSE = 0.0 LIMITRANGE = (RANGE .GT. 0.0 .AND. RANGE .LT. 360.0) MAXRIN = NUMR(3,NRING) RANGE = COS(RANGE*DGR_TO_RAD) NIMALCG = NUMREF 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(7,NUMEXP), EXPDIR(3,NUMEXP), STAT=IRTFLG) IF (IRTFLG .NE. 0) THEN MWANT = 10*NUMEXP CALL ERRT(46,'ANGEXP....',MWANT) RETURN ENDIF CALL AP_GETANGA(IEXPLIST,NUMEXP,0,EXPANGDOC,EXPPAT, & INPIC,INANG,7,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(7,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),STAT=IRTFLG) IF (IRTFLG .NE. 0) THEN WRITE(6,*) 'MRQLI_PS: FAILED TO ALLOCATE PSIZE' RETURN ENDIF ALLOCATE(NBASE(NPROCS), STAT=IRTFLG) IF (IRTFLG .NE. 0) THEN WRITE(6,*) 'MRQLI_PS: FAILED TO ALLOCATE NBASE' RETURN ENDIF CALL SETPART(NUMEXP, PSIZE, NBASE) NEXPLOC = PSIZE(MYPID+1) ALLOCATE(ALOC(NSAM,NROW,NEXPLOC),STAT=IRTFLG) IF (IRTFLG .NE. 0) THEN WRITE(6,*) ' MRQLI_PS: FAILED TO ALLOCATE ALOC' RETURN END IF ALOC = 0.0 #ifdef MPI_DEBUG WRITE(6,111) NBASE(MYPID+1), MYPID CALL FLUSHFILE(6) 111 FORMAT(' MRQLI_PS: NBASE = ', I5, ' MYPID = ', I5) #endif C C === PROCESS 0 READS ALL IMAGES AND DISTRIBUTE C THEM TO DIFFERENT PROCESSORS. === ALLOCATE(ABUF(NSAM,NROW,PSIZE(1)),STAT=IRTFLG) IF (IRTFLG .NE. 0) THEN WRITE(6,*) ' MRQLI_PS: FAILED TO ALLOCATE ABUF' RETURN ENDIF #ifdef MPI_DEBUG TCOM0 = MPI_WTIME() #endif 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(DLISTLOC(7,NEXPLOC),STAT=IRTFLG) IF (IRTFLG .NE. 0) THEN WRITE(6,*) ' MRQLI_PS: FAILED TO ALLOCATE DLISTLOC' 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 ALLOCATE(PARTAB(15,NUMEXP), PARTABLOC(15,NEXPLOC), & STAT=IRTFLG) IF (IRTFLG .NE. 0) THEN CALL ERRT(46,'PARTAB',NUMEXP) RETURN 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),RANGE, & DLISTLOC(2,IT),DLISTLOC(3,IT), & DLISTLOC(4,IT),DLISTLOC(5,IT),DLISTLOC(6,IT), & NIMALCG,CKMIRROR,LIMITRANGE, FFTW_PLANS) C OUTPUT (IN DLIST POSITION IS INCREASED BY 1, NO.1 IS THE KEY). 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 C 6 - INPUT IMAGE NUMBER. IMGEXP = IEXPLIST(IGLB) C DLIST(2,IT IS LIST NUMBER OF MOST SIMILAR REF. IMAGE C (<0 IF MIRRORED, 0 IF NONE ) IREF = INT(DLISTLOC(2,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(3,IT) RANGNEW = DLISTLOC(4,IT) XSHNEW = DLISTLOC(5,IT) YSHNEW = DLISTLOC(6,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,LSAM,LROW,CCROT,PEAKV, & RANGNEW,XSHNEW,YSHNEW,MIRRORNEW,EXPPAT,REFPAT, & NIMALCG, CTYPE, A, 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(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 NIMALCG (OUTPUT) C C-********************************************************************** SUBROUTINE APRQ2D(EXPBUF,CIRCREF,TT,NUMR, & NSAM,NROW,ISHRANGE,ISTEP, & LCIRC,NRING,NUMREF,MODE, & REFDIR,EXPDIR,RANGE, & DIREF,CCROT,RANGNEW,XSHSUM,YSHSUM,NIMALCG, & 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 :: TOTA,TMTA 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) ELSE C RESTRICTED RANGE SEARCH ALLOCATE(LCG(NUMREF),STAT=IRTFLG) IF (IRTFLG .NE. 0) THEN CALL ERRT(46,'LCG',NUMREF) RETURN ENDIF C write(6,*) 'range: ',range,' nima:',nima C MAKE LIST OF NEARBY REFERENCE IMAGES NIMALCG = 0 DO IMI=1,NUMREF 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 MIRRORED OR NON-MIRRORED IS WITHIN RANGE IF (CKMIRROR .OR. DT .GT. 0) THEN C DO NOT DISCARD IF NOT MIRRORED OR WANT MIRRORED NIMALCG = NIMALCG + 1 LCG(NIMALCG) = IMI IF (DT .LT. 0) LCG(NIMALCG) = -IMI C write(6,*) 'imi: ',imi,' dt: ',dt,' dtabs: ',dtabs, C & refdir(1,imi),refdir(2,imi),refdir(3,imi) ENDIF ENDIF ENDDO IF (NIMALCG .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 = NIMALCG 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 COORDINATES, CREATE FOURIER OF: A_CIRC C NO WEIGHTING OF RINGS 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, & TOTA,TOT, & 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, & TOTA,TOT, TMTA,TMT, & 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 (TOTA .GE. CCROTD) THEN C GOOD MATCH WITH TOTA (MIRRORED OR NOT) POSITION CCROTD = TOTA IBE = IR ISX = IT ISY = JT RANGNEW = ANG_N(TOT,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 (TMTA .GE. CCROTD) THEN C GOOD MATCH WITH MIRRORED POSITION CCROTD = TMTA IBE = IR ISX = IT ISY = JT RANGNEW = ANG_N(TMT,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 DEBUG 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 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 DEBUG 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 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 DEBUG 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 DEBUG 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