C++*********************************************************************
C
C FQ_Q.F 
C 12/22/94  
C                RDPRAF REMOVED                    DEC 2005 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 FQ_Q.F(LUN,LUNO,B,LSD,N2S,N2R,NSAM,NROW,IOPT)
C
C PURPOSE: QUICK FILTERING OF REAL-SPACE FILE BY IN-CORE FFT
C
C PARAMETERS:
C        LUN         LOGICAL UNIT NUMBER OF REAL-SPACE FILE TO BE FILTERED
C        LUNO        LOGICAL UNIT NUMBER OF REAL-SPACE OUTPUT FILE 
C        IOPT        TYPE OF THE FILTER
C        B           BUFFER 
C        NSAM,NROW   DIMENSIONS OF REAL-SPACE FILE
C        N2S=2*NSAM  AT LEAST
C        N2R=2*NROW  "    "
C
C IMAGE_PROCESSING_ROUTINE
C
C23456789012345678901234567890123456789012345678901234567890123456789012  
C--*******************************************************************

        SUBROUTINE FQ_Q(LUN,LUNO,B,LSD,N2S,N2R,NSAM,NROW,IOPT)
        
	INCLUDE 'CMBLOCK.INC'
	
	DIMENSION         B(LSD,N2R)
	DIMENSION         BFPS(4)
	DOUBLE PRECISION  AVE
	
	
C       To set them to something.
	PARM1 = 0.0
	PARM2 = 0.0

C       READ  IMAGE
	DO I=1,NROW
 	   CALL  REDLIN(LUN,B(1,I),NSAM,I)
	ENDDO

C       BORDER PADDING
	IF (N2S.NE.NSAM .AND. N2R.NE.NROW)  THEN
	   AVE = (SUM(B(1:NSAM,1))+SUM(B(1:NSAM,NROW))
     &	       +SUM(B(1,2:NROW-1))+SUM(B(NSAM,2:NROW-1)) )
     &		/REAL(2*(NSAM+NROW)-4)
c$omp      parallel do private(i,j)
	   DO J=1,N2R
	      DO I=NSAM+1,N2S
	         B(I,J)=AVE
	      ENDDO
	   ENDDO
c$omp      parallel do private(i,j)
	   DO J=NROW+1,N2R
	      DO I=1,NSAM
	         B(I,J)=AVE
	      ENDDO
	   ENDDO
	ENDIF

	INV=1
	CALL FMRS_2(B,N2S,N2R,INV)
	IF (INV.EQ.0) THEN
	   IOPT = -1
	   RETURN
	ENDIF

C       BUTTERWORTH FILTER ***********************************

	IF (IOPT.EQ.7 .OR. IOPT.EQ.8)  THEN
	   NMAX = 4
	   BFPS = 0.0
	   CALL RDPRA('PASS-BAND FREQUENCY & STOP-BAND FREQUENCY',NMAX,
     &        0,.FALSE.,BFPS,NGOT,IRTFLG)

	   EPS = 0.882
	   AA  = 10.624
	   ORD = 2.*ALOG10(EPS/SQRT(AA**2-1.0))
	   IF (BFPS(3).EQ.0.0 .AND. BFPS(4).EQ.0.0) THEN
	      ORD   = ORD/ALOG10(BFPS(1)/BFPS(2))
	      PARM1 = BFPS(1)/(EPS)**(2./ORD)
	   ELSE
C             BUTTERWORTH FILTER ELLIPTIC FILTER:
C             LOW-PASS  IOPT=9,  HIGH-PASS IOPT=10
	      IOPT = IOPT + 2
           ENDIF

	ELSE
  	   CALL RDPRM2(PARM1,PARM2,NOT_USED,'FILTER RADIUS')
	   IF (PARM1.LT.0.0.OR.PARM1.GT.0.5)  PARM1=0.5*PARM1/(NSAM/2)
	   IF (PARM2.EQ.0.0)  PARM2=PARM1
	   IF (PARM2.LT.0.0.OR.PARM2.GT.0.5)  PARM1=0.5*PARM2/(NROW/2)
	   IF (IOPT.EQ.5.OR.IOPT.EQ.6)  THEN

