C++********************************************************************* C C ORMD_P.F C 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 ORMD_P(NSAM,NROW,LSAM,LROW, C NRING,LCIRC,MAXRIN,NUMR,X,NPEAK,MODE,LUNEXP,LUNREF) C C PARAMETERS: NSAM,NROW USED SIZE C LSAM,LROW ACTUAL SIZE C C LUNEXP,LUNREF,LUNDOC IO UNITS INPUT C C C23456789012345678901234567890123456789012345678901234567890123456789012 C--********************************************************************* SUBROUTINE ORMD_P(NSAM,NROW, LSAM,LROW, & NRING,LCIRC,MAXRIN,NUMR,X,NPEAK,MODE,FFTW_PLANS, & LUNEXP,LUNREF,LUNDOC) INCLUDE 'CMBLOCK.INC' INTEGER, PARAMETER :: NDLI=2 REAL :: DLIST(NDLI) REAL :: X(NSAM,NROW),BUFIN(LSAM) REAL, ALLOCATABLE, DIMENSION(:) :: CIRC,CIRCREF INTEGER :: NUMR(3,NRING) REAL :: TEMP(MAXRIN+2) REAL :: PEAK(2,NPEAK) DOUBLE PRECISION :: T7(-3:3) CHARACTER(LEN=1) :: MODE REAL :: WR(NRING) #ifdef SP_LIBFFT C TT FOR SGI FFT USE REAL, ALLOCATABLE, DIMENSION(:) :: TT #endif DOUBLE PRECISION :: QMAX REAL :: POS_MAX JACUP = 0 C FIND WEIGHTS FOR TRANSFORMED CIRC RINGS CALL RINGWE_NEW(WR,NUMR,NRING,MAXRIN) ALLOCATE (CIRC(LCIRC), CIRCREF(LCIRC), STAT=IRTFLG) IF (IRTFLG.NE.0) THEN CALL ERRT(46,'ORMD_P, CIRC',2*LCIRC) RETURN ENDIF C FIND CENTRAL WINDOW SIZE LQ = LROW/2+1 LR1 = (NROW-1)/2 LR2 = LQ+LR1 LR1 = LQ-LR1 LQ = LSAM/2+1 LS1 = (NSAM-1)/2 LS2 = LQ+LS1 LS1 = LQ-LS1 C READ IN CENTRAL WINDOW FROM THE WHOLE EXP. IMAGE DO K2=LR1,LR2 CALL REDLIN(LUNEXP,BUFIN,LSAM,K2) DO K3=LS1,LS2 X(K3-LS1+1,K2-LR1+1) = BUFIN(K3) ENDDO ENDDO CLOSE(LUNEXP) 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, X, .FALSE., & MODE,NUMR,NRING,LCIRC, WR,FFTW_PLANS, & CIRC,IRTFLG) C READ IN CENTRAL WINDOW FROM THE WHOLE REFERENCE IMAGE DO K2=LR1,LR2 CALL REDLIN(LUNREF,BUFIN,LSAM,K2) DO K3=LS1,LS2 X(K3-LS1+1,K2-LR1+1) = BUFIN(K3) ENDDO ENDDO CLOSE(LUNREF) C EXTRACT REF. IMAGE RINGS, NORMALIZE & FFT THEM WR(1) = 0.0 ! ONLY WEIGHT ONCE CALL APRINGS_ONE_NEW(NSAM,NROW, CNS2,CNR2, X, .FALSE., & MODE,NUMR,NRING,LCIRC, WR,FFTW_PLANS, & CIRCREF,IRTFLG) #ifdef SP_LIBFFT C CREATE TT FOR SGI FFT USE LENTT = NUMR(3,NRING) + 15 ALLOCATE (TT(LENTT), STAT=IRTFLG) IF (IRTFLG.NE.0) THEN CALL ERRT(46,'ORMD_P, TT ',LENTT) RETURN ENDIF CALL DZFFT1DUI(LENTT-15,TT) #endif CALL CROSRNG_NEW(CIRC,CIRCREF,LCIRC, NRING,MAXRIN,NUMR, & TEMP,QMAX,POS_MAX, TT,.FALSE., FFTW_PLANS,MODE) IF (NPEAK .LE. 1) THEN RANGNEW = ANGMOR_NEW(POS_MAX,MODE,MAXRIN) FQMAX = 2.0 * QMAX / MAXRIN / NRING WRITE(NOUT,2700)RANGNEW,FQMAX 2700 FORMAT(' Angle = ',F10.4,' Peak Height = ',G12.5) C ONE PEAK, FIRST REGISTER IS ANGLE, SECOND IS PEAK VALUE CALL REG_SET_NSEL(1,2,RANGNEW,FQMAX, 0.0,0.0,0.0,IRTFLG) ELSE C FIND SPECIFIED NUMBER OF PEAKS DO K2=1,NPEAK PEAK(1,K2) = -HUGE(PEAK) PEAK(2,K2) = -1.0 ENDDO DO J=1,MAXRIN J2 = MOD(J+MAXRIN-1,MAXRIN)+1 J1 = MOD(J+MAXRIN-2,MAXRIN)+1 J3 = MOD(J+MAXRIN,MAXRIN)+1 IF (TEMP(J2) .GT. TEMP(J1) .AND. & TEMP(J2) .GT. TEMP(J3)) THEN DO K2=1,NPEAK IF (TEMP(J2) .GT. PEAK(1,K2)) THEN IF (NPEAK .GT. 1) THEN C MOVE THIS PEAK IN LIST DO K3=NPEAK,K2+1,-1 PEAK(1,K3) = PEAK(1,K3-1) PEAK(2,K3) = PEAK(2,K3-1) ENDDO ENDIF PEAK(1,K2) = TEMP(J2) PEAK(2,K2) = J EXIT ENDIF ENDDO ENDIF ENDDO C CONVERT TO ANGLES AND DO THE INTERPOLATION DO K2=1,NPEAK JTOT = PEAK(2,K2) NPT = K2 IF (JTOT .EQ. -1) EXIT DO K3=-3,3 J = MOD(JTOT+K3+MAXRIN-1,MAXRIN) + 1 T7(K3) = TEMP(J) ENDDO C SUB-PIXEL INTERPOLATION CALL PRB1D(T7,7,POS) QMAX = 2.0 * PEAK(1,K2) / MAXRIN/ MAXRIN / NRING ! HEIGHT RANG = PEAK(2,K2) + POS ! LOCATION C CONVERT PEAK LOCATION TO ANGLE PEAK(2,K2) = ANGMOR_NEW(RANG,MODE,MAXRIN) !write(6,90) j,pos,rang,peak(2,k2) 90 format( ' j,pos,rang,angle: ', i5, f8.2, f8.2, f8.2) DLIST(1) = PEAK(2,K2) ! ANGLE DLIST(2) = QMAX ! HEIGHT CALL LUNDOCWRTDAT(LUNDOC,K2,DLIST,NDLI,IRTFLG) WRITE(NOUT,2701) PEAK(2,K2),QMAX 2701 FORMAT(' Angle = ',F10.4,' Peak Height = ',G12.5) ENDDO C MULTIPLE PEAKS, FIRST REGISTER IS THE NUMBER OF PEAKS SAVED CALL REG_SET_NSEL(1,1,FLOAT(NPT), 0.0,0.0,0.0,0.0,IRTFLG) ENDIF 9999 IF (ALLOCATED(CIRC)) DEALLOCATE (CIRC) IF (ALLOCATED(CIRCREF)) DEALLOCATE (CIRCREF) END C -------------------- ANGMOR_NEW ------------------------------ REAL FUNCTION ANGMOR_NEW(RKK,MODE,MAXRIN) REAL, INTENT(IN) :: RKK CHARACTER(LEN=1), INTENT(IN) :: MODE INTEGER, INTENT(IN) :: MAXRIN IF (MODE .EQ. 'H') THEN ANGMOR_NEW = (RKK - 1.0) / MAXRIN * 180.0 IF (ANGMOR_NEW .GT. 0.0) ANGMOR_NEW = 360.0 - ANGMOR_NEW ANGMOR_NEW = MOD(ANGMOR_NEW + 180.0,180.0) ELSEIF (MODE .EQ. 'F') THEN ANGMOR_NEW = (RKK - 1.0) / MAXRIN * 360.0 IF (ANGMOR_NEW .GT. 0.0) ANGMOR_NEW = 360.0 - ANGMOR_NEW ANGMOR_NEW = MOD(ANGMOR_NEW+360.0, 360.0) ENDIF END