C++********************************************************************* C C APSCC NEW FEB 08 ARDEAN LEITH C FLIPPED EXP & REF BUG AUG 09 ARDEAN LEITH C REMOVED FROM APMASTER 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 APSCC_DO(IREFLST,NUMREF, IEXPLST,NUMEXP, C NSAM,NROW,NSLICE, ISHX,ISHY,ISHZ, NORMIT, C REFPAT,EXPPAT, LUNREF,LUNEXP,LUNDOC) C C PURPOSE: 2D & 3D PADDED, CROSS CORRELATION MULTI-REFERENCE C SHIFT ALIGNMENT C C PARAMETERS: C IREFLST LIST OF REF. IMAGE FILE NUMBERS (INPUT) C NUMREF NO. OF IMAGES (INPUT) C IEXPLST LIST OF EXP. IMAGE FILE NUMBERS (INPUT) C NUMEXP NO. OF IMAGES (INPUT) C NSAM,NROW,NSLICE ACTUAL (NOT-WINDOWED) IMAGE SIZE (INPUT) C ISHX,ISHY,ISHZ SHIFT SEARCH RANGE (INPUT) C NORMIT NORMALIZATION WANTED FLAG (INPUT) C REFPAT REF. IMAGE SERIES FILE TEMPLATE (INPUT) C EXPPAT EXP. IMAGE SERIES FILE TEMPLATE (INPUT) C LUNREF,LUNEXP,LUNDOC I/O UNITS (INPUT) C C23456789012345678901234567890123456789012345678901234567890123456789012 C--********************************************************************* SUBROUTINE APSCC() INCLUDE 'CMLIMIT.INC' INCLUDE 'CMBLOCK.INC' INTEGER, ALLOCATABLE, DIMENSION(:) :: IMGLST INTEGER, ALLOCATABLE, DIMENSION(:,:) :: NUMR CHARACTER (LEN=MAXNAM) :: FILNAM,REFPAT,EXPPAT,OUTANG CHARACTER(LEN=1) :: MODE,NULL,YN CHARACTER(LEN=160) :: COMMEN LOGICAL :: NEWFILE,NORMIT CALL SET_MPI(ICOMM,MYPID,MPIERR) ! SETS ICOMM AND MYPID DATA LUNREF,LUNEXP/50,51/ DATA INPIC,INANG,NDOC,NSCF/77,78,55,50/ !USED IN CALLED ROUTINE NULL = CHAR(0) C ASK FOR TEMPLATE AND NUMBERS FOR REFERENCE IMAGES NILMAX = NIMAX CALL FILELIST(.TRUE.,LUNREF,REFPAT,NLET,INUMBR,NILMAX,NUMREF, & 'ENTER TEMPLATE FOR REFERENCE IMAGES',IRTFLG) IF (IRTFLG .NE. 0) RETURN C NUMREF - TOTAL NUMBER OF REF. IMAGES IF (NUMREF .LE. 0) THEN CALL ERRT(101,'No reference images',IDUM) GOTO 9999 ENDIF IF (MYPID .LE. 0) WRITE(NOUT,90) NUMREF 90 FORMAT(' Number of reference images: ',I7) C GET FIRST REFERENCE IMAGE TO DETERMINE DIMENSIONS NLET = 0 CALL FILGET(REFPAT,FILNAM,NLET,INUMBR(1),INTFLG) MAXIM = 0 CALL OPFILEC(0,.FALSE.,FILNAM,LUNREF,'O',IFORM,NSAM,NROW,NSLICE, & MAXIM,' ',.FALSE.,IRTFLG) IF (IRTFLG .NE. 0) GOTO 9999 CLOSE(LUNREF) C MULTI-SHIFT CROSS-CORRELATION ISHX = 0 ISHY = 0 ISHZ = 0 CALL RDPRI3S(ISHX,ISHY,ISHZ,NOT_USED, & 'TRANSLATION SEARCH RANGE IN X,Y,& Z (ZERO FOR ALL)' ,IRTFLG) IF (IRTFLG .NE. 0) GOTO 9999 C POSSIBLE SHIFT AMOUNTS FOR SEARCH IF (ISHX .EQ. 0) ISHX = NSAM / 2 IF (ISHY .EQ. 0) ISHY = NROW / 2 IF (ISHZ .EQ. 0) ISHZ = NSLICE / 2 CALL RDPRMC(YN,NLET,.TRUE.,'NORMALIZE PEAK HEIGHT (Y/N)', & NULL,IRTFLG) IF (IRTFLG .NE. 0) GOTO 9999 NORMIT = (YN .EQ. 'Y') C GET LIST OF EXPERIMENTAL IMAGES TO BE ALIGNED ALLOCATE(IMGLST(NILMAX),STAT=IRTFLG) IF (IRTFLG .NE. 0) THEN CALL ERRT(46,'IMGLST',NILMAX) GOTO 9999 ENDIF CALL FILELIST(.TRUE.,LUNEXP,EXPPAT,NLEP, & IMGLST,NILMAX,NUMEXP, & 'ENTER TEMPLATE FOR IMAGE SERIES TO BE ALIGNED',IRTFLG) IF (IRTFLG .NE. 0) GOTO 9999 IF (MYPID .LE. 0) WRITE(NOUT,91) NUMEXP 91 FORMAT(' Number of experimental images: ',I6) C GET NAME FOR OUTPUT DOC FILE CALL REG_GET_USED(NSEL_USED) C OPEN OUTPUT DOC FILE (FOR APPENDING) NOUTANG = NDOC CALL OPENDOC(OUTANG,.TRUE.,NLET,NDOC,NOUTANG,.TRUE., & 'OUTPUT ALIGNMENT DOCUMENT',.FALSE.,.TRUE.,.TRUE., & NEWFILE,IRTFLG) IF (IRTFLG .NE. 0) GOTO 9999 IF (IRTFLG .EQ. -1) THEN C DO NOT WANT OUTPUT DOC FILE NOUTANG = 0 ELSE C WANT OUTPUT DOC FILE COMMEN = ' ' // & 'EXP#, REF#, SX, SY, '// & 'SZ, PEAK' CALL LUNDOCPUTCOM(NOUTANG,COMMEN,IRTFLG) ENDIF NORMIT = .TRUE. CALL APSCC_DO(INUMBR,NUMREF ,IMGLST,NUMEXP, & NSAM,NROW,NSLICE, ISHX,ISHY,ISHZ, NORMIT, & REFPAT,EXPPAT,LUNREF,LUNEXP,NOUTANG) 9999 IF (ALLOCATED(IMGLST)) DEALLOCATE(IMGLST) IF (ALLOCATED(NUMR)) DEALLOCATE(NUMR) CLOSE(NDOC) END C ------------------- APSCC_DO ---------------------------------- SUBROUTINE APSCC_DO(IREFLST,NUMREF, IEXPLST,NUMEXP, & NSAM,NROW,NSLICE, ISHX,ISHY,ISHZ, NORMIT, & REFPAT,EXPPAT, LUNREF,LUNEXP,LUNDOC) INCLUDE 'CMLIMIT.INC' INCLUDE 'CMBLOCK.INC' INTEGER, DIMENSION(NUMREF) :: IREFLST INTEGER, DIMENSION(NUMEXP) :: IEXPLST CHARACTER (LEN=*) :: REFPAT,EXPPAT LOGICAL :: NORMIT,SPIDER_SIGN CHARACTER (LEN=MAXNAM) :: FILNAM C ALLOCATABLE ARRAYS REAL, ALLOCATABLE, DIMENSION(:) :: BUFI,BUFR,BUFE C AUTOMATIC ARRAYS REAL :: DLIST(6) CALL REG_GET_USED(NSEL_USED) NSAMP = NSAM * 2 NROWP = NROW * 2 NSLICEP = NSLICE * 2 NSAMP1 = NSAM + 1 NROWP1 = NROW + 1 NSLICEP1 = NSLICE + 1 C EXTRA SPACE FOR FOURIER TRANSFORM LSE = NSAMP + 2 - MOD(NSAMP,2) C MAKE BUFI & BUFR FOR PADDED IMAGES WITH FOURIER ROW LENGTH ALLOCATE(BUFR(LSE*NROWP*NSLICEP), & BUFI(LSE*NROWP*NSLICEP), & STAT=IRTFLG) IF (IRTFLG .NE. 0) THEN MWANT = 2 * LSE * NROWP * NSLICEP CALL ERRT(46,'BUFI & BUFR...',MWANT) GOTO 9999 ENDIF IF (NUMREF .GT. 1) THEN C MUST SAVE BUFI SINCE IT IS OVERWRITTEN BY APCC ALLOCATE(BUFE(LSE*NROWP*NSLICEP), & STAT=IRTFLG) IF (IRTFLG .NE. 0) THEN MWANT = 1 * LSE * NROWP * NSLICEP CALL ERRT(46,'BUFE...',MWANT) GOTO 9999 ENDIF ENDIF DO IEXP = 1,NUMEXP NLET = 0 CALL FILGET(EXPPAT,FILNAM,NLET,IEXPLST(IEXP),IRTFLG) IF (IRTFLG .NE. 0) GOTO 9999 MAXIM = 0 CALL OPFILEC(0,.FALSE.,FILNAM,LUNEXP,'O',IFORM, & LSAM,LROW,LSLICE,MAXIM,' ',.FALSE.,IRTFLG) IF (IRTFLG .NE. 0) GOTO 9999 IF (LSAM.NE.NSAM.OR.LROW.NE.NROW.OR.LSLICE.NE.NSLICE)THEN CALL ERRT(102,'INCONSISTENT IMAGE SIZE',LSAM) IRTFLG = 1 GOTO 9999 ENDIF IF (IMAMI .NE. 1) & CALL NORM3(LUNEXP,NSAM,NROW,NSLICE,FMAX,FMIN,AV) AVI = AV ! NOT IMPLEMENTED AVI = 0.0 C LOAD & PAD EXP. IMAGE TO DOUBLE SIZE ILOC = 1 IREC = 1 DO ISLICE = 1,NSLICE DO IROW = 1,NROW CALL REDLIN(LUNEXP,BUFI(ILOC),NSAM,IREC) IREC = IREC + 1 C FILL REMAINING PADDING COLS BUFI(ILOC+NSAM: ILOC+LSE-1) = AVI BUFI(ILOC+LSE*NROW:ILOC+LSE*NROW+LSE-1) = AVI ILOC = ILOC + LSE ENDDO ILOC = ILOC + LSE * NROW ENDDO C FILL REMAINING PADDING SLICES BUFI(LSE*NROW*NSLICE+1:LSE*NROWP*NSLICEP) = AVI CLOSE(LUNEXP) IF (NORMIT) CALL NRMS(BUFI,LSE,NSAMP,NROWP,NSLICEP,SIGI) !write(6,*)'sigif,sigi: ',sigif,sigi IF (NUMREF .GT. 1) THEN C MUST SAVE BUFI SINCE IT IS OVERWRITTEN BY APCC BUFE = BUFI ENDIF DO IREF = 1,NUMREF NLET = 0 CALL FILGET(REFPAT,FILNAM,NLET,IREFLST(IREF),IRTFLG) IF (IRTFLG .NE. 0) RETURN MAXIM = 0 CALL OPFILEC(0,.FALSE.,FILNAM,LUNREF,'O',IFORM, & LSAM,LROW,LSLICE,MAXIM,' ',.FALSE.,IRTFLG) IF (IRTFLG .NE. 0) RETURN IF (LSAM.NE.NSAM.OR.LROW.NE.NROW.OR.LSLICE.NE.NSLICE)THEN CALL ERRT(102,'INCONSISTENT IMAGE SIZE',LSAMT) IRTFLG = 1 RETURN ENDIF IF (IMAMI .NE. 1) & CALL NORM3(LUNREF,NSAM,NROW,NSLICE,FMAX,FMIN,AV) AVR = AV ! NOT IMPLEMENTED AVR = 0.0 C LOAD & PAD REFERENCE IMAGE TO DOUBLE SIZE ILOC = 1 IREC = 1 DO ISLICE = 1,NSLICE DO IROW = 1,NROW CALL REDLIN(LUNREF,BUFR(ILOC),NSAM,IREC) IREC = IREC + 1 C FILL REMAINING PADDING COLS BUFR(ILOC+NSAM:ILOC+LSE-1) = AVR BUFR(ILOC+LSE*NROW:ILOC+LSE*NROW+LSE-1) = AVR ILOC = ILOC + LSE ENDDO ILOC = ILOC + LSE * NROW ENDDO C FILL REMAINING PADDING SLICES BUFR(LSE*NROW*NSLICE+1:LSE*NROWP*NSLICEP) = AVR CLOSE(LUNREF) IF (NORMIT) CALL NRMS(BUFR,LSE,NSAMP,NROWP,NSLICEP,SIGR) ! write(6,*)'sigrf,sigr: ',sigrf,sigr C CROSS CORRELATION --------------------------------------- CC C APCC RETURNS PEAK IMAGE IN: BUFR, BUFI(0,0) SET TO 0,0 SPIDER_SIGN = .FALSE. CALL APCC(LSE, NSAMP,NROWP,NSLICEP, BUFI,BUFR, & .FALSE.,NORMIT,SPIDER_SIGN, ISHX,ISHY,ISHZ, & XSHNEW,YSHNEW,ZSHNEW, PEAKV,IRTFLG) IF (IRTFLG .NE. 0) THEN CALL ERRT(101,' CC ERROR',NE) GOTO 9999 ENDIF !write(6,*)'peakv,sigi,sigr: ',peakv,sigi,sigr IF (NORMIT) THEN C NORMALIZATION PEAKV = PEAKV/FLOAT(NSAMP*NROWP*NSLICEP-1)/SIGI/SIGR ENDIF DLIST(1) = IEXPLST(IEXP) DLIST(2) = IREFLST(IREF) DLIST(3) = XSHNEW DLIST(4) = YSHNEW DLIST(5) = ZSHNEW DLIST(6) = PEAKV IF (LUNDOC .GT. 0) THEN C SAVE IN ALIGNMENT DOC FILE C IMG#, REF#, SX,SY,SZ, PEAK-HEIGHT CALL LUNDOCWRTDAT(LUNDOC,IEXP,DLIST,6,IRTFLG) IF (IRTFLG .NE. 0) GOTO 9999 ENDIF IF (NSEL_USED .GT. 0) THEN C OUTPUT TO SPIDER'S REGISTERS CALL REG_SET_NSEL(1,5, DLIST(1),DLIST(2),DLIST(3), & DLIST(4),DLIST(5),IRTFLG) CALL REG_SET_NSEL(6,1, DLIST(6),0.0,0.0, & 0.0,0.0,IRTFLG) ENDIF IF (NUMREF .GT. 1) THEN C MUST RETRIEVE BUFI SINCE IT IS OVERWRITTEN BY APCC BUFI = BUFE ENDIF ENDDO ENDDO 9999 IF (ALLOCATED(BUFI)) DEALLOCATE(BUFI) IF (ALLOCATED(BUFR)) DEALLOCATE(BUFR) IF (ALLOCATED(BUFE)) DEALLOCATE(BUFE) END C ----------------- UNUSED BELOW HERE -------------------------- #ifdef NEVER #ifdef DEBUG write(6,*)' ' write(6,*)'peakv: ',peakv !PEAKV = PEAKV / FLOAT(NSAMP*NROWP*NSLICEP) !write(6,*)' ' !write(6,*)'PEAKV /(NSAMP*NROWP*NSLICEP): ',peakv write(6,*)'peakv/sigif/sigrf: ',peakv/sigif/sigrf write(6,*)'peakv/sigi/sigr: ',peakv/sigi/sigr nnn = nsam*nrow*nslice nnnp = nsamp*nrowp*nslicep write(6,*)'peakv/float(nnn)/sigif/sigrf: ', & peakv/float(nnn)/sigif/sigrf write(6,*)'peakv/float(nnnp)/sigif/sigrf: ', & peakv/float(nnnp)/sigif/sigrf write(6,*)'peakv/float(nnn)/sigi/sigr: ', & peakv/float(nnn)/sigi/sigr write(6,*)'peakv/float(nnnp)/sigi/sigr: ', & peakv/float(nnnp)/sigi/sigr write(6,*)' ' write(6,*)'final peakv: ',peakv #endif c---------------------- maxim = 0 itype = 3 call opfilec(0,.false.,'jnkexppad',98,'U',itype, & lse,nrowp,nslicep,maxim,' ',.false.,irtflg) call wrtvol(98,lse,nrowp, 1,nslicep, bufi,irtflg) close(98) maxim = 0 itype = 3 call opfilec(0,.false.,'jnkrefpad',98,'U',itype, & lse,nrowp,nslicep,maxim,' ',.false.,irtflg) call wrtvol(98,lse,nrowp,1,nslicep, bufr,irtflg) close(98) c----------------------- #endif