C++******************************************************************** C C DSGR_PM.F C ADDED APSHIFT NOV 03 ARDEAN LEITH C AP_END FEB 04 ARDEAN LEITH C CROSRNG_E SPEEDS UP AUG 04 ARDEAN LEITH C AP_END CALL HAS PARLIST OCT 04 ARDEAN LEITH C CCROT STATISTICS FEB 05 ARDEAN LEITH C APRINGS_ONE MAR 05 ARDEAN LEITH C CROSRNG_E REWRITE MAR 08 ARDEAN LEITH C IREFLIST VAR. NAME MAR 08 ARDEAN LEITH C FMRS_PLAN MAR 08 ARDEAN LEITH C APRINGS_ONE_NEW MAR 08 ARDEAN LEITH C ********************************************************************** 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_PM(IREFLIST,NUMREF,IEXPLIST,NUMEXP, C NSAM,NROW,NR,LENTT,RANGE, C NRING,LCIRC,NUMR,CIRCREF, 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_E. C SINCE Y WERE PRE-WEIGHTED THE RESULTS IS ALREADY CORRECT. C C PARAMETERS: C IREFLIST LIST OF REF. IMAGE FILE NUMBERS (INPUT) C NUMREF NO. OF IMAGES (INPUT) C IEXPLIST LIST OF EXP. IMAGE FILE NUMBERS (INPUT) C NUMEXP NO. OF IMAGES (INPUT) C NSAM,NROW ACTUAL (NOT-WINDOWED) IMAGE SIZE (INPUT) C REFANGDOC REF. ANGLES FILE NAME (INPUT) C EXPANGDOC EXP. ANGLES FILE NAME (INPUT) C REFPAT REF. IMAGE SERIES FILE TEMPLATE (INPUT) C EXPPAT EXP. IMAGE SERIES FILE TEMPLATE (INPUT) C C OPERATIONS: 'AP MD', 'AP REF', 'AP RD', 'AP RN' C C23456789 123456789 123456789 123456789 123456789 123456789 123456789 12 C--********************************************************************* SUBROUTINE DSGR_PM(IREFLIST,NUMREF,IEXPLIST,NUMEXP, & NSAM,NROW,ANGDIFTHR,LENTT,RANGE, & NRING,LCIRC,NUMR,CIRCREF, & MODE, REFANGDOC,EXPANGDOC,SCRFILE,FFTW_PLANS, & REFPAT,EXPPAT,CKMIRROR,CTYPE,ISHRANGE,LUNDOC) #ifdef SP_LIBFFTW3 USE TYPE_KINDS INTEGER(KIND=I_8) :: IPLAN !STRUCTURE POINTER #endif INCLUDE 'CMLIMIT.INC' INCLUDE 'CMBLOCK.INC' INTEGER, DIMENSION(NUMREF) :: IREFLIST INTEGER, DIMENSION(NUMEXP) :: IEXPLIST INTEGER, DIMENSION(3,NRING) :: NUMR REAL, DIMENSION(LCIRC,NUMREF) :: CIRCREF CHARACTER(LEN=1) :: MODE,NULL CHARACTER (LEN=*) :: REFPAT,EXPPAT,SCRFILE CHARACTER (LEN=*) :: REFANGDOC,EXPANGDOC LOGICAL :: CKMIRROR CHARACTER (LEN=*) :: CTYPE C ALLOCATABLE ARRAYS REAL, ALLOCATABLE, DIMENSION(:,:,:) :: XBUF REAL, ALLOCATABLE, DIMENSION(:) :: CCLIST REAL, ALLOCATABLE, DIMENSION(:) :: RANGNLIST INTEGER, ALLOCATABLE, DIMENSION(:) :: NPROJ INTEGER, ALLOCATABLE, DIMENSION(:) :: ILIST LOGICAL, ALLOCATABLE, DIMENSION(:) :: MIRLIST !integer, allocatable, dimension(:) :: iptr LOGICAL :: CIRCREF_IN_CORE LOGICAL :: MIRRORNEW LOGICAL :: GOTREFANG DOUBLE PRECISION :: TOTMIN_DUM C AUTOMATIC ARRAYS REAL, DIMENSION(6) :: DLIST DOUBLE PRECISION, DIMENSION(LENTT) :: TT REAL, DIMENSION(8,NUMEXP) :: ANGEXP REAL, DIMENSION(3,NUMEXP) :: EXPDIR REAL, DIMENSION(3,NUMREF) :: REFDIR,ANGREF 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/ NULL = CHAR(0) C FIND NUMBER OF OMP THREADS CALL GETTHREADS(NUMTH) #ifdef SP_LIBFFT IF (LENTT .GT. 15) THEN C CREATE TT FOR SGI FFT USE CALL DZFFT1DUI(LENTT-15,TT) ENDIF #endif #ifdef DEBUGD NR = NUMR(1,NRING) ! OUTER RING NUMBER MAXRIN = NUMR(3,NRING) write(88,*) ' DEBUG OUTPUT FROM DSGR_PM' write(88,902) NR,MAXRIN 902 format(' OUTER RING NUMBER: ',I10,' RING LENGTH:',I10) do i = 1, nring write(88,901) NUMR(1,I),NUMR(2,I),NUMR(3,I) 901 format(3i10) enddo #endif C THIS ALTERS RANGE! RANGE = COS(RANGE*DGR_TO_RAD) C MAKE XBUF BIG ENOUGH FOR ALL THREADS MWANTX = NSAM * NROW * NUMTH ALLOCATE(XBUF(NSAM,NROW,NUMTH), & ILIST(NUMTH), & CCLIST(NUMTH), & RANGNLIST(NUMTH), & NPROJ(NUMTH), & MIRLIST(NUMTH), & STAT=IRTFLG) IF (IRTFLG .NE. 0) THEN MWANT = MWANTX + 5*NUMTH CALL ERRT(46,'XBUF...',MWANT) GOTO 9999 ENDIF C DUMMY CALL TO INITIALIZE REV. FFTW3 PLAN FOR USE WITHIN OMP || CALL FMRS_PLAN(.FALSE.,DUM,NSAM,NROW,1, 1, -1, IPLAN, IRTFLG) CALL FMRS_PLAN(.FALSE.,DUM,NSAM,NROW,1, 1, +1, IPLAN, IRTFLG) C ZERO DLIST ARRAY DLIST = 0.0 GOTREFANG = .FALSE. IF (CTYPE(1:2) .NE. 'MD') THEN C GET REF. PROJ. ANGLES (ANGREF) FROM DOC. FILE (REFANGDOC) OR C IMAGE FILE (REFPAT) HEAD CALL AP_GETANGA(IREFLIST,NUMREF,0,REFANGDOC,REFPAT, & INPIC,INANG,3,ANGREF,NGOTREF,IRTFLG) IF (IRTFLG .NE. 0) GOTO 9999 GOTREFANG = .TRUE. C CONVERT REF. ANGLES TO UNITARY DIRECTIONAL VECTORS (REFDIR). CALL AP_GETSATA(ANGREF,REFDIR,3,NUMREF,IRTFLG) IF (IRTFLG .NE. 0) GOTO 9999 ENDIF NGOTPAR = 0 IF (EXPANGDOC .NE. NULL) THEN C LOAD EXP. PROJ. ANGLES & ALIGNMENT PARAMETERS (ANGEXP) C FROM DOC. FILE (EXPANGDOC) CALL AP_GETANGA(IEXPLIST,NUMEXP,0,EXPANGDOC,EXPPAT, & INPIC,INANG,8,ANGEXP,NGOTPAR,IRTFLG) IF (IRTFLG .NE. 0) GOTO 9999 C CONVERT EXP. ANGLES TO UNITARY DIRECTIONAL VECTORS (EXPDIR). CALL AP_GETSATA(ANGEXP,EXPDIR,8,NUMEXP,IRTFLG) IF (IRTFLG .NE. 0) GOTO 9999 ELSE ANGEXP = 0.0 ENDIF C READ REFERENCE IMAGES INTO REFERENCE RINGS (CIRCREF) ARRAY CIRCREF_IN_CORE = .TRUE. CALL APRINGS_NEW(IREFLIST,NUMREF, NSAM,NROW, & NRING,LCIRC,NUMR,MODE,FFTW_PLANS, & REFPAT,INPIC,CIRCREF,CIRCREF_IN_CORE, & LUNRING,SCRFILE,IRTFLG) IF (IRTFLG .NE. 0) GOTO 9999 C INITIALIZE CCROT STATISTICS CCROTAVG = 0.0 IMPROVCCROT = 0 CCROTIMPROV = 0.0 IWORSECCROT = 0 CCROTWORSE = 0.0 ANGDIFAVG = 0.0 IBIGANGDIF = 0 DO IEXPT=1,NUMEXP,NUMTH C LOOP OVER ALL EXPERIMENTAL (SAMPLE) IMAGES C LOAD NUMTH EXPERIMENTAL IMAGES INTO ARRAY XBUF IEND = MIN(NUMEXP,IEXPT+NUMTH-1) CALL AP_GETDATS(IEXPLIST,NUMEXP,NSAM,NROW, & NUMTH,EXPPAT,INPIC, IEXPT,IEND, & XBUF, IRTFLG) IF (IRTFLG .NE. 0) GOTO 9999 c$omp parallel do private(IEXP,IT) DO IEXP=IEXPT,IEND IT = IEXP-IEXPT+1 CALL DSGR_2(XBUF(1,1,IT),NSAM,NROW,NUMR,NRING, & MODE,CIRCREF,LCIRC,NUMREF,TT,LENTT, & REFDIR,EXPDIR(1,IEXP),RANGE, & ILIST(IT),CCLIST(IT),RANGNLIST(IT),MIRLIST(IT), & NPROJ(IT),CKMIRROR,FFTW_PLANS) ENDDO C LOOP OVER CURRENT SET OF EXPERIMENTAL (SAMPLE) IMAGES DO IEXP=IEXPT,MIN(NUMEXP,IEXPT+NUMTH-1) IT = IEXP-IEXPT+1 C ILIST(IT) IS LIST NUMBER OF MOST SIMILAR REF. IMAGE IREF = ILIST(IT) IF (IREF .LE. 0) THEN C NO NEARBY REFERENCE IMAGE IMGREF = 0 C IREFT IS FOR REFDIR INDEX IREFT = 1 ELSE IMGREF = IREFLIST(IREF) C IREFT IS FOR REFDIR INDEX IREFT = IREF ENDIF IMGEXP = IEXPLIST(IEXP) CCROT = CCLIST(IT) RANGNEW = RANGNLIST(IT) MIRRORNEW = MIRLIST(IT) PEAKV = 0.0 XSHNEW = 0.0 YSHNEW = 0.0 CALL AP_END(IEXP,IMGEXP,IMGREF, & ANGREF(1,IREFT),REFDIR(1,IREFT), & ANGEXP(1,IEXP), EXPDIR(1,IEXP),ISHRANGE, & GOTREFANG, NGOTPAR,NSAM,NROW,CCROT,PEAKV, & RANGNEW,XSHNEW,YSHNEW,MIRRORNEW,EXPPAT,REFPAT, & NPROJ(IT), CTYPE, XBUF,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) CCROTAVG = CCROTAVG + CCROT IF (CCROT .GE. CCROTLAS) THEN IMPROVCCROT = IMPROVCCROT + 1 CCROTIMPROV = CCROTIMPROV + CCROT ELSE IWORSECCROT = IWORSECCROT + 1 CCROTWORSE = CCROTWORSE + CCROT ENDIF ENDIF 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 IF (ALLOCATED(XBUF)) DEALLOCATE(XBUF) IF (ALLOCATED(ILIST)) DEALLOCATE(ILIST) IF (ALLOCATED(CCLIST)) DEALLOCATE(CCLIST) IF (ALLOCATED(RANGNLIST)) DEALLOCATE(RANGNLIST) IF (ALLOCATED(MIRLIST)) DEALLOCATE(MIRLIST) IF (ALLOCATED(NPROJ)) DEALLOCATE(NPROJ) END C++************************** DSGR_2 ********************************** SUBROUTINE DSGR_2(XIM,NSAM,NROW, NUMR,NRING,MODE, & CIRCREF,LCIRC,NUMREF,TT,LENTT, & SA,TA,RANGE, & IMGREFL,CCROT,RANGNEW,MIRNEW,NPROJ, & CKMIRROR,FFTW_PLANS) C NOTE: RUNS WITHIN OMP PARALLEL SECTION OF CODE! REAL, DIMENSION(NSAM,NROW) :: XIM INTEGER, DIMENSION(3,NRING) :: NUMR CHARACTER(LEN=1) :: MODE REAL, DIMENSION(LCIRC,NUMREF) :: CIRCREF DOUBLE PRECISION, DIMENSION(LENTT) :: TT REAL, DIMENSION(3,NUMREF) :: SA REAL, DIMENSION(3) :: TA REAL :: RANGE LOGICAL :: CKMIRROR,ONLYMIRROR LOGICAL :: MIRNEW,LIMITRANGE LOGICAL :: MIRRORED INTEGER *8 :: FFTW_PLANS(*) DOUBLE PRECISION :: CCROTD,TOTMIN,TOTMIR C AUTOMATIC ARRAYS INTEGER, DIMENSION(NUMREF) :: LCG REAL, DIMENSION(LCIRC) :: CIRCEXP !integer, dimension(lcirc/2) :: iptr MAXRIN = NUMR(3,NRING) #ifdef SP_LIBFFTW3 MAXRIN = NUMR(3,NRING) - 2 #endif WR = 0.0 NPROJ = NUMREF LIMITRANGE = (RANGE .LT. 1.0) IF (LIMITRANGE) THEN C DETERMINE WHICH ONES ARE TO BE COMPARED NPROJ = 0 DO IMI=1,NUMREF C ABS - DIRECTIONS AT 180 DEGREES ARE THE SAME (MIRRORED) C DT NEAR 1.0 = NOT-MIRRORED, DT NEAR -1.0 = MIRRORED DT = (TA(1) * SA(1,IMI) + & TA(2) * SA(2,IMI) + & TA(3) * SA(3,IMI)) DTABS = ABS(DT) IF (DTABS .GE. RANGE) THEN C EITHER MIRRORED OR NON-MIRRORED IS WITHIN RANGE NPROJ = NPROJ + 1 LCG(NPROJ) = IMI IF (DT .LT. 0) LCG(NPROJ) = -IMI ENDIF ENDDO IF (NPROJ .LE. 0) THEN C NO REF. IMAGE WITHIN COMPARISON ANGLE IMGREFL = 0 CCROT = -1.0 RANGNEW = 0.0 MIRNEW = .FALSE. RETURN ENDIF ENDIF C CALCULATE DIMENSIONS FOR NORMALIZING IN APRINGS_ONE CNS2 = NSAM/2+1 CNR2 = NROW/2+1 C EXTRACT EXP. IMAGE RINGS, NORMALIZE & FFT THEM CALL APRINGS_ONE_NEW(NSAM,NROW, CNS2,CNR2, XIM, .FALSE., & MODE,NUMR,NRING,LCIRC, WR,FFTW_PLANS, & CIRCEXP,IRTFLG) CCROTD = -1.0D20 IMGREFL = 0 DO IMIL=1,NPROJ IMI = IMIL IF (LIMITRANGE) IMI = ABS(LCG(IMIL)) IF ((CKMIRROR .AND. LIMITRANGE) .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,TOT, & 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,TOT, TOTMIR,TMT, & TT, FFTW_PLANS(1)) ENDIF IF (TOTMIN .GE. CCROTD) THEN C GOOD MATCH WITH TOTA (MIRRORED OR NOT) POSITION CCROTD = TOTMIN RANGNEW = TOT MIRNEW = (LIMITRANGE .AND. (LCG(IMIL) .LT. 0)) IMGREFL = IMI ENDIF IF (CKMIRROR .AND. .NOT. LIMITRANGE) THEN C HAVE TO COMPARE WITH MIRRORED POSITION IF (TOTMIR .GE. CCROTD .AND. & TOTMIR .GT. TOTMIN) THEN C GOOD MATCH WITH MIRRORED POSITION CCROTD = TOTMIR RANGNEW = TMT MIRNEW = .TRUE. IMGREFL = IMI ENDIF ENDIF C END OF: DO IMIL=1,NPROJ ENDDO IF (MODE .EQ. 'F') THEN RANGNEW = (RANGNEW-1.0) / MAXRIN*360.0 ELSE RANGNEW = (RANGNEW-1.0) / MAXRIN*180.0 ENDIF CCROT = CCROTD END cXXXXXXXXXXXXXXXXXXXXX UNUSED BELOW HERE XXXXXXXXXXXXXXXXXXXXXXXXXXX #ifdef NEVER allocate(iptr(lcirc/2),stat=irtflg) if (irtflg .ne. 0) then call errt(46,'iptr...',lcirc/2) goto 9999 endif c load iptr with location in q do i=1,nring igo = (numr(2,i)/2) + 1 nval = numr(3,i) j = 1 do iloc=igo,igo+(nval/2)-1 iptr(iloc) = j j = j + 1 enddo ccc write(6,*) ' iptr(',igo,'): ',iptr(igo),igo,nval enddo write(6,*) ' calling crosrng_e OK' call crosrng_e(circref(1,imi),circexp,lcirc,nring, & maxrin,numr, totmin,tot, & tt, .false.) call crosrng_new_e(circref(1,imi),circexp,lcirc,nring, & maxrin,numr, totmin,tot, & tt, .false.) igot = 1 do i=1,lcirc,2 circs(igot) = circref(i,imi) circs(igot+1) = circref(i+1,imi) circs(igot+2) = circexp(i) circs(igot+3) = circexp(i+1) igot = igot + 4 enddo call crosrng_new_e2(circs,lcirc,nring, & maxrin,numr, totmin,tot, & tt, .false.) call crosrng_new_c(circref(1,imi),circexp,lcircd2,nring, & maxrin,numr, totmin,tot, & tt, .false.,iptr) #endif