C             FERMI DISTRIBUTION FILTER ********************
	      CALL RDPRM(TEMP,NOT_USED,'TEMPERATURE(0=CUTOFF)')

C             EXPONENTIAL FOR HIGH-PASS OPTION
	      IF (IOPT.EQ.6) TEMP=-TEMP
	   ENDIF
	ENDIF

	NR2    = N2R/2
	X1     = FLOAT(N2S/2)**2
	Y1     = FLOAT(NR2)**2
	PARM   = PARM1**2
	PARM22 = PARM2**2
C       KEEP ZERO TERM FOR HIGH PASS OPTIONS
	AVG = B(1,1)

c$omp   parallel do private(i,j,ix,iy,f,fpe,fse,ordt,parmt)
	DO J=1,N2R
	   IY = (J-1)
	   IF (IY.GT.NR2) IY=IY-N2R
	   DO I=1,LSD,2
	      IX = (I-1)/2
	      IF (IOPT.EQ.1) THEN
C                LOWPASS *************************************
                 IF(
     &  0.25*(FLOAT(IX*IX)/X1/PARM+FLOAT(IY*IY)/Y1/PARM22).GT.1.0) THEN
	            B(I,J)   = 0.0
	            B(I+1,J) = 0.0
	         ENDIF
	      ELSEIF (IOPT.EQ.2) THEN	
C                HIGH PASS ******************
	         IF( (IX.NE.0 .OR. IY.NE.0) .AND.
     &  0.25*(FLOAT(IX*IX)/X1/PARM+FLOAT(IY*IY)/Y1/PARM22).LE.1.0) THEN
	            B(I,J)=0.0
	            B(I+1,J)=0.0
	         ENDIF
	      ELSEIF(IOPT.EQ.3)  THEN
C                GAUSSIAN LOW PASS ***************************
	         F=0.125*(FLOAT(IX*IX)/X1/PARM+FLOAT(IY*IY)/Y1/PARM22)
	         IF (F.LT.16.0)  THEN
	            F        = EXP(-F)
                    B(I,J)   = B(I,J)*F
                    B(I+1,J) = B(I+1,J)*F
	         ELSE
                    B(I,J)   = 0.0
                    B(I+1,J) = 0.0
	         ENDIF
	      ELSEIF (IOPT.EQ.4)  THEN	
C                GAUSSIAN HIGH PASS ********************************
	         IF (IX.NE.0 .OR. IY.NE.0)  THEN
	  F=0.125*(FLOAT(IX*IX)/X1/PARM+FLOAT(IY*IY)/Y1/PARM22)
	            IF(F.LT.16.0)  THEN
	               F=1.0-EXP(-F)
                       B(I,J)=B(I,J)*F
                       B(I+1,J)=B(I+1,J)*F
	            ENDIF
	         ENDIF

	      ELSEIF (IOPT.EQ.5.OR.IOPT.EQ.6)  THEN
C                FERMI DISTRIBUTION FILTER *************************	      
	      F=(0.5*SQRT(FLOAT(IX*IX)/X1+FLOAT(IY*IY)/Y1)-PARM1)/TEMP
	         F=AMIN1(AMAX1(F,-10.0),10.0)
                 F=(1.0/(1.0+EXP(F)))
                 B(I,J)=B(I,J)*F
                 B(I+1,J)=B(I+1,J)*F

	      ELSEIF (IOPT.EQ.7)THEN
C                BUTTERWORTH LOWPASS FILTER ************************
	         F=0.5*SQRT(FLOAT(IX*IX)/X1+FLOAT(IY*IY)/Y1)
	         F=SQRT(1.0/(1.0+(F/PARM1)**ORD))
                 B(I,J)=B(I,J)*F
                 B(I+1,J)=B(I+1,J)*F

	      ELSEIF (IOPT.EQ.8)THEN
