C++********************************************************************* C C CORR1.F REWRITTEN AUG 96 PP C ADDED 'CC H' MAR 02 ARDEAN LEITH C OPFILEC FEB 03 ARDEAN LEITH C 'CC MS' bug OCT 03 ARDEAN LEITH C FMRS USED, UNUSED ALLOCS REMOVED FEB 08 ARDEAN LEITH C PHASE REMOVED FEB 08 ARDEAN LEITH C MOD PGI COMPILER BUG FEB 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 CORR1 C C 'CC N' CALL TREE: C CORR1 ---> READV C ---> NRMS C ---> FMRS C C ---> READV C ---> NRMS C ---> FMRS C C ---> CCRS C C23456789012345678901234567890123456789012345678901234567890123456789012 C--********************************************************************* SUBROUTINE CORR1 #ifdef SP_LIBFFTW3 USE TYPE_KINDS INTEGER(KIND=I_8) :: IPLAN = 0 !STRUCTURE POINTER #endif INCLUDE 'CMBLOCK.INC' INCLUDE 'CMLIMIT.INC' CHARACTER(LEN=MAXNAM) :: FILNAM1,FILNAM2,FILNAMM COMPLEX, ALLOCATABLE, DIMENSION(:,:,:) :: QK1 COMPLEX, ALLOCATABLE, DIMENSION(:,:,:) :: QK2 COMPLEX, ALLOCATABLE, DIMENSION(:) :: QKB CHARACTER(LEN=1) :: NULL PARAMETER (NFUNC=3) CHARACTER(LEN=2),DIMENSION(NFUNC) :: FUNC LOGICAL :: ACASE, BOTH_INC LOGICAL :: SPIDER_SIGN LOGICAL :: SPIDER_SCALE DATA FUNC/'AC', 'CC', 'CN'/ DATA LUN1,LUN2,LUN3/21,22,23/ NULL = CHAR(0) IRTFLG = 0 C DETERMINE IFUNC DO IFUNC = 1,NFUNC IF (FCHAR(1:2) .EQ. FUNC(IFUNC)(1:2)) THEN GOTO 1111 ENDIF ENDDO C OPERATION NOT HERE, RETURN TO CALLER RETURN 1111 CONTINUE C CATCH EXCEPTIONS C ---------------------------------------------------------- CC C IF (FCHAR(1:4) .EQ. 'CC C') THEN MAXIM = 0 CALL OPFILEC(0,.TRUE.,FILNAM1,LUN1,'O',IFORM, & NSAM1,NROW1,NSLICE1,MAXIM,'INPUT1',.TRUE.,IRTFLG) IF (IRTFLG .NE. 0) RETURN MAXIM = 0 CALL OPFILEC(0,.TRUE.,FILNAM2,LUN2,'O',IFORM, & NSAM2,NROW2,NSLICE2,MAXIM,'INPUT1',.TRUE.,IRTFLG) IF (IRTFLG .NE. 0) THEN GOTO 998 ELSEIF (NSAM1 .NE. NSAM2 .OR. & NROW1 .NE. NROW2 .OR. & NSLICE1 .NE. NSLICE2) THEN CALL ERRT(1,'CORR1',NE) GOTO 998 ENDIF MAXIM = 0 CALL OPFILEC(0,.TRUE.,FILNAMM,LUN3,'O',IFORM, & NSAMM,NROWM,NSLICEM,MAXIM,'MASK',.FALSE.,IRTFLG) IF (IRTFLG .NE. 0) GOTO 998 IF (NSAM1 .NE. NSAMM .OR. & NROW1 .NE. NROWM .OR. & NSLICE1 .NE. NSLICEM) THEN CALL ERRT(1,'CCC',NE) GOTO 998 ENDIF CALL CCC(LUN1,FILNAM1,NSAM1,NROW1,NSLICE1, & LUN2,FILNAM2,NSAM2,NROW2,NSLICE2, & LUN3,FILNAMM) GOTO 998 C ----------------------------------------------------- CC MS ELSEIF(FCHAR(1:5) .EQ. 'CC MS') THEN C CROSS CORRELATION - MASKED AND NORMALIZED CALL MCCF RETURN C ------------------------------------------------------ AC MS ELSEIF (FCHAR(1:5) .EQ. 'AC MS') THEN IF (FCHAR(6:6) .EQ. 'S') THEN C AUTO CORRELATION - MASKED AND NORMALIZED CALL MACF('S') ELSE C AUTO CORRELATION - MASKED AND NORMALIZED CALL MACF(' ') ENDIF RETURN C ------------------------------------------------------ CC P ELSEIF(FCHAR(1:4).EQ.'CC P') THEN C THIS APPEARS TO BE UNUSED ELSEWHERE IN SPIDER?? CALL POLAR_CC RETURN C ------------------------------------------------------ CC H ELSEIF (FCHAR(4:4) .EQ. 'H') THEN CALL ERRT(46,'OPERATION -CC H- NO LONGER SUPPORTED',IER) RETURN ENDIF C -------------------------------------------------------- OTHER C OPEN FIRST INPUT FILE, FOURIER INPUT ALLOWED MAXIM = 0 CALL OPFILEC(0,.TRUE.,FILNAM1,LUN1,'O',IFORM, & NSAM1,NROW1,NSLICE1,MAXIM,'INPUT',.TRUE.,IRTFLG) IF (IRTFLG .NE. 0) RETURN IF (FCHAR(4:4) .EQ. 'N' .AND. IFORM .LT. 0) THEN C NO STATISTICS IN FILE, CAN NOT NORMALIZE OUTPUT CALL ERRT(101, & 'CAN NOT NORMALIZE OUTPUT - FILE LACKS STATISTICS.',NE) GOTO 998 ELSEIF (FCHAR(4:4) .EQ. 'N') THEN C FOR 'CC N' IMAMI1 = IMAMI SIG1 = SIG ENDIF IFORM1 = IFORM C CALCULATE DIMENSIONS IF (IFORM1 .GT. 0) THEN C REAL SPACE INPUT LS1 = NSAM1+2-MOD(NSAM1,2) LREC1 = NSAM1 IFORM3 = IFORM1 NSAM3 = NSAM1 ELSE C FOURIER SPACE INPUT LS1 = NSAM1 LREC1 = NSAM1 NSAM1 = NSAM1-MOD(-IFORM1,10) IF (IFORM1 .GT. -20) THEN IFORM3 = 1 ELSE IFORM3 = 3 ENDIF NSAM3 = NSAM1 ENDIF IF (FCHAR(1:2) .EQ. 'AC') THEN C AUTO CORRELATION WANTED, NO SECOND FILE ACASE = .TRUE. IFORM2 = IFORM1 ELSE C GET NAME FOR SECOND INPUT FILE CALL FILERD(FILNAM2,NLET,NULL,'REFERENCE',IRTFLG) IF (IRTFLG .NE. 0) GOTO 998 IF (FILNAM1 .EQ. FILNAM2) THEN C FILENAMES ARE SAME, AUTOCORRELATION WANTED ACASE = .TRUE. IFORM2 = IFORM1 ELSE C OPEN SECOND INPUT FILE, FOURIER INPUT ALLOWED MAXIM = 0 CALL OPFILEC(0,.FALSE.,FILNAM2,LUN2,'O',IFORM, & NSAM2,NROW2,NSLICE2,MAXIM,NULL,.TRUE.,IRTFLG) IF (IRTFLG .NE. 0) GOTO 998 IFORM2 = IFORM ACASE = .FALSE. IF (FCHAR(4:4) .EQ. 'N' .AND. IFORM .LT. 0) THEN C NO STATISTICS IN FILE, CAN NOT NORMALIZE OUTPUT CALL ERRT(101, & 'CAN NOT NORMALIZE OUTPUT - FILE LACKS STATISTICS.',I) GOTO 998 ELSEIF(FCHAR(4:4) .EQ. 'N') THEN C FOR 'CC N' IMAMI2 = IMAMI SIG2 = SIG ENDIF C CALCULATE DIMENSIONS IF (IFORM2 .GT. 0) THEN C REAL SPACE INPUT LS2 = NSAM2+2-MOD(NSAM2,2) ELSE C FOURIER SPACE INPUT LS2 = NSAM2 NSAM2 = NSAM2-MOD(-IFORM2,10) ENDIF C CHECK THAT DIMENSIONS ARE THE SAME FOR BOTH FILES IF (NSAM1 .NE. NSAM2 .OR. NROW1 .NE. NROW2 .OR. & NSLICE1 .NE. NSLICE2) THEN C ERROR. IMAGES DO NOT HAVE SAME DIMENSIONS CALL ERRT(1,'CORR1',NE) GOTO 998 ENDIF ENDIF ENDIF C -------------------------------------------------- FIRST INPUT ALLOCATE(QK1(LS1/2,NROW1,NSLICE1), STAT=IRTFLG) IF (IRTFLG .NE. 0) THEN MWANT = LS1/2*NROW1*NSLICE1 CALL ERRT(46,'CORR1, QK1',MWANT) GOTO 996 ENDIF CALL READV(LUN1,QK1,LS1,NROW1,LREC1,NROW1,NSLICE1) IF (FCHAR(4:4) .EQ. 'N' .AND. IMAMI1 .NE. 1) & CALL NRMS(QK1,LS1,NSAM1,NROW1,NSLICE1,SIG1) IF (IFORM1 .GT. 0) THEN C REAL SPACE INPUT, TRANSFORM IT TO FOURIER INV = +1 SPIDER_SIGN = .TRUE. SPIDER_SCALE = .TRUE. CALL FMRS(QK1, NSAM1,NROW1,NSLICE1,IPLAN, & SPIDER_SIGN,SPIDER_SCALE, INV,IRTFLG) IF (IRTFLG .NE. 0) THEN CALL ERRT(101,FOURIER TRANSFORM FAILED,NE) GOTO 996 ENDIF ENDIF IF (FCHAR(4:4) .EQ. 'N') THEN QK1(1,1,1) = (0.0,0.0) ENDIF C ------------------------------------------------- SECOND INPUT IF (.NOT. ACASE .AND. IFORM2 .GT. 0) THEN C CROSS-CORRELATION WITH REAL IMAGES, USE IN-CORE BOTH_INC = .TRUE. ALLOCATE(QK2(LS2/2,NROW2,NSLICE2), STAT=IRTFLG) IF (IRTFLG .NE. 0) THEN MWANT = LS2 / 2 * NROW2 * NSLICE2 CALL ERRT(46,'CORR1, QK2',MWANT) GOTO 998 ENDIF CALL READV(LUN2,QK2,LS2,NROW2,NSAM2,NROW2,NSLICE2) IF (FCHAR(4:4) .EQ. 'N' .AND. IMAMI2 .NE. 1) & CALL NRMS(QK2,LS2,NSAM2,NROW2,NSLICE2,SIG2) C REAL SPACE INPUT, TRANSFORM IT TO FOURIER INV = +1 SPIDER_SIGN = .TRUE. SPIDER_SCALE = .TRUE. CALL FMRS(QK2, NSAM2,NROW2,NSLICE2, IPLAN, & SPIDER_SIGN,SPIDER_SCALE, INV,IRTFLG) IF (IRTFLG .NE. 0) THEN CALL ERRT(101,FOURIER TRANSFORM FAILED,NE) GOTO 996 ENDIF ELSE BOTH_INC = .FALSE. ALLOCATE (QKB((NSAM1+2-MOD(NSAM1,2))/2), STAT=IRTFLG) IF (IRTFLG .NE. 0) THEN MWANT = NSAM1+2-MOD(NSAM1,2)/2 CALL ERRT(46,'CORR1, QKB',MWANT) GOTO 996 ENDIF ENDIF C -------------------------------------------------- CORRELATION LS = NSAM1+2-MOD(NSAM1,2) IF (ACASE) THEN C 'AC' AUTO CORRELATION IF (NSLICE1.LE.1) THEN IF (FCHAR(4:4) .EQ. 'S' .OR. FCHAR(4:5).EQ. 'NS') THEN CALL ACRS_2S(QK1,QK1, LS,NSAM1,NROW1) ELSE CALL ACRS_2(QK1,QK1, LS,NSAM1,NROW1) ENDIF ELSE IF (FCHAR(4:4) .EQ. 'S' .OR. FCHAR(4:5) .EQ. 'NS') THEN CALL ACRS_3S(QK1,QK1, LS,NSAM1,NROW1,NSLICE1) ELSE CALL ACRS_3(QK1,QK1, LS,NSAM1,NROW1,NSLICE1) ENDIF ENDIF ELSEIF (FCHAR(1:2) .EQ. 'CC') THEN C 'CC' CROSS CORRELATION IF (BOTH_INC) THEN SPIDER_SIGN = .TRUE. SPIDER_SCALE = .TRUE. CALL CCRS(QK1,QK2,QK1, LS,NSAM1,NROW1,NSLICE1, & SPIDER_SIGN,SPIDER_SCALE, IRTFLG) IF (IRTFLG .NE. 0) THEN CALL ERRT(101,'CCRS FAILED',NE) GOTO 996 ENDIF ELSEIF (NSLICE1 .LE. 1) THEN CALL CCRD_2(LUN2,QK1,QKB,QK1, LS,NSAM1,NROW1) ELSE CALL CCRD_3(LUN2,QK1,QKB,QK1, LS,NSAM1,NROW1,NSLICE1) ENDIF ELSEIF (FCHAR(1:2) .EQ. 'CN') THEN C 'CN' CONVOLUTION, NOT CORRELATION IF (NSLICE1.LE.1) THEN IF (BOTH_INC) THEN CALL CNRS_2(QK1,QK2,QK1, LS,NSAM1,NROW1) ELSE CALL CNRD_2(LUN2,QK1,QKB,QK1, LS,NSAM1,NROW1) ENDIF ELSE IF (BOTH_INC) THEN CALL CNRS_3(QK1,QK2,QK1, LS,NSAM1,NROW1,NSLICE1) ELSE CALL CNRD_3(LUN2,QK1,QKB,QK1, LS,NSAM1,NROW1,NSLICE1) ENDIF ENDIF ENDIF C ------------------------------------------------------ OUTPUT IF (ACASE) SIG2 = SIG1 IF (FCHAR(4:4) .EQ. 'N') THEN C NORMALIZATION ROUTINE HERE FAN = 1.0 / (NSAM1*NROW1*FLOAT(NSLICE1)-1.0) / SIG1 / SIG2 QK1 = QK1 * FAN ENDIF C SAVE RESULTS IN OUTPUT FILE IFORM = IFORM3 MAXIM = 0 CALL OPFILEC(LUN1,.TRUE.,FILNAM1,LUN3,'U',IFORM3, & NSAM3,NROW1,NSLICE1,MAXIM,'OUTPUT',.TRUE.,IRTFLG) IF (IRTFLG .NE. 0) GOTO 996 C THIS ONLY WRITES FIRST LS1 VALUES FROM EACH ROW CALL WRITEV(LUN3,QK1,LS1,NROW1,NSAM1,NROW1,NSLICE1) 996 IF (ALLOCATED(QKB)) DEALLOCATE (QKB) IF (ALLOCATED(QK1)) DEALLOCATE(QK1) IF (ALLOCATED(QK2)) DEALLOCATE(QK2) 998 CLOSE(LUN3) CLOSE(LUN2) CLOSE(LUN1) RETURN END