C++*********************************************************************
C
C    RFACTSDO.F              ADDED FOURIER INPUT JULY 2000 ARDEAN LEITH
C                            OPFILEC             FEB  03   ARDEAN LEITH
C
C **********************************************************************
C=* FROM: SPIDER - MODULAR IMAGE PROCESSING SYSTEM.   AUTHOR: J.FRANK  *
C=* Copyright (C) 1985-2005  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    RFACTSDO: DIFFERENTIAL R-FACTOR AND PHASE RESIDUAL COMPARISON OF 
C              TWO FOURIER TRANSFORMS ON A SERIES OF RINGS.  FOLLOWS
C              PHILOSOPHY OF FRANK ET AL. SCIENCE 214 (1981) 1353-1355. 
C              IN THE CURRENT ROUTINE, ALL RINGS ARE COMPUTED AT ONCE,
C              BUT SCALE SEARCH IS DONE ON EACH RING SEPARATELY.  THUS THE
C              RESULT IS EQUIVALENT TO THE RESULT OF APPLYING A SERIES
C              OF CALLS TO "RF S" WITH SUCCESSIVE RINGS.
C NOTE THAT THIS APPROACH WILL LEAD TO UNREASONABLE PHASE RESIDUAL RESULTS 
C IF THE TWO FOURIER TRANSFORMS HAVE STRONGLY DIFFERENT RADIAL BEHAVIORS.
C
C IMAGE_PROCESSING_ROUTINE
C
C23456789012345678901234567890123456789012345678901234567890123456789012
C-************************************************************************

        SUBROUTINE RFACTSDO(LUN1,LUN2,NSAM,NROW,NSLICE,
     &                               NSAM2,NROW2,NSLICE2,
     &                               ITYPE1,ITYPE2)

        INCLUDE 'CMBLOCK.INC'      
 
        REAL, ALLOCATABLE, DIMENSION(:) ::    AIMG,BIMG,CSUM1,CSUM
        REAL, ALLOCATABLE, DIMENSION(:,:) ::  PR,AMP,AVSUM
        REAL, ALLOCATABLE, DIMENSION(:) ::    CSUM2
        INTEGER, ALLOCATABLE, DIMENSION(:) :: LR

        CHARACTER*1       SER

        PARAMETER(NSCALE=20)

        CALL RDPRM(WI,NOT_USED,'RING WIDTH')

        CALL RDPRM2(SCALE1,SCALE2,NOT_USED,
     &             'SCALE FACTOR(LOWER,UPPER)')

        IF (ITYPE1 .LT. 0) THEN
C          FOURIER INPUT FILES
           LSD  = NSAM
           NSAM = NSAM - MOD(-ITYPE2,10)
        ELSE
           LSD = NSAM+2-MOD(NSAM,2)
        ENDIF

        SER    = 'C'
        SSANG  = 90.0

C       FACTOR FOR NOISE COMPARISON
        FACT   = 3.0
        DSCALE = (SCALE2-SCALE1)/FLOAT(NSCALE-1)
        Y1     = FLOAT(MAX0(NSAM,NROW))
        INC    = INT(Y1/WI)/2+1

        ALLOCATE (AIMG(LSD*NROW),BIMG(LSD*NROW), STAT=IRTFLG)
        IF (IRTFLG.NE.0) THEN 
           CALL ERRT(46,'RF, AIMG & BIMG',IER)
           GOTO 9999
        ENDIF

        IF (ITYPE1 .EQ. 1) THEN
C          INPUT FILES ARE REAL SPACE, NOT FOURIER
           DO  I = 1,NROW
              CALL REDLIN(LUN1,AIMG(1+(I-1)*LSD),NSAM,I)
           ENDDO

           DO  I = 1,NROW
              CALL REDLIN(LUN2,BIMG(1+(I-1)*LSD),NSAM,I)
           ENDDO

C          CONVERT TO FOURIER SPACE
           INV = +1
           CALL FMRS_2(AIMG,NSAM,NROW,INV)
           IF (INV.EQ.0) THEN
              CALL ERRT(38,'RF ',NE)           
              GOTO 9999
           ENDIF

           INV = +1
           CALL FMRS_2(BIMG,NSAM,NROW,INV)
           IF (INV.EQ.0)THEN
              CALL ERRT(38,'RF ',NE)
              GOTO 9999
           ENDIF
        ELSE
C          INPUT IS FOURIER ALREADY
           CALL READV(LUN1,AIMG,LSD,NROW,NSAM,NROW,NSLICE)
           CALL READV(LUN2,BIMG,LSD,NROW,NSAM,NROW,NSLICE)

        ENDIF

        ALLOCATE(PR(NSCALE,INC), AMP(NSCALE,INC),CSUM1(INC),
     &           LR(INC),CSUM(INC),CSUM2(INC), AVSUM(NSCALE,INC),
     &           STAT=IRTFLG)
        IF (IRTFLG.NE.0) THEN 
           CALL ERRT(46,'RF,  ARRAYS',IER)
           GOTO 9999
        ENDIF

        CALL PR3DB(AIMG,BIMG,PR,AMP,CSUM1,LR,CSUM,CSUM2,
     &      AVSUM,LSD,NSAM,NROW,NSLICE,DSCALE,NSCALE,SCALE1,
     &      SSANG,INC,Y1,WI,SER)

C     WRITE RESULT INTO DOC FILE AND RESULT FILE
      WRITE(NOUT,*) '2D PHASE RESIDUE AND FOURIER SHELL CORRELATION'

      WRITE(NOUT,5600) WI
5600  FORMAT(1X,'RING WIDTH = ',G12.5)

      WRITE(NOUT,5700)
5700  FORMAT(10X,'|NUMBER|','|RING RADIUS',
     &           '|DPH|','|FRC|','|FRCCRIT|','|PIXELS|')          

      WRITE(NOUT,5800)
5800  FORMAT(10X,'        ',' NORMALIZED|')             

        CALL  RFACTSD2(PR,AMP,CSUM1,LR,CSUM,CSUM2,AVSUM,
     &                 NSCALE,INC,WI,FACT,NOUT)

9999    CLOSE(LUN1)
        CLOSE(LUN2)
        IF (ALLOCATED(AIMG))  DEALLOCATE (AIMG)
        IF (ALLOCATED(BIMG))  DEALLOCATE (BIMG)
        IF (ALLOCATED(PR))    DEALLOCATE (PR)
        IF (ALLOCATED(AMP))   DEALLOCATE (AMP)
        IF (ALLOCATED(CSUM1)) DEALLOCATE (CSUM1)
        IF (ALLOCATED(LR))    DEALLOCATE (LR)
        IF (ALLOCATED(CSUM))  DEALLOCATE (CSUM)
        IF (ALLOCATED(CSUM2)) DEALLOCATE (CSUM2)
        IF (ALLOCATED(AVSUM)) DEALLOCATE (AVSUM)

        RETURN
        END