C++********************************************************************* C C APRINGS.F USED CMLIMIT AUG 00 ARDEAN LEITH C ADDED REF_CIRC FILE APR 01 ARDEAN LEITH C NORMASS -> NORMAS OCT 01 ARDEAN LEITH C PROMPTS JAN 02 ARDEAN LEITH C AUTO REF RINGS FILE DEC 04 ARDEAN LEITH C SPIDER REF RINGS FILE FEB 05 ARDEAN LEITH C INULL = lnblnkn(SCRFILE) APR 05 ARDEAN LEITH C REWRITE FOR SPEED MAR 08 ARDEAN LEITH C BCAST_MPI NOV 08 ARDEAN LEITH C OUTPUT TO NOUT AUG 09 ARDEAN LEITH 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 APRINGS(ILIST,NUMREF,NSAM,NROW,NDUM,NDUM, C NRING,LCIRC,NUMR,MODE, C REFPAT,LUNREF,CIRCREF,CIRCREF_IN_CORE, C LUNRING,SCRFILE,IRTFLG) C C PARAMETERS: ILIST LIST OF IMAGE FILE NUMBERS (INPUT) C NIMREF NO. OF IMAGES (INPUT) C NSAM,NROW IMAGE DIMENSIONS (INPUT) C REFPAT IMAGE SERIES FILE TEMPLATE (INPUT) C LUNREF IMAGE FILE IO UNIT (INPUT) C CIRCREF OUTPUT ARRAY (OUTPUT) C C CIRCREF_IN_CORE NO OUTPUT ARRAY FLAG (OUTPUT) C LUNRING REF-RINGS FILE IO UNIT (INPUT) C SCRFILE REF-RINGS FILE (INPUT) C IRTFLG ERROR FLAG (OUTPUT) C C NOTE: MOST MEMORY DEMAND DEPENDENT ON LCIRC & NUMREF. LCIRC IS THE C TOTAL LENGTH OF ARRAY THAT HOLDS THE CIRCULAR RINGS, SO IT IS C DEPENDENT ON NUMBER OF RINGS AND THEIR RADIUS. ARRAY ALLOCATED C IS: CIRCREF(LCIRC,NUMREF) ANOTHER SMALL ALLOCATED ARRAY IS: C A(NSAM,NROW,NUMTH) C C23456789012345678901234567890123456789012345678901234567890123456789012 C--********************************************************************* C -------------------- APRINGS_NEW ---------------------------- SUBROUTINE APRINGS_NEW(ILIST,NUMREF, & NSAM,NROW, & NRING,LCIRC,NUMR,MODE,FFTW_PLANS, & REFPAT,LUNREF,CIRCREF,CIRCREF_IN_CORE, & LUNRING,SCRFILE,IRTFLG) INCLUDE 'CMBLOCK.INC' INCLUDE 'CMLIMIT.INC' CHARACTER(LEN=1) :: MODE INTEGER, DIMENSION(NUMREF) :: ILIST INTEGER, DIMENSION(3,NRING) :: NUMR REAL, DIMENSION(LCIRC,*) :: CIRCREF CHARACTER (LEN=MAXNAM) :: FILNAM CHARACTER (LEN=*) :: SCRFILE,REFPAT LOGICAL :: USEREFFILE,CIRCREF_IN_CORE LOGICAL :: ISOPEN LOGICAL :: WINDOW C FFTW_PLANS CONTAINS POINTERS TO STRUCTURES INTEGER*8 :: FFTW_PLANS(*) C ALLOCATABLE ARRAYS REAL, ALLOCATABLE, DIMENSION(:,:) :: XIM CALL SET_MPI(ICOMM,MYPID,MPIERR) ! SETS ICOMM AND MYPID C SEE IF THIS FILE EXISTS, (RETURNS EX, ISOPEN, LUNOP) C SEE IF SCRATCH "FILE" EXISTS (MAY BE INCORE FILE) ILOCAT = INDEX(SCRFILE,'@') INULL = lnblnkn(SCRFILE) IF (INULL .GT. 0 .AND. SCRFILE(1:1) .NE. '_' .AND. & ILOCAT .EQ. 0) THEN C ADD EXTENSION TO PHYSICAL FILENAME CALL FILNAMANDEXT(SCRFILE,DATEXC,FILNAM,NLET,.TRUE.,IRTFLG) IF (IRTFLG .NE. 0) RETURN ELSE FILNAM = SCRFILE NLET = lnblnk(FILNAM) ENDIF #ifdef USE_MPI IF (MYPID .LE. 0) THEN INQUIRE(FILE=FILNAM,EXIST=USEREFFILE,OPENED=ISOPEN, & NUMBER=LUNOP,IOSTAT=IRTFLG) ENDIF CALL BCAST_MPI('APRINGS','USEREFFILE',USEREFFILE,1, 'L',ICOMM) #else USEREFFILE = .FALSE. IF (INULL .GT. 0) THEN CALL INQUIREIF1(LUNRING,FILNAM,DUM,NDUM,USEREFFILE,ISOPEN, & LUNOP,INLNED,IMGNUM,IRTFLG) ENDIF #endif IF (USEREFFILE .AND. MYPID .LE. 0) THEN WRITE(NOUT,*) ' Using existing reference rings file: ', & FILNAM(1:NLET) ELSEIF (MYPID .LE. 0) THEN WRITE(NOUT,*) ' No existing reference rings file: ', & FILNAM(1:NLET) ENDIF C NUMBER OF OMP THREADS CALL GETTHREADS(NUMTH) IF (CIRCREF_IN_CORE .AND. .NOT. USEREFFILE) THEN C CALCULATE REF. RINGS DATA AND FILL CIRCREF ARRAY WITH IT CALL APRINGS_FILL_NEW(ILIST,NUMREF, NSAM,NROW,NUMTH, & NRING,LCIRC,NUMR,MODE,FFTW_PLANS, & REFPAT,LUNREF, & CIRCREF,NUMREF,0,IRTFLG) IF (MYPID .LE. 0) & WRITE(NOUT,*) ' Created incore reference rings' ELSEIF (CIRCREF_IN_CORE .AND. USEREFFILE) THEN C READ EXISTING REF RINGS FILE AND FILL CIRCREF ARRAY WITH ITS DATA C OPEN EXISTING REFERENCE RINGS FILE MAXIM = 0 CALL OPFILEC(0,.FALSE.,SCRFILE,LUNRING,'O',IFORM, & LCIRCT,NUMREFT,NSLICE,MAXIM,' ', .TRUE.,IRTFLG) IF (IRTFLG .NE. 0) RETURN IF (LCIRCT .NE. LCIRC .OR. NUMREFT .NE. NUMREF) THEN CALL ERRT(101,'REF. RINGS FILE HAS WRONG SIZE',NE) IRTFLG = 1 RETURN ENDIF C FILL CIRCREF WITH EXISTING RINGS DATA FROM FILE AND RETURN. DO IREF=1,NUMREF CALL REDLIN(LUNRING,CIRCREF(1,IREF),LCIRC,IREF) ENDDO IF (MYPID .LE. 0) THEN WRITE(NOUT,*) ' Loaded reference rings file incore' ENDIF ELSEIF (.NOT. CIRCREF_IN_CORE .AND. USEREFFILE) THEN C OPEN EXISTING REF. RINGS FILE TO READ REF RINGS DATA MAXIM = 0 CALL OPFILEC(0,.FALSE.,SCRFILE,LUNRING,'O',IFORM, & LCIRCT,NUMREFT,NSLICE,MAXIM,' ', .FALSE.,IRTFLG) IF (IRTFLG .NE. 0) RETURN IF (LCIRCT .NE. LCIRC .OR. NUMREFT .NE. NUMREF) THEN CALL ERRT(101,'REF. RINGS FILE HAS DIFFERENT SIZE',NE) IRTFLG = 1 RETURN ENDIF ELSEIF (.NOT. CIRCREF_IN_CORE .AND. .NOT. USEREFFILE) THEN C CALCULATE REF. RINGS DATA AND FILL REF RINGS FILE WITH DATA C CREATE REFERENCE RINGS FILE FOR OUTPUT NSL = 1 CALL OPFILEC(0,.FALSE.,SCRFILE,LUNRING,'B',IFORM, & LCIRC,NUMREF,NSL,MAXIM,' ', .FALSE.,IRTFLG) IF (IRTFLG .NE. 0) RETURN IF (MYPID .LE. 0) WRITE (NOUT,90) SCRFILE 90 FORMAT (' Created reference rings file: ',A ) C FILL REF_CIRC FILE WITH RINGS DATA CALL APRINGS_FILL_NEW(ILIST,NUMREF, NSAM,NROW,NUMTH, & NRING,LCIRC,NUMR,MODE,FFTW_PLANS, & REFPAT,LUNREF, & CIRCREF,NUMTH,LUNRING,IRTFLG) IF (MYPID .LE. 0) WRITE (NOUT,92) SCRFILE 92 FORMAT (' Filled reference rings file: ',A ) ENDIF END C --------------------- APRINGS_FILL_NEW -------------------------- SUBROUTINE APRINGS_FILL_NEW(ILIST,NUMREF, & NSAM,NROW,NUMTH, & NRING,LCIRC,NUMR,MODE,FFTW_PLANS, & REFPAT,LUNREF, & CIRCREF,ICORE,LUNRING,IRTFLG) INCLUDE 'CMBLOCK.INC' INCLUDE 'CMLIMIT.INC' CHARACTER(LEN=1) :: MODE CHARACTER (LEN=*) :: REFPAT LOGICAL :: WINDOW C EXTERNAL ARRAYS INTEGER,DIMENSION(NUMREF) :: ILIST INTEGER,DIMENSION(3,NRING) :: NUMR REAL,DIMENSION(LCIRC,ICORE) :: CIRCREF C FFTW_PLANS CONTAINS POINTERS TO STRUCTURES INTEGER*8 :: FFTW_PLANS(*) C AUTOMATIC ARRAYS REAL,DIMENSION(NRING) :: WR C ALLOCATABLE ARRAYS REAL, ALLOCATABLE, DIMENSION(:,:,:) :: XIM CALL SET_MPI(ICOMM,MYPID,MPIERR) ! SETS ICOMM AND MYPID IRTFLG = 0 IF (MYPID .LE. 0) & WRITE(NOUT,*) ' Filling reference rings with threads:',NUMTH MAXRIN = NUMR(3,NRING) #ifdef SP_LIBFFTW3 MAXRIN = NUMR(3,NRING) - 2 #endif C RINGWE RETURNS WR WEIGHTS CALL RINGWE_NEW(WR,NUMR,NRING,MAXRIN) IF (MODE .EQ. 'H') WR = WR * 0.5 ALLOCATE(XIM(NSAM,NROW,NUMTH), STAT=IRTFLG) IF (IRTFLG.NE.0) THEN CALL ERRT(46,'XIM',NSAM*NROW*NUMTH) RETURN ENDIF C CALCULATE CENTER FOR APRINGS CNS2 = NSAM / 2 + 1 CNR2 = NROW / 2 + 1 C PREPARE CIRCULAR RINGS DATA FOR ALL REFERENCE IMAGES C LOOP OVER ALL REF. IMAGES, CHUNKS OF NUMTH DO IGO=1,NUMREF,NUMTH IEND = MIN(IGO+NUMTH-1,NUMREF) C LOAD SET OF REFERENCE IMAGES (THERE ARE NUMTH IMAGES IN SET) CALL AP_GETDAT(ILIST,NUMREF,NSAM,NROW,NSAM,NROW, & NUMTH,REFPAT,LUNREF, IGO,IEND, & 1,NROW,1,NSAM, XIM, IRTFLG) IF (IRTFLG .NE. 0) GOTO 9999 c$omp parallel do private(IMI,IPT,IT) DO IMI=IGO,IEND IT = IMI - IGO + 1 ! XIM INDEX IPT = IMI ! CIRCREF INDEX IF (LUNRING .GT. 0) IPT = IT C CONVERT XIM TO POLAR RINGS, FFT, & WEIGHT THE RINGS CALL APRINGS_ONE_NEW(NSAM,NROW, CNS2,CNR2, XIM(1,1,IT), & .FALSE.,MODE,NUMR,NRING,LCIRC,WR,FFTW_PLANS, & CIRCREF(1,IPT),IRTFLG) c write(6,*) 'circref(1,',ipt,'):',CIRCREF(1,IPT),it ENDDO IF (LUNRING .GT. 0) THEN C SAVE CIRCREF IN FILE OPENED ON LUNRING DO IMI=IGO,IEND IV = IMI - IGO + 1 CALL WRTLIN(LUNRING,CIRCREF(1,IV),LCIRC,IMI) ENDDO ENDIF ENDDO C DEALLOCATE ARRAY 9999 IF (ALLOCATED(XIM)) DEALLOCATE(XIM) END C --------------------- APRINGS_ONE_NEW -------------------------- SUBROUTINE APRINGS_ONE_NEW(NSAM,NROW, CNS2,CNR2, XIM, OMP_FRNG, & MODE,NUMR,NRING,LCIRC,WR,FFTW_PLANS, & CIRC,IRTFLG) C EXTERNAL ARRAYS REAL, INTENT(IN) :: XIM(NSAM,NROW) REAL, INTENT(OUT) :: CIRC(LCIRC) INTEGER, INTENT(IN) :: NUMR(3,NRING) REAL, INTENT(IN) :: WR(NRING) LOGICAL, INTENT(IN) :: OMP_FRNG CHARACTER(LEN=1) :: MODE C FFTW_PLANS CONTAINS POINTERS TO STRUCTURES INTEGER*8, INTENT(IN) :: FFTW_PLANS(*) DOUBLE PRECISION :: AVT,VRINV LOGICAL, PARAMETER :: NEWFFTW = .TRUE. LOGICAL, PARAMETER :: SPIDER_SIGN = .FALSE. C FIND PARAMETERS TO NORMALIZE UNDER THE MASK, C TRIED DOING THIS COMPLETELY ON THE C POLAR RINGS BUT IT GIVES SOME DIFFERENT REF. CHOICES. al C CALCULATE DIMENSIONS FOR NORMALIZING MASK NSB = -CNS2 NSE = NSB + NSAM - 1 NRB = -CNR2 NRE = NRB + NROW - 1 CALL NORMASC(XIM, NSB,NSE,NRB,NRE, NUMR,NUMR(1,NRING), & AVT,VRINV) C INTERPOLATE INTO POLAR COORDINATES & APPLY NORMALIZATION C CREATING CIRC (RADIAL IMAGE CIRCLES) FOR THIS IMAGE POSITION CALL ALRQ_MS_NEW(XIM,NSAM,NROW, CNS2,CNR2, & NUMR,CIRC,LCIRC, NRING,MODE, NEWFFTW,OMP_FRNG, & AVT,VRINV) C FOURIER TRANSFORM CIRC RINGS CALL FRNGS_NEWT(CIRC,LCIRC, NUMR,NRING, SPIDER_SIGN, & FFTW_PLANS, OMP_FRNG) IF (WR(1) .GT. 0.0) THEN C WEIGHT TRANSFORMED CIRC RINGS USING FACTORS FROM: WR CALL APPLYWS_NEW(CIRC,LCIRC,NUMR,WR,NRING) ENDIF IRTFLG = 0 END C --------------- APRINGS_INIT_PLANS ----------------------------- SUBROUTINE APRINGS_INIT_PLANS(NUMR,NRING, & FFTW_PLANS,NPLANS,IRTFLG) C INITIALIZE REV. FFTW3 PLAN(S) FOR USE WITHIN OMP || INTEGER, INTENT(IN) :: NUMR(3,NRING) INTEGER, INTENT(IN) :: NRING INTEGER, INTENT(IN) :: NPLANS INTEGER, INTENT(OUT) :: IRTFLG #ifdef SP_LIBFFTW3 C FFTW_PLANS IS AN ARRAY OF POINTERS TO STRUCTURES INTEGER*8, INTENT(OUT) :: FFTW_PLANS(NPLANS) FFTW_PLANS = 0 IRTFLG = 1 NUMTH = 1 LENO = -1 DO I=1,NRING LEN = NUMR(3,I) - 2 ! LENGTH OF RING IF (LEN .GT. LENO) THEN C CREATE PLAN FOR THIS NEW RING LENGTH INDX = LOG2(LEN) - 1 !write(6,*) 'plans; i,len,indx: ',i,len,indx IF (INDX .LT. 2 .OR. INDX .GT. NPLANS) THEN WRITE(NOUT,*) ' CREATING FFTW_PLANS(',INDX,')' CALL ERRT(101,'FFTW_PLANS LIMITED TO 2...',NPLANS) RETURN ENDIF CALL FFTW3_MAKEPLAN(LEN,1,1, NUMTH, & FFTW_PLANS(INDX),+1,IRTFLG) IF (IRTFLG .NE. 0) RETURN LENO = LEN ENDIF ENDDO C CAN USE INDX = 1 FOR REVERSE PLAN AS IT IS NEVER NEEDED ELSEWISE C INITIALIZE REV. FFTW3 PLAN FOR USE WITHIN OMP || MAXRIN = NUMR(3,NRING) - 2 CALL FFTW3_MAKEPLAN(MAXRIN,1,1, NUMTH, & FFTW_PLANS(1),-1,IRTFLG) IF (IRTFLG .NE. 0) RETURN #endif IRTFLG = 0 END C --------------------- APRINGS_ONE -------------------------- C still used in: writpro_n.f should be replaced !!!!!!!!!!!!!!!!! SUBROUTINE APRINGS_ONE(NSAM,NROW,WR, & XIM,CIRC,MODE,NUMR,NRING,LCIRC, & IRTFLG) INCLUDE 'CMLIMIT.INC' INCLUDE 'CMBLOCK.INC' C EXTERNAL ARRAYS REAL, DIMENSION(NSAM,NROW) :: XIM REAL, DIMENSION(LCIRC) :: CIRC INTEGER, DIMENSION(3,NRING) :: NUMR REAL,DIMENSION(NRING) :: WR CHARACTER(LEN=1) :: MODE DOUBLE PRECISION :: PI,DFI C DEFAULT RETURN FLAG IRTFLG = 1 C CALCULATE DIMENSIONS FOR NORMASS NSB = -NSAM/2 NSE = -NSB-1+MOD(NSAM,2) NRB = -NROW/2 NRE = -NRB-1+MOD(NROW,2) C NORMALIZE UNDER THE MASK, TRIED DOING THIS ON THE C POLAR RINGS BUT IT GIVES DIFFERENT ANSWERS. al CALL NORMASS(XIM,NSB,NSE,NRB,NRE,NUMR,NUMR(1,NRING)) MAXRIN = NUMR(3,NRING) NS2 = NSAM / 2 + 1 NR2 = NROW / 2 + 1 FNR2 = NR2 FNS2 = NS2 PI = 2 * DATAN(1.0D0) C CONVERT WINDOW FROM IMAGE INTO POLAR COORDINATES IF (MODE .EQ. 'F') THEN LTF = 4 DO I=1,NRING C RADIUS OF THE RING INR = NUMR(1,I) YQ = INR LT = NUMR(3,I) / LTF LT2 = LT + LT LT3 = LT2 + LT DFI = PI / LT KCIRC = NUMR(2,I) CIRC(KCIRC) = XIM(NS2, NR2+INR) CIRC(LT+KCIRC) = XIM(NS2+INR, NR2) CIRC(LT2+KCIRC) = XIM(NS2, NR2-INR) CIRC(LT3+KCIRC) = XIM(NS2-INR, NR2) DO J=1,LT - 1 FI = DFI * J X = SIN(FI) * YQ Y = COS(FI) * YQ JT = J + KCIRC CIRC(JT) = QUADRI(FNS2+X, FNR2+Y, NSAM,NROW,XIM) CIRC(JT+LT) = QUADRI(FNS2+Y, FNR2-X, NSAM,NROW,XIM) CIRC(JT+LT2) = QUADRI(FNS2-X, FNR2-Y, NSAM,NROW,XIM) CIRC(JT+LT3) = QUADRI(FNS2-Y, FNR2+X, NSAM,NROW,XIM) ENDDO ENDDO ELSEIF (MODE .EQ. 'H') THEN LTF = 2 DO I=1,NRING C RADIUS OF THE RING INR = NUMR(1,I) YQ = INR LT = NUMR(3,I) / LTF DFI = PI / LT KCIRC = NUMR(2,I) CIRC(KCIRC) = XIM(NS2, NR2+INR) CIRC(LT+KCIRC) = XIM(NS2+INR, NR2) DO J=1,LT - 1 FI = DFI * J X = SIN(FI) * YQ Y = COS(FI) * YQ JT = J + KCIRC CIRC(JT) = QUADRI(FNS2+X, FNR2+Y, NSAM,NROW,XIM) CIRC(JT+LT) = QUADRI(FNS2+Y, FNR2-X, NSAM,NROW,XIM) ENDDO ENDDO ENDIF CALL FRNGS(CIRC,LCIRC,NUMR,NRING) C WEIGHT CIRC USING WR IF (WR(1) .GT. 0.0) THEN CALL APPLYWS(CIRC,LCIRC,NUMR,WR,NRING,MAXRIN) ENDIF IRTFLG = 0 END #ifdef DEBUGNEVER write(6,*) ' circ(1..50) ------ ' write(6,901) (circ(ii),ii=1,50) 901 format(5f11.2) write(6,*) ' circ(1000..1050) ------ ' write(6,901) (circ(ii),ii=1000,1050) #endif #ifdef NEVER C FOURIER OF CIRC DO I=1,NRING L = LOG2(NUMR(3,I)) CALL FFTR_Q(CIRC(NUMR(2,I)),L) ENDDO C WEIGHT CIRC USING WR IF (WR(1) .GT. 0.0) THEN CALL APPLYWS(CIRC,LCIRC,NUMR,WR,NRING,MAXRIN) ENDIF IF (WR(1) .GT. 0.0) THEN C WEIGHT CIRC USING WR DO I=1,NRING NUMR3I = NUMR(3,I) NUMR2I = NUMR(2,I) W = WR(I) CIRC(NUMR2I) = CIRC(NUMR2I)*W IF (NUMR3I .EQ. MAXRIN) THEN CIRC(NUMR2I+1) = CIRC(NUMR2I+1) * 0.5*W ELSE CIRC(NUMR2I+1) = CIRC(NUMR2I+1) * W ENDIF DO J=3,NUMR3I JC = J + NUMR2I-1 CIRC(JC) = CIRC(JC)*W ENDDO ENDDO ENDIF c write(6,*) ' unweighted (1-100):',circ(1),circ(100) c avtt = 0 c vrtt = 0 c do i =1,lcirc c avtt = avtt + circ(i) c vrtt = vrtt + circ(i) * circ(i) c enddo c write(6,*) 'n,avtt,vrtt:',n,avtt,vrtt c write(6,*) 'n,sumavt,sumvr:',n,avt,vr c avtn = avt / n c vrn = dsqrt((vr - n * avtn * avtn) / (n-1)) c tt = vr - n * avt * avt c write(6,*) 'n,avtn,vrn:',n,avtn,vrn,tt cc AVT = AVT / N C MULTIPLICATION IS FASTER cc VRINV = 1.0 / (DSQRT((VR - N * AVT * AVT) / (N-1))) C NORMALIZE CIRC cc CIRC = (CIRC - AVT) * VRINV C unused development idea C --------------------- APRINGS_TWO -------------------------- SUBROUTINE APRINGS_TWO(NSAM,NROW,WR, & XIM,CIRC,MODE,NUMR,NRING,LCIRC, & IRTFLG) C EXTERNAL ARRAYS REAL, DIMENSION(NSAM,NROW) :: XIM REAL, DIMENSION(LCIRC) :: CIRC INTEGER, DIMENSION(3,NRING) :: NUMR REAL,DIMENSION(NRING) :: WR CHARACTER(LEN=1) :: MODE DOUBLE PRECISION :: PI,DFI C DEFAULT RETURN FLAG IRTFLG = 1 C CALCULATE DIMENSIONS FOR NORMASS NSB = -NSAM/2 NSE = -NSB-1+MOD(NSAM,2) NRB = -NROW/2 NRE = -NRB-1+MOD(NROW,2) C NORMALIZE UNDER THE MASK, TRIED DOING THIS ON THE C POLAR RINGS BUT IT GIVES DIFFERENT ANSWERS. al CALL NORMASS(XIM,NSB,NSE,NRB,NRE,NUMR,NUMR(1,NRING)) MAXRIN = NUMR(3,NRING) NS2 = NSAM / 2 + 1 NR2 = NROW / 2 + 1 FNR2 = NR2 FNS2 = NS2 PI = 2 * DATAN(1.0D0) C CONVERT WINDOW FROM IMAGE INTO POLAR COORDINATES IF (MODE .EQ. 'F') THEN LTF = 4 DO I=1,NRING C RADIUS OF THE RING INR = NUMR(1,I) YQ = INR LT = NUMR(3,I) / LTF LT2 = LT + LT LT3 = LT2 + LT DFI = PI / LT KCIRC = NUMR(2,I) CIRC(KCIRC) = XIM(NS2, NR2+INR) CIRC(LT+KCIRC) = XIM(NS2+INR, NR2) CIRC(LT2+KCIRC) = XIM(NS2, NR2-INR) CIRC(LT3+KCIRC) = XIM(NS2-INR, NR2) DO J=1,LT - 1 FI = DFI * J X = SIN(FI) * YQ Y = COS(FI) * YQ JT = J + KCIRC CIRC(JT) = QUADRI(FNS2+X, FNR2+Y, NSAM,NROW,XIM) CIRC(JT+LT) = QUADRI(FNS2+Y, FNR2-X, NSAM,NROW,XIM) CIRC(JT+LT2) = QUADRI(FNS2-X, FNR2-Y, NSAM,NROW,XIM) CIRC(JT+LT3) = QUADRI(FNS2-Y, FNR2+X, NSAM,NROW,XIM) ENDDO ENDDO ELSEIF (MODE .EQ. 'H') THEN LTF = 2 DO I=1,NRING C RADIUS OF THE RING INR = NUMR(1,I) YQ = INR LT = NUMR(3,I) / LTF DFI = PI / LT KCIRC = NUMR(2,I) CIRC(KCIRC) = XIM(NS2, NR2+INR) CIRC(LT+KCIRC) = XIM(NS2+INR, NR2) DO J=1,LT - 1 FI = DFI * J X = SIN(FI) * YQ Y = COS(FI) * YQ JT = J + KCIRC CIRC(JT) = QUADRI(FNS2+X, FNR2+Y, NSAM,NROW,XIM) CIRC(JT+LT) = QUADRI(FNS2+Y, FNR2-X, NSAM,NROW,XIM) ENDDO ENDDO ENDIF C FOURIER OF CIRC CALL FRNGS_NEW(CIRC,LCIRC,NUMR,NRING,MAXRIN) C WEIGHT CIRC USING WR IF (WR(1) .GT. 0.0) THEN CALL APPLYWS(CIRC,LCIRC,NUMR,WR,NRING,MAXRIN) ENDIF IRTFLG = 0 END #endif