C++********************************************************************* C C MRQLI_SS.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 'OR MQ' SUPPORT SEP 03 ARDEAN LEITH C AP_OUT USAGE FEB 04 ARDEAN LEITH C 'MQ' NOT READ EXP. ANGLES APR 04 ARDEAN LEITH C BAD PEAK IF INTERP ON EDGE AUG 04 ARDEAN LEITH C CROSRNG_E SPEEDS UP AUG 04 ARDEAN LEITH C AP_END CALL HAS PARLIST OCT 04 ARDEAN LEITH C LIMITRANGE BUG OCT 04 ARDEAN LEITH C PEAKV = 1 JAN 05 ARDEAN LEITH C DISCARD MIRROR ... JUN 06 ARDEAN LEITH C AP_STAT CALL JAN 07 ARDEAN LEITH C USE FFTW3 IN APRINGS MAR 08 ARDEAN LEITH C APRINGS_ONE_NEW APR 08 ARDEAN LEITH C AP_END CALL ALTERED NOV 08 ARDEAN LEITH C AP_STAT_ADD NOV 08 ARDEAN LEITH C CIRCREF NOT INCORE AUG 09 ARDEAN LEITH C 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_SS(IREFLIST,NUMREF,IEXPLIST,NUMEXP, C LSAM,LROW,NR,LENTT,ISHRANGE,ISTEP, C NRING,LCIRC,NUMR,CIRCREF,CIRCREF_IN_CORE, C MODE,REFANGDOC,EXPANGDOC,SCRFILE,FFTW_PLANS, C REFPAT,EXPPAT,RANGE,CKMIRROR,CTYPE,LUNDOC) C C PURPOSE: FIND ROTATIONAL AND SHIFT PARAMETERS TO ALIGN A SERIES OF C REFERENCE IMAGES WITH SERIES OF SAMPLE IMAGES C C VERSION FOR MP AND A SMALL NUMBER OF IMAGES TO BE ALIGNED. C NOT FOR MPI USE. C C OPERATIONS: 'AP SH' C C PARAMETERS: C IREFLIST LIST OF REF. IMAGE FILE NUMBERS (INPUT) C NUMREF NO. OF REF. IMAGES (INPUT) C IEXPLIST LIST OF EXP. IMAGE FILE NUMBERS (INPUT) C NUMEXP NO. OF EXP 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 C23456789012345678901234567890123456789012345678901234567890123456789012 C--********************************************************************* SUBROUTINE MRQLI_SS(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 'CMBLOCK.INC' INCLUDE 'CMLIMIT.INC' INTEGER, DIMENSION(NUMREF) :: IREFLIST INTEGER, DIMENSION(NUMEXP) :: IEXPLIST INTEGER, DIMENSION(3,NRING) :: NUMR REAL, DIMENSION(LCIRC,NUMREF) :: CIRCREF LOGICAL :: CIRCREF_IN_CORE,CKMIRROR LOGICAL :: MIRRORNEW LOGICAL :: GOTREFANG LOGICAL :: LIMITRANGE CHARACTER(LEN=1) :: MODE,NULL CHARACTER (LEN=*) :: REFANGDOC,EXPANGDOC CHARACTER (LEN=*) :: SCRFILE CHARACTER (LEN=*) :: REFPAT,EXPPAT CHARACTER (LEN=*) :: CTYPE DOUBLE PRECISION :: FITP(-1:1,-1:1) DOUBLE PRECISION :: CCROTD,PEAK,CCROTD_INTERP INTEGER *8 :: FFTW_PLANS(*) #ifdef __ia64 C INTEL PARALLEL COMPILER BUGGY MUST ALLOCATE THESE ARRAYS IF || C DOES NOT HANDLE NEGATIVE SUBSCRIPTS CORRECTLY DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:,:,:) :: CCOA,CCMA REAL, ALLOCATABLE, DIMENSION(:,:,:) :: RANGOA,RANGMA DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:,:) :: FIT REAL, ALLOCATABLE, DIMENSION(:,:) :: ROTMP #else C AUTOMATIC ARRAYS DOUBLE PRECISION CCOA(NUMREF,-ISHRANGE:ISHRANGE, & -ISHRANGE:ISHRANGE) DOUBLE PRECISION CCMA(NUMREF,-ISHRANGE:ISHRANGE, & -ISHRANGE:ISHRANGE) DIMENSION RANGOA(NUMREF, -ISHRANGE:ISHRANGE, & -ISHRANGE:ISHRANGE), & RANGMA(NUMREF, -ISHRANGE:ISHRANGE, & -ISHRANGE:ISHRANGE) DOUBLE PRECISION FIT( -ISTEP:ISTEP,-ISTEP:ISTEP) DIMENSION ROTMP(-ISTEP:ISTEP,-ISTEP:ISTEP) #endif DOUBLE PRECISION, DIMENSION(LENTT) :: TT INTEGER, DIMENSION(NUMREF) :: LCG REAL, DIMENSION(3) :: ANGOUT REAL, DIMENSION(3) :: EXPDIR LOGICAL :: MIRRORED C ALLOCATED ARRAYS REAL, ALLOCATABLE, DIMENSION(:,:,:) :: CIRCEXP REAL, ALLOCATABLE, DIMENSION(:,:) :: A REAL, ALLOCATABLE,DIMENSION(:,:) :: REFDIR REAL, ALLOCATABLE,DIMENSION(:,:) :: ANGREF,ANGEXP INTEGER, PARAMETER :: NLISTMAX = 15 REAL, DIMENSION(NLISTMAX) :: PARLIST PARAMETER (QUADPI = 3.1415926535897932384626) PARAMETER (DGR_TO_RAD = (QUADPI/180)) DATA INPIC/77/,INANG/78/,LUNRING/50/ NULL = CHAR(0) C INITIALIZE CCROT STATISTICS COUNTERS ANGDIFTHR = 0.0 CALL AP_STAT_ADD(-1,CCROT,ANGDIF,ANGDIFTHR,CCROTLAS, & CCROTAVG,IBIGANGDIF,ANGDIFAVG,IMPROVCCROT, & CCROTIMPROV,IWORSECCROT,CCROTWORSE) C FLAG FOR RESTRICTED PROJECTION RANGE LIMITRANGE = RANGE .GT. 0 .AND. RANGE .LT. 360 MAXRIN = NUMR(3,NRING) #ifdef SP_LIBFFTW3 MAXRIN = NUMR(3,NRING) -2 ! ACTUAL LENGTH OF LONGEST RING #endif C THIS ALTERS RANGE! RANGE = COS(RANGE*DGR_TO_RAD) WR = 0.0 ! DUMMY VALUE FLAG FOR APRINGS CALL 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 OR C CREATE REFERENCE RINGS FILE FOR LATER READING 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 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 NSIS = MAX(ISHRANGE/ISTEP, ISTEP) ALLOCATE(CIRCEXP(LCIRC,-NSIS:NSIS,-NSIS:NSIS), & A(NSAM,NROW), & STAT=IRTFLG) IF (IRTFLG .NE. 0) THEN NSIST = 2 * NSIS + 1 MWANT = LCIRC*NSIST*NSIST + NSAM*NROW CALL ERRT(46,' CIRCEXP,....',MWANT) GOTO 9999 ENDIF NUMREFLCG = NUMREF IEND = NUMREF NGOTPAR = 0 IF (LIMITRANGE .OR. CTYPE(1:2) .EQ. 'SH') THEN C RESTRICTED ANGULAR SEARCH FOR 'RQ' OR 'SH' ALLOCATE(ANGREF(3,NUMREF), STAT=IRTFLG) IF (IRTFLG .NE. 0) THEN MWANT = 3*NUMREF CALL ERRT(46,'ANGREF',MWANT) RETURN ENDIF C LOAD REF. PROJ. ANGLES () FROM DOC. FILE (REFANGDOC) OR C IMAGE FILE (REFPAT) HEAD INTO ANGREF CALL AP_GETANGA(IREFLIST,NUMREF,0,REFANGDOC,REFPAT, & INPIC,INANG,3,ANGREF,NGOTREF,IRTFLG) IF (IRTFLG .NE. 0) GOTO 9999 ENDIF IF (CTYPE(1:2) .NE. 'MQ' .AND. & EXPANGDOC .NE. CHAR(0)) THEN C READ EXP. ANGLES INTO ANGEXP IF (REFANGDOC .NE. NULL) THEN ALLOCATE(REFDIR(3,NUMREF),STAT=IRTFLG) IF (IRTFLG .NE. 0) THEN MWANT = 3*NUMREF CALL ERRT(46,'REFDIR',MWANT) RETURN ENDIF C CONVERT REF. ANGLES TO UNITARY DIRECTIONAL VECTORS (REFDIR). CALL AP_GETSATA(ANGREF,REFDIR,3,NUMREF,IRTFLG) ENDIF C LOAD EXP. PROJ. ANGLES & ALIGNMENT PARAMETERS (ANGEXP) C FROM DOC. FILE (EXPANGDOC) OR IMAGE FILE (REFPAT) HEAD ALLOCATE(ANGEXP(8,NUMEXP), STAT=IRTFLG) IF (IRTFLG .NE. 0) THEN MWANT = 8*NUMEXP CALL ERRT(46,'ANGEXP',MWANT) RETURN ENDIF C THIS RETURNS NGOTPAR CALL AP_GETANGA(IEXPLIST,NUMEXP,0,EXPANGDOC,EXPPAT, & INPIC,INANG,8,ANGEXP,NGOTPAR,IRTFLG) IF (IRTFLG .NE. 0) GOTO 9999 GOTREFANG = .TRUE. ENDIF C LOOP OVER EXP. IMAGES TO BE ALIGNED DO IEXP=1,NUMEXP IMGEXP = IEXPLIST(IEXP) IF (CTYPE(1:2) .EQ. 'RQ' .OR. & (LIMITRANGE .AND. EXPANGDOC .NE. NULL)) THEN C CONVERT EXP. ANGLE TO UNITARY DIRECTIONAL VECTORS (EXPDIR). CALL AP_GETSATA(ANGEXP(1,IEXP),EXPDIR,8,1,IRTFLG) C MAKE LIST OF NEARBY REFERENCE IMAGES NUMREFLCG = 0 DO IMIR=1,NUMREF C LOOP OVER REFERENCE IMAGES C DT NEAR 1.0 = NOT-MIRRORED, DT NEAR -1.0 = MIRRORED DT = (EXPDIR(1) * REFDIR(1,IMIR) + & EXPDIR(2) * REFDIR(2,IMIR) + & EXPDIR(3) * REFDIR(3,IMIR)) 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 NUMREFLCG = NUMREFLCG + 1 LCG(NUMREFLCG) = IMIR IF (DT .LT. 0) LCG(NUMREFLCG) = -IMIR ENDIF IEND = NUMREFLCG ENDIF ENDDO IF (NUMREFLCG .LE. 0) THEN C REPORT THAT THERE ARE NO NEARBY REFERENCE IMAGES IMGREF = 0 PEAKV = 0.0 CCROT = 0.0 XSHNEW = 0.0 YSHNEW = 0.0 MIRRORNEW = .FALSE. CALL AP_END(IEXP,IMGEXP,IMGREF, & ANGREF,REFDIR, & ANGEXP(1,IEXP),EXPDIR,ISHRANGE, & GOTREFANG, NGOTPAR, CCROT,PEAKV, & RANGNEW,XSHNEW,YSHNEW,MIRRORNEW,REFPAT, & NUMREFLCG, CTYPE, LUNDOC,PARLIST) CYCLE ENDIF ENDIF C LOAD EXP. IMAGE DATA FOR THIS IMAGE IMITT = IEXP IF (CTYPE(1:2) .EQ. 'OR') IMITT = 0 CALL AP_GETDATS(IEXPLIST,NUMEXP,LSAM,LROW, & NUMTH,EXPPAT,INPIC,IMITT,IMITT, & A, AVI,SIGI, IRTFLG) IF (IRTFLG .NE. 0) GOTO 9999 #ifdef __ia64 ALLOCATE(CCOA(NUMREF,-ISHRANGE:ISHRANGE, -ISHRANGE:ISHRANGE), & CCMA(NUMREF,-ISHRANGE:ISHRANGE, -ISHRANGE:ISHRANGE), & RANGOA(NUMREF, -ISHRANGE:ISHRANGE, -ISHRANGE:ISHRANGE), & RANGMA(NUMREF, -ISHRANGE:ISHRANGE, -ISHRANGE:ISHRANGE), & FIT( -ISTEP:ISTEP, -ISTEP:ISTEP), & ROTMP(-ISTEP:ISTEP, -ISTEP:ISTEP), & STAT=IRTFLG) IF (IRTFLG .NE. 0) THEN MWANT1 = 4 * (NUMREF * (ISHRANGE*2+1) * (ISHRANGE*2+1)) + & 2 * (ISTEP*2+1) CALL ERRT(46,' CCOA,....',MWANT1) GOTO 9999 ENDIF #endif IF ( CIRCREF_IN_CORE) THEN C USE CIRCREF ARRAY FOR REFERENCE RINGS !write(6,*) ' incore, old parallel ' C LOOP OVER ALL SHIFTS IN || c$omp parallel do private(jt,it,imil,imi,cnr2,cns2,mirrored) DO JT=-ISHRANGE,ISHRANGE,ISTEP C LOOP OVER SHIFTED CENTERS IN Y CNR2 = NROW/2+1+JT DO IT=-ISHRANGE,ISHRANGE,ISTEP C LOOP OVER SHIFTED CENTERS IN X CNS2 = NSAM/2+1+IT C NORMALIZE IMAGE VALUES UNDER THE MASK OVER VARIANCE C RANGE INTERPOLATE TO POLAR COORDINATES, CREATE C FOURIER OF: CIRCEXP. NO WEIGHTING OF RINGS CALL APRINGS_ONE_NEW(NSAM,NROW, CNS2,CNR2, & A,.FALSE., & MODE,NUMR,NRING,LCIRC,WR,FFTW_PLANS, & CIRCEXP(1,IT/ISTEP,JT/ISTEP),IRTFLG) DO IMIL=1,IEND C LOOP OVER REFERENCE IMAGES 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(1,IT/ISTEP,JT/ISTEP), & LCIRC,NRING,MAXRIN,NUMR, & CCOA(IMI,IT,JT),RANGOA(IMI,IT,JT), & TT,MIRRORED, FFTW_PLANS(1)) ELSE C CHECK BOTH NON-MIRRORED & MIRRORED POSITIONS CALL CROSRNG_MSP_NEW(CIRCREF(1,IMI), & CIRCEXP(1,IT/ISTEP,JT/ISTEP), & LCIRC,NRING,MAXRIN,NUMR, & CCOA(IMI,IT,JT),RANGOA(IMI,IT,JT), & CCMA(IMI,IT,JT),RANGMA(IMI,IT,JT), & TT,FFTW_PLANS(1)) ENDIF ENDDO ! END OF: DO IMIL=1,IEND ENDDO ! END OF: DO IT=-ISHRANGE,ISHRANGE,ISTEP ENDDO ! END OF: DO JT=-ISHRANGE,ISHRANGE,ISTEP c$omp end parallel do C END OF THE OMP PARALLEL SECTION ELSE !write(6,*) ' not incore, new parallel ' C USE REFERENCE RINGS FILE (MIGHT BE AN INCORE FILE?) ----- DO JT=-ISHRANGE,ISHRANGE,ISTEP C LOOP OVER SHIFTED CENTERS IN Y CNR2 = NROW/2+1+JT DO IT=-ISHRANGE,ISHRANGE,ISTEP C LOOP OVER SHIFTED CENTERS IN X CNS2 = NSAM/2+1+IT C NORMALIZE IMAGE VALUES UNDER THE MASK OVER VARIANCE C RANGE INTERPOLATE TO POLAR COORDINATES, CREATE C FOURIER OF: CIRCEXP. NO WEIGHTING OF RINGS CALL APRINGS_ONE_NEW(NSAM,NROW, CNS2,CNR2, & A,.FALSE., & MODE,NUMR,NRING,LCIRC,WR,FFTW_PLANS, & CIRCEXP(1,IT/ISTEP,JT/ISTEP),IRTFLG) c$omp parallel do private(imil,imi,ithread,mirrored) c$omp& schedule(static,1) DO IMIL=1,IEND C LOOP OVER ALL REFERENCE IMAGES IN || 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(1,IT/ISTEP,JT/ISTEP), & LCIRC,NRING,MAXRIN,NUMR, & CCOA(IMI,IT,JT),RANGOA(IMI,IT,JT), & TT,MIRRORED, FFTW_PLANS(1)) ELSE C CHECK BOTH NON-MIRRORED & MIRRORED POSITIONS CALL CROSRNG_MSP_NEW(CIRCREF(1,ITHREAD), & CIRCEXP(1,IT/ISTEP,JT/ISTEP), & LCIRC,NRING,MAXRIN,NUMR, & CCOA(IMI,IT,JT),RANGOA(IMI,IT,JT), & CCMA(IMI,IT,JT),RANGMA(IMI,IT,JT), & TT,FFTW_PLANS(1)) ENDIF ENDDO ! END OF: DO IMIL=1,IEND ENDDO ! END OF: DO IT=-ISHRANGE,ISHRANGE,ISTEP ENDDO ! END OF: DO JT=-ISHRANGE,ISHRANGE,ISTEP ENDIF C LOCATE BEST CCROT MATCH FROM ALL THE VALUES ACCUMULATED ABOVE CCROTD = -1.0D23 DO JT=-ISHRANGE,ISHRANGE,ISTEP DO IT=-ISHRANGE,ISHRANGE,ISTEP DO IRR=1,IEND C LOOP OVER REFERENCE IMAGES IR = IRR IF (LIMITRANGE) IR = ABS(LCG(IRR)) IF (CCOA(IR,IT,JT) .GE. CCROTD) THEN C BETTER POSITION (MAY BE MIRRORED IF LIMITRANGE) CCROTD = CCOA(IR,IT,JT) IREF = IR ISX = IT ISY = JT RANGNEW = ANG_N(RANGOA(IR,IT,JT),MODE,MAXRIN) MIRRORNEW = (LIMITRANGE .AND. (LCG(IRR) .LT. 0)) ENDIF IF (CKMIRROR .AND. .NOT. LIMITRANGE) THEN C HAVE TO COMPARE WITH MIRRORED POSITION IF (CCMA(IR,IT,JT) .GE. CCROTD) THEN C BETTER MIRRORED POSITION IREF = IR CCROTD = CCMA(IR,IT,JT) ISX = IT ISY = JT RANGNEW = ANG_N(RANGMA(IR,IT,JT),MODE,MAXRIN) MIRRORNEW = .TRUE. ENDIF ENDIF ENDDO ENDDO ENDDO C write(6,*) 'iref: ',iref, ' ccrotd: ',ccrotd SX = ISX SY = ISY CCROT = CCROTD IMGREF = IREFLIST(IREF) IF (CIRCREF_IN_CORE) THEN IREFT = IREF ELSE C FILL CIRCREF (1,1) FROM REFERENCE RINGS FILE IREFT = 1 CALL REDLIN(LUNRING,CIRCREF(1,IREFT),LCIRC,IREF) ENDIF #ifdef DEBUG write(6,921) imgref,isx,isy,ccrotd,rangnew 921 format(' 1 ',i5,' (',i3,',',i3,'): ',f14.4,' ',2f8.2,f6.1) #endif C CHECK LOCATIONS WITHIN ISHRANGE AROUND MAX ------------------ IF (ABS(ISX) .NE. ISHRANGE .AND. & ABS(ISY) .NE. ISHRANGE) THEN C NOT ON BOUNDARY, HAVE TO FIND NEIGHBOURING VALUES FIT(0,0) = CCROTD ROTMP(0,0) = RANGNEW c$omp parallel do private(jt,it,cnr2,cns2) DO JT=-ISTEP,ISTEP DO IT=-ISTEP,ISTEP CNR2 = NROW / 2 + 1 + JT + ISY IF (IT.NE.0 .OR. JT.NE.0) THEN CNS2 = NSAM / 2 + 1 + IT + ISX 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, & A,.FALSE., & MODE,NUMR,NRING,LCIRC,WR,FFTW_PLANS, & CIRCEXP(1,IT,JT),IRTFLG) CALL CROSRNG_EP_NEW(CIRCREF(1,IREFT), & CIRCEXP(1,IT,JT), & LCIRC,NRING,MAXRIN,NUMR, & FIT(IT,JT),ROTMP(IT,JT), & TT, MIRRORNEW,FFTW_PLANS(1)) C RECORD BEST ANGLE IN ROTMP ROTMP(IT,JT) = ANG_N(ROTMP(IT,JT),MODE,MAXRIN) ENDIF ENDDO ! END OF: DO IT=-ISTEP,ISTEP ENDDO ! END OF: DO JT=-ISTEP,ISTEP c$omp end parallel do 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 ENDDO ENDIF ! END OF: IF (ISTEP .GT. 1) C write(6,*) ((fit(i,j),i=-1,1),j=-1,1) C TEMP VARIABLE OVERCOMES COMPILER BUG ON OPT 64 PGI 6.0 RANGNEW = RANGNEWT 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 THE PEAK CCROTD = AFIT 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)imgref,isx,isy,ccrotd,rangnew,cns2,cnr2,sx,sy 905 format(' 2 ',i5,' (',i3,',',i3,'): ',f12.4,' ',5f8.2) #endif C SUB-PIXEL INTERPOLATION ------------------------------ C FIND PEAK BY FITTING PARABOLA TO 3x3 NEIGHBORHOOD CALL PARABLD(FITP,SSX,SSY,PEAK) IF (ABS(SSX) .LT. 1.0 .AND. ABS(SSY) .LT. 1.0) THEN C INTERPOLATED LOCATION IS NOT ON BOUNDARY 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, A,.FALSE., & MODE,NUMR,NRING,LCIRC,WR,FFTW_PLANS, & CIRCEXP,IRTFLG) CALL CROSRNG_EP_NEW(CIRCREF(1,IREFT),CIRCEXP,LCIRC, & NRING,MAXRIN,NUMR,CCROTD_INTERP, & RANGNEW_INTERP,TT,MIRRORNEW,FFTW_PLANS(1)) #ifdef DEBUG rt1 = ang_n(rangnew_interp,mode,maxrin) write(6,904) imgref,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 ! END OF: IF (ABS(SX) .LT. 1.0 ..... ENDIF ! END OF: IF ( ABS(ISX) .NE. ISHRANGE ...... 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) XSHNEW = SX*CO - SY*SO YSHNEW = SX*SO + SY*CO PEAKV = 1.0 C ALMOST ZERO IS LIKELY TO BE ZERO IF (ABS(XSHNEW) .LT. 0.04) XSHNEW = 0.0 IF (ABS(YSHNEW) .LT. 0.04) YSHNEW = 0.0 IF (ABS(RANGNEW) .LT. 0.04) RANGNEW = 0.0 CALL AP_END(IEXP,IMGEXP,IMGREF, & ANGREF,REFDIR(1,IREF), & ANGEXP(1,IEXP),EXPDIR,ISHRANGE, & GOTREFANG, NGOTPAR, CCROT,PEAKV, & RANGNEW,XSHNEW,YSHNEW,MIRRORNEW,REFPAT, & NUMREFLCG, CTYPE, LUNDOC,PARLIST) CALL AP_STAT_ADD(NGOTPAR,CCROT,PARLIST(10), & ANGDIFTHR,ANGEXP(8,IEXP), & CCROTAVG,IBIGANGDIF,ANGDIFAVG,IMPROVCCROT, & CCROTIMPROV,IWORSECCROT,CCROTWORSE) #ifdef DEBUG cns2 = nsam / 2 + 1 - sx cnr2 = nrow / 2 + 1 - sy write(6,906) imgref,isx,isy, ccrotd,rangnew, & cns2,cnr2, xshsum,yshsum 906 format(' 4 ',i5,' (',i3,',',i3,'): ',f12.4,' ',5f8.2) write(6,*) ' ------------------------------------- ' write(6,*) ' ' #endif ENDDO CALL REG_GET_USED(NSEL_USED) IF (NSEL_USED .GT. 0 .AND. CTYPE(1:2) .EQ. 'OR') THEN C OUTPUT TO REGISTER NOT TO DOC FILE (FOR 'OR') DMR = 0 IF (MIRRORNEW) DMR = 1 CALL REG_SET_NSEL(1,5,RANGNEW,XSHNEW,YSHNEW, & DMR,CCROT,IRTFLG) ENDIF IF (LUNDOC .GT. 0 .AND. 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 IF (ALLOCATED(CIRCEXP)) DEALLOCATE(CIRCEXP) IF (ALLOCATED(A)) DEALLOCATE(A) IF (ALLOCATED(ANGREF)) DEALLOCATE(ANGREF) IF (ALLOCATED(ANGEXP)) DEALLOCATE(ANGEXP) IF (ALLOCATED(REFDIR)) DEALLOCATE(REFDIR) END