C++*********************************************************************
C
C FILTDOC.F     ALLOCATION                        MAY 2000 ARDEAN LEITH
C               OPENDOC PARAMETERS                DEC 2000 ARDEAN LEITH
C               LUNDOCREDSEQ RETURNS MAXY         APR 2003 ARDEAN LEITH
C               INCORE OPENDOC                    JUL 2003 ARDEAN LEITH
C               FOURIER INPUT AND OUTPUT FILE     OCT 2003 BIMAL RATH
C               NGOTY BUG                         FEB 2004 ARDEAN LEITH
C **********************************************************************
C=*                                                                    *
C=* This file is part of:   SPIDER - Modular Image Processing System.  *
C=* SPIDER System Authors:  Joachim Frank & ArDean Leith               *
C=* Copyright 1985-2010  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=* You should have received a copy of the GNU General Public License  *
C=* along with this program. If not, see <http://www.gnu.org/licenses> *
C=*                                                                    *
C **********************************************************************
C
C  FILTDOC(LUNINT,LUNOUT,NSAM,NROW,NSLICE)
C
C  PURPOSE:      FILTERS AN IMAGE USING COEEFICIENTS FROM A
C                DOC. FILE
C
C  PARAMETERS:   LUNINT         INPUT UNIT FOR IMAGE & DOC       (SENT)
C                LUNOUT         OUTPUT UNIT FOR IMAGE            (SENT)
C                NSAM,NROW & NSLICE                              (SENT)
C
C--********************************************************************

	SUBROUTINE FILTDOC(LUNINT,LUNOUT,NSAM,NROW,NSLICE,IFORMIN)

        INCLUDE 'CMBLOCK.INC' 
        INCLUDE 'CMLIMIT.INC' 

        COMPLEX, ALLOCATABLE, DIMENSION(:,:,:) :: G
        REAL, ALLOCATABLE, DIMENSION(:)        :: COEFFS
        CHARACTER(LEN=MAXNAM)                  :: DOCFIL
        LOGICAL                                :: LDUM
  
C       SET VALUE OF ISPACE AS PER REAL OR FT INPUT IMAGE
        IF (IFORMIN .GT. 0) THEN
	   ISPACE = (NSAM + 2 - MOD(NSAM,2))/2
        ELSE
           ISPACE = NSAM
        ENDIF

	NSLICE=MAX0(1,NSLICE)
        ALLOCATE (G(ISPACE,NROW,NSLICE), STAT=IRTFLG)
        IF (IRTFLG .NE. 0) THEN
           CALL ERRT(46,'FILTDOC, G',NDUM)
           RETURN
        ENDIF

C       LOAD IMAGE/VOLUME DATA (CAN NOT USE REDVOL HERE)
        DO K=1,NSLICE
           DO J=1,NROW
              NR = J+(K-1)*NROW
              CALL REDLIN(LUNINT,G(1,J,K),NSAM,NR)
	   ENDDO
        ENDDO
        CLOSE(LUNINT)

C       FOR 2D   (SQRT(2)/2) =~ .75
	LBS  = 0.75 * MAX0(ISPACE,NROW,NSLICE)
C       FOR 3D   (SQRT(3)/2) =~0.9
        IF (NSLICE .GT. 1) LBS  = 0.9 * MAX0(ISPACE,NROW,NSLICE)
        MAXX = 1
        MAXY = LBS

        ALLOCATE (COEFFS(MAXY), STAT=IRTFLG)
        IF (IRTFLG .NE. 0) THEN
           CALL ERRT(46,'FILTDOC, COEFFS',NDUM)
           IF (ALLOCATED(G)) DEALLOCATE(G)
           RETURN
        ENDIF

         CALL OPENDOC(DOCFIL,.TRUE.,NLET,LUNINT,LUNIN,.TRUE.,'DOCUMENT',
     &            .TRUE.,.FALSE.,.TRUE.,LDUM,IRTFLG)
        IF (IRTFLG .NE. 0) GOTO 9999

C       LOAD COEFFS DATA USING SEQUENTIAL READ FROM FIRST REGISTER
        CALL LUNDOCREDSEQ(LUNIN,COEFFS,MAXX,MAXY,MAXY,NGOTY,IRTFLG)
        CLOSE(LUNINT)
        IF (IRTFLG .EQ. 1) GOTO 9999

C       ZERO REMAINING COEFFS
	IF (NGOTY .LT. LBS) COEFFS(NGOTY+1:LBS) = 0.0

C       DO FT FOR REAL INPUT IMAGES
        IF (IFORMIN .GT. 0) THEN
 	   INV = +1
 	   IF (NSLICE .LE. 1)  THEN
              CALL FMRS_2(G,NSAM,NROW,INV)
 	   ELSE
              CALL FMRS_3(G,NSAM,NROW,NSLICE,INV)
	   ENDIF
        ENDIF

c$omp   parallel do private(k,j,i,iz,iy,ir,rr,dd,fr)
	DO K=1,NSLICE
           IZ = K-1
           IF (IZ .GT. NSLICE/2)  IZ = IZ - NSLICE
	   DO J=1,NROW
              IY = J-1
              IF (IY .GT. NROW/2) IY = IY - NROW
	      DO I=1,ISPACE
	         RR       = SQRT(FLOAT((I-1)**2 + IY*IY + IZ*IZ))
	         IR       = IFIX(RR)
	         DD       = RR - FLOAT(IR)
	         FR       = (1.0-DD)*COEFFS(IR+1)+DD*COEFFS(IR+2)
	         G(I,J,K) = G(I,J,K) * FR	
	      ENDDO
	   ENDDO
	ENDDO

C       DO INVERSE FT FOR REAL INPUT IMAGES
        IF (IFORMIN .GT. 0) THEN
	   INV = -1
 	   IF (NSLICE.LE.1)  THEN
              CALL FMRS_2(G,NSAM,NROW,INV)
 	   ELSE
              CALL FMRS_3(G,NSAM,NROW,NSLICE,INV)
	   ENDIF
        ENDIF

        DO K=1,NSLICE
          DO J=1,NROW
             NR = J+(K-1)*NROW
             CALL WRTLIN(LUNOUT,G(1,J,K),NSAM,NR)
	  ENDDO
	ENDDO

C       CLOSE OUTPUT FILE
9999    CLOSE(LUNOUT)

C       DEALLOCATE RUN-TIME MEMORY
        IF (ALLOCATED(COEFFS)) DEALLOCATE(COEFFS)
        IF (ALLOCATED(G))      DEALLOCATE(G)

        RETURN
	END