C++*********************************************************************
C
C    THRESH.F   
C                     ADDED NREPL                 OCT 2007 ARDEAN LEITH  
C    
C **********************************************************************
C=* FROM: SPIDER - MODULAR IMAGE PROCESSING SYSTEM.   AUTHOR: J.FRANK  *
C=* Copyright (C) 1985-2007  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    THRESH(LUN1,LUN2,NSAM,NROW,NSLICE)
C
C    PURPOSE: THRESHOLD AN IMAGE FILE
C
C--********************************************************************

	SUBROUTINE THRESH(LUN1,LUN2,NSAM,NROW,NSLICE)

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

        COMMON /IOBUF/ BUF(NBUFSIZ)

        CHARACTER   ANS,NULL
	REAL        B(2,2)

	DATA B/0.0,1.0,1.0,0.0/

        NREPL = 0

        IF (FCHAR(4:4) .EQ. 'M') THEN
           CALL RDPRMC(ANS,NCHAR,.TRUE.,
     &     'BLANK OUT  (A)BOVE OR (B)ELOW THRESHOLD? (A/B)',NULL,IRT)
           IF (IRT .NE. 0) RETURN

           THR = 1.0
	   CALL RDPRM1S(THR,NOT_USED,'THRESHOLD',IRT)
           IF (IRT .NE. 0) RETURN

 	   ISW = 1
	   IF (ANS .EQ. 'A') ISW = 2

           DO  J=1, NSLICE     
	      DO  I=1,NROW
	         CALL REDLIN(LUN1,BUF,NSAM,(J-1)*NROW+I)
	         DO  K=1,NSAM
	            IF (BUF(K) .LT. THR) THEN 
	               BUF(K) = B(1,ISW)
	            ELSE 
	               BUF(K) = B(2,ISW)
	            ENDIF
	         ENDDO
                 CALL WRTLIN(LUN2,BUF,NSAM,(J-1)*NROW+I)
	      ENDDO
	   ENDDO
           FREPL = FLOAT(NREPL)
           CALL REG_SET_NSEL(1,1, FREPL,0.0, 0.0, 0.0, 0.0,IRTFLG)

           WRITE(NOUT,*) ' PIXELS INSIDE MASK: ',NREPL
 	   RETURN
        ENDIF



        CALL RDPRMC(ANS,NCHAR,.TRUE.,
     &   'ALTER (A)BOVE THRESHOLD, (B)ELOW, OR (C) BOTH SIDES (A/B/C)',
     &     NULL,IRT)
        IF (IRT .NE. 0) RETURN

	IF (ANS .EQ. 'C') THEN
           CALL RDPRM2S(TH1,TH2,NOT_USED,'UPPER, LOWER THRESHOLD',IRT)
           IF (IRT .NE. 0) RETURN

        ELSE
           CALL RDPRM1S(THR,NOT_USED,'THRESHOLD',IRT)
           IF (IRT .NE. 0) RETURN

	   IF (ANS .EQ. 'B')  THEN
	      TH2 = THR
	      TH1 = HUGE(TH1)
	   ELSE
	      TH1 = THR
	      TH2 = -HUGE(TH1)
	   ENDIF
        ENDIF

 	FIX1 = TH1
	FIX2 = TH2
	IF (FCHAR(4:4) .EQ. 'F')  THEN
	   CALL RDPRM1S(FIX1,NOT_USED,'FIXUP DENSITY',IRT)
           IF (IRT .NE. 0) RETURN
	   FIX2 = FIX1
	ENDIF

        NREPL = 0

	DO  J=1,NSLICE
	   DO  I=1,NROW
	      CALL REDLIN(LUN1,BUF,NSAM,(J-1)*NROW+I)
	      DO  K=1,NSAM
	         IF (BUF(K) .GE. TH1)  THEN
	            BUF(K) = FIX1
                    NREPL  = NREPL + 1

	         ELSEIF (BUF(K) .LE. TH2)  THEN
	            BUF(K) = FIX2
                    NREPL  = NREPL + 1
	         ENDIF
	      ENDDO
	      CALL WRTLIN(LUN2,BUF,NSAM,(J-1)*NROW+I)
	   ENDDO
	ENDDO

        FREPL = FLOAT(NREPL)
        CALL REG_SET_NSEL(1,1, FREPL,0.0, 0.0, 0.0, 0.0,IRTFLG)

        WRITE(NOUT,*) ' PIXELS REPLACED: ',NREPL
 
	END