C                BUTTERWORTH HIGHPASS FILTER ***********************	
                 IF(IX.NE.0 .OR. IY.NE.0) THEN
	            F=0.5*SQRT(FLOAT(IX*IX)/X1+FLOAT(IY*IY)/Y1)
	            F=(1.0-SQRT(1.0/(1.0+(F/PARM1)**ORD)))
                    B(I,J)=B(I,J)*F
                    B(I+1,J)=B(I+1,J)*F
	         ENDIF

	      ELSEIF (IOPT.EQ.9)THEN
C                BUTTERWORTH ELLIPTIC LOWPASS FILTER **************	
C  CALCULATE EFFECTIVE FP AND FS IN A GIVEN DIRECTION ON THE PLANE
                 IF(IX.NE.0 .OR. IY.NE.0) THEN
	            FPE=ATAN2(BFPS(1)*SQRT(FLOAT(IY*IY)/Y1),
     &                  BFPS(3)*SQRT(FLOAT(IX*IX)/X1))
          FPE=SQRT((BFPS(1)*COS(FPE))**2+(BFPS(3)*SIN(FPE))**2)
	            FSE=ATAN2(BFPS(2)*SQRT(FLOAT(IY*IY)/Y1),
     &                  BFPS(4)*SQRT(FLOAT(IX*IX)/X1))
          FSE=SQRT((BFPS(2)*COS(FSE))**2+(BFPS(4)*SIN(FSE))**2)
	            ORDT=ORD/ALOG10(FPE/FSE)
	            PARMT=FPE/(EPS)**(2./ORDT)
	            F=0.5*SQRT(FLOAT(IX*IX)/X1+FLOAT(IY*IY)/Y1)
	            F=SQRT(1.0/(1.0+(F/PARMT)**ORDT))
                    B(I,J)=B(I,J)*F
                    B(I+1,J)=B(I+1,J)*F
	         ENDIF

	      ELSEIF(IOPT.EQ.10)THEN
C                BUTTERWORTH ELLIPTIC HIGHPASS FILTER *************	
                 IF (IX.NE.0 .OR. IY.NE.0) THEN
	            FPE=ATAN2(BFPS(1)*SQRT(FLOAT(IY*IY)/Y1),
     &        BFPS(3)*SQRT(FLOAT(IX*IX)/X1))
              FPE=SQRT((BFPS(1)*COS(FPE))**2+(BFPS(3)*SIN(FPE))**2)
	            FSE=ATAN2(BFPS(2)*SQRT(FLOAT(IY*IY)/Y1),
     &                  BFPS(4)*SQRT(FLOAT(IX*IX)/X1))
          FSE=SQRT((BFPS(2)*COS(FSE))**2+(BFPS(4)*SIN(FSE))**2)
	            ORDT=ORD/ALOG10(FPE/FSE)
	            PARMT=FPE/(EPS)**(2./ORDT)
	            F=0.5*SQRT(FLOAT(IX*IX)/X1+FLOAT(IY*IY)/Y1)
	            F=(1.0-SQRT(1.0/(1.0+(F/PARMT)**ORDT)))
                    B(I,J)=B(I,J)*F
                    B(I+1,J)=B(I+1,J)*F
	         ENDIF
	      ENDIF
	   ENDDO
	ENDDO

C       RESTORE ZERO TERM FOR HIGH PASS OPTIONS
	IF(IOPT.EQ.2.OR.IOPT.EQ.4.OR.IOPT.EQ.6.OR.IOPT.EQ.8)B(1,1)=AVG

C       WRITE  IMAGE
	INV = -1
	CALL FMRS_2(B,N2S,N2R,INV)
	DO I=1,NROW
 	   CALL  WRTLIN(LUNO,B(1,I),NSAM,I)
	ENDDO

        END