C++********************************************************************* C C APCC.F NEW FROM APSHIFT FOR SPEEDUP FEB 08 ARDEAN LEITH C CCRS INOUT APR 09 ARDEAN LEITH C ********************************************************************** C=* FROM: SPIDER - MODULAR IMAGE PROCESSING SYSTEM. AUTHOR: J.FRANK * C=* Copyright (C) 1985-2009 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 APCC(LSC,NSAM,NROW,NSLICE, BUFI,BUFR, C SKIP_PEAK,NORMIT,SPIDER_SIGN, C ISHRANGEX,ISHRANGEY,ISHRANGEZ, C XSHNEW,YSHNEW,ZSHNEW,PEAKV,IRTFLG) C C PURPOSE: CROSS-CORRELATE WITH REFERENCE IMAGE, OPTIONALLY FIND C CC PEAK LOCATION. C C PARAMETERS: NSAM... PADDED INPUT IMAGE DIMENSIONS SENT C NSAM... PADDED INPUT IMAGE DIMENSIONS SENT C BUFI SAMPLE IMAGE --> PEAK IMAGE SENT/RET. C BUFR REFERENCE IMAGE -- ALTERED SENT/RET. C SKIP_PEAK SKIP PEAK LOCATION FLAG SENT C NORMIT NORM LOCATION 0 FLAG SENT C SPIDER_SIGN USE SPIDERS FFT SIGN CONVENTION SENT C ISHRANGEX.. POSSIBLE IMAGE SHIFT SENT C XSHNEW.. SUB PIXEL IMAGE SHIFT RET. C PEAKV PEAK HEIGHT RET. C IRTFLG ERROR FLAG RET. C C NOTE: BUFI and BUFR SHOULD CONTAIN IMAGES INSIDE A ROW C LENGTH: NSAM+2-MOD(NSAM,2)). FOR BEST ACCURACY IMAGES C SHOULD BE PADDED TO 2X ORIGINAL SIZE, THEN ROW LENGTH WILL C BE: 2*NSAM+2 (WHERE NSAM IS ORIGINAL IMAGE ROW LENGTH) C C--********************************************************************* SUBROUTINE APCC(LSC, NSAM,NROW,NSLICE, BUFI,BUFR, & SKIP_PEAK,NORMIT,SPIDER_SIGN, & ISHRANGEX,ISHRANGEY,ISHRANGEZ, & XSHNEW,YSHNEW,ZSHNEW,PEAKV,IRTFLG) USE TYPE_KINDS INTEGER(KIND=I_8) :: IPLAN = 0 !STRUCTURE POINTER INCLUDE 'CMBLOCK.INC' REAL :: BUFI(LSC,NROW,NSLICE) REAL :: BUFR(LSC,NROW,NSLICE) INTEGER, INTENT(IN) :: NSAM,NROW,NSLICE LOGICAL, INTENT(IN) :: SKIP_PEAK,NORMIT,SPIDER_SIGN INTEGER, INTENT(IN) :: ISHRANGEX,ISHRANGEY,ISHRANGEZ REAL, INTENT(OUT) :: XSHNEW,YSHNEW,ZSHNEW,PEAKV INTEGER, INTENT(OUT) :: IRTFLG INTEGER, DIMENSION(3) :: ILOCS LOGICAL :: SPIDER_SCALE IRTFLG = 0 SPIDER_SCALE = .FALSE. C CROSS CORRELATION --------------------------------------- CC INV = +1 CALL FMRS(BUFI, NSAM,NROW,NSLICE, IPLAN, & SPIDER_SIGN,SPIDER_SCALE, INV,IRTFLG) IF (IRTFLG .NE. 0) THEN CALL ERRT(101,'APCC FFT ERROR ON BUFI',NE) RETURN ENDIF INV = +1 CALL FMRS(BUFR, NSAM,NROW,NSLICE, IPLAN, & SPIDER_SIGN, SPIDER_SCALE, INV,IRTFLG) IF (IRTFLG .NE. 0) THEN CALL ERRT(101,'APCC FFT ERROR ON BUFR',NE) RETURN ENDIF IF (NORMIT) THEN ! WARNING: BUFI IS RETURNED!! BUFI(1,1,1) = 0.0 BUFI(2,1,1) = 0.0 ENDIF CALL CCRS(BUFI,BUFR, LSC,NSAM,NROW,NSLICE, & SPIDER_SIGN,SPIDER_SCALE, IRTFLG) IF (IRTFLG .NE. 0) THEN CALL ERRT(101,'CCRS CROSS CORRELATION ERROR',NE) RETURN ENDIF IF (SKIP_PEAK) RETURN ! DO NOT WANT TO FIND PEAK LOCATION C PEAK SEARCH ON THE WINDOWED CC IMAGE -------------------- PEAK C WINDOW SIZE & WINDOW CORNER IGOX = NSAM/2 - ISHRANGEX +1 IENDX = IGOX + 2 * ISHRANGEX IGOY = NROW/2 - ISHRANGEY +1 IENDY = IGOY + 2 * ISHRANGEY IGOZ = NSLICE/2 - ISHRANGEZ +1 IENDZ = IGOZ + 2 * ISHRANGEZ C PEAK SEARCH WITHIN SEARCH RANGE WINDOW ILOCS = MAXLOC(BUFI(IGOX:IENDX, IGOY:IENDY, IGOZ:IENDZ)) IX = IGOX + ILOCS(1) - 1 IY = IGOY + ILOCS(2) - 1 IZ = IGOZ + ILOCS(3) - 1 PEAKV = BUFI(IX,IY,IZ) C FFT FASTER IF SCALING SKIPPED THERE, SO MUST REVALUE OUTPUT PEAKV = PEAKV / FLOAT(NSAM * NROW * NSLICE) IXSH = -(ILOCS(1) - ISHRANGEX - 1) IYSH = -(ILOCS(2) - ISHRANGEY - 1) IZSH = -(ILOCS(3) - ISHRANGEZ - 1) XSHNEW = IXSH ! SET HERE IN CASE OF BOUNDARY RETURN YSHNEW = IYSH ZSHNEW = IZSH IF (NSLICE .GT. 1) THEN !-------------------------- VOLUME IF (IX .LE. 1 .OR. IX .GE. NSAM .OR. & IY .LE. 1 .OR. IY .GE. NROW .OR. & IZ .LE. 1 .OR. IZ .GE. NSLICE) RETURN C BINARY INTERPOLATION TO FIND SUB-PIXEL PEAK QMAX = BUFI(IX,IY,IZ) ! CENTRAL ELEMENT C X, Y & Z SUB-PIXEL SHIFT CALL PKSR3_SUB(QMAX,BUFI(IX-1,IY,IZ),BUFI(IX+1,IY,IZ),1,RX) CALL PKSR3_SUB(QMAX,BUFI(IX,IY-1,IZ),BUFI(IX,IY+1,IZ),1,RY) CALL PKSR3_SUB(QMAX,BUFI(IX,IY,IZ-1),BUFI(IX,IY,IZ+1),1,RZ) XSHNEW = IXSH - RX YSHNEW = IYSH - RY ZSHNEW = IZSH - RZ ELSEIF (NROW .GT. 1) THEN !------------------------- IMAGE IF (IX .LE. 1 .OR. IX .GE. NSAM .OR. & IY .LE. 1 .OR. IY .GE. NROW) RETURN C PARABOLIC FIT TO THE 3X3 NEIGHBORHOOD OF PEAK CALL PARABL(BUFI(IX-1:IX+1, & IY-1:IY+1, 1), XSHNEW,YSHNEW, PEAKV) XSHNEW = IXSH - XSHNEW YSHNEW = IYSH - YSHNEW ZSHNEW = 0.0 ELSE ! ----------------------------- ROW IF (IX .LE. 1 .OR. IX .GE. NSAM) RETURN C BINARY INTERPOLATION TO FIND SUB-PIXEL PEAK QMAX = BUFI(IX,IY,IZ) ! CENTRAL ELEMENT CALL PKSR3_SUB(QMAX,BUFI(IX-1,IY,IZ),BUFI(IX+1,IY,IZ),1,RX) XSHNEW = IXSH - RX YSHNEW = 0.0 ZSHNEW = 0.0 ENDIF C NORMALIZATION DONE IN CALLER IF DESIRED IRTFLG = 0 RETURN END C ------------------------- UNUSED BELOW HERE ------------------------- #ifdef NEVER #ifdef DEBUG write(6,*)' ' write(6,*)'IGOX,IENDX,IGOY,IENDY ',IGOX,IENDX,IGOY,IENDY write(6,*)' ' write(6,*)'ILOCS: ',ilocs write(6,*)'IX,IY,IZ: ',ix,iy write(6,*)'IXSH,IYSH: ',ixsh,iysh write(6,*)'XSHNEW,YSHNEW:',xshnew,yshnew write(6,*)'PEAKV: ',peakv #endif #ifdef DEBUG write(6,*)' ' write(6,*)'XSHNEW,YSHNEW:',xshnew,yshnew #endif write(6,*)' ' write(6,*)'IGOX,IENDX,IGOY,IENDY,IGOZ,IENDZ ' write(6,*) IGOX,IENDX,IGOY,IENDY,IGOZ,IENDZ write(6,*)' ' write(6,*)'ILOCS: ',ilocs write(6,*)'IX,IY,IZ: ',ix,iy,iz write(6,*)'IXSH,IYSH,IZSH: ',ixsh,iysh,izsh write(6,*)' ' write(6,*)'RX,RY,RZ:',rx,ry,rz write(6,*)'XSHNEW,YSHNEW,ZSHNEW:',xshnew,yshnew,zshnew write(6,*)' ' write(6,*)'peakv: ',peakv nnn = nsam*nrow*nslice lnn = lsc*nrow*nslice write(6,*)'nsam*nrow*nslice:',nnn write(6,*)'lsc*nrow*nslice:',lnn pdnnn = peakv / float(nnn) pdlnn = peakv / float(lnn) write(6,*)'pdnnn,pdlnn:',pdnnn,pdlnn #endif #ifdef NEVER call opfilec(0,.false.,'jnk002',98,'O',itype, & nsam,nrow,nslice,maxim,' ',.true.,irtflg) call redvol(98,nsam,nrow, 1,nslice, bufs,irtflg) close(98) write(6,*) ' bufs filled with jnk002 ' maxim = 0 itype = -21 if (mod(nsam,2).eq.0) itype = -22 call opfilec(0,.false.,'jnkapsccffti',98,'U',itype, & lsc,nrow,nslice,maxim,' ',.true.,irtflg) call wrtvol(98,lsc,nrow, 1,nslice, bufi,irtflg) close(98) itype = -21 if (mod(nsam,2).eq.0) itype = -22 call opfilec(0,.false.,'jnkapsccfftr',98,'U',itype, & lsc,nrow,nslice,maxim,' ',.true.,irtflg) call wrtvol(98,lsc,nrow,1,nslice, bufr,irtflg) close(98) maxim = 0 itype = 3 call opfilec(0,.false.,'jnknewpkpad',98,'U',itype, & nsam,nrow,nslice,maxim,' ',.false.,irtflg) c this only writes first lsc values from each row call writev(98,bufi,lsc,nrow, nsam,nrow,nslice) close(98) #endif