C @(#)smosub2.for	19.1 (ES0-DMD) 02/25/03 14:00:30
C===========================================================================
C Copyright (C) 1995 European Southern Observatory (ESO)
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 
C the 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
C GNU General Public License for more details.
C
C You should have received a copy of the GNU General Public 
C License along with this program; if not, write to the Free 
C Software Foundation, Inc., 675 Massachusetss Ave, Cambridge, 
C MA 02139, USA.
C
C Corresponding concerning ESO-MIDAS should be addressed as follows:
C	Internet e-mail: midas@eso.org
C	Postal address: European Southern Observatory
C			Data Management Division 
C			Karl-Schwarzschild-Strasse 2
C			D 85748 Garching bei Muenchen 
C			GERMANY
C===========================================================================
C
      SUBROUTINE AVER2A(METH,A,B,NX,NY,RADIUS,TR,ALL,V,SUM,CUTS)
C
C ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C.IDENTIFICATION
C  subroutine AVER2A  		version 1.00	861216
C  K. Banse                     ESO - Garching
C  1.10	870227	  1.20	870302		1.50	880224
C  1.60	910521
C
C.KEYWORDS
C  smoothing, neighbourhood
C
C.PURPOSE
C remove spot noise
C
C.ALGORITHM
C  a) smooth image by taking the average A(x,y)of all points in neighbourhood
C     of x,y with given radius, and replace original value I(x,y) by A(x,y),
C     if average differs more than THRESH*I(x,y) from original value
C  b) as a) but take median instead of average
C
C.INPUT/OUTPUT
C  call as   AVER2A(METH,A,B,NX,NY,RADIUS,TR,ALL,V,SUM)
C
C input par:
C  METH:	char.exp.	(1:1) = method to use: always A for average
C                        	(2:2) = A for absolute,
C					R for relative threshold,
C					F for fixed threshold interval
C  A:		R*4 array	input array
C  B:		R*4 array	output array
C  NX:		I*4		no. of pixels per line
C  NY:		I*4		no. of lines in input array
C  RADIUS:	I*4 array	"radius" of neighbourhood in x and y,
C				i.e. no. of pixels in horizontal and vertical
C				direction to include
C  TR:		R*4 (array)	threshold(s) for center point to determine
C				replacement
C  ALL:		I*4		flag for inclusion of center point in averaging
C				= 1: include center point
C				= 0: do not include center point
C  V:		R*8		working buff. of initial rowsums (RADIUS(2)*2+1)
C  SUM:		R*8		working buff. of sums (NPIX(1))
C
C.VERSIONS
C  1.00       has been built on 2.80 of old module AVERAGE
C  1.10       fix bug with index DIFX in relative threshold option
C             (it was 1 too low...)
C  1.20       add fixed threshold interval option
C  1.50       move to FORTRAN 77
C  1.60       also filter the borders (input frame larger now)
C
C --------------------------------------------------------------------
C
      IMPLICIT NONE
C
      INTEGER      NX,NY,ALL
      INTEGER      DIFX,ICOUNT
      INTEGER      INDX,BINDX
      INTEGER      N,NN,NL,NO,NOZ
      INTEGER      NP,NRL,NXMRX
      INTEGER      OFF,RX,RX1,RXM1,RY,RYNX,TOP
      INTEGER      RADIUS(2),NW1,NW2,NOSUM
C
      REAL         A(*),B(*),TR(*),CUTS(*)
      REAL         W,TEST,COUNT,FACTOR,THRESH
C
      DOUBLE PRECISION   V(*),SUM(NX)
      DOUBLE PRECISION   VNEW
C
      CHARACTER*(*)      METH
C
C  initialize
C
      BINDX = 1
      RX = RADIUS(1)
      RY = RADIUS(2)
      NW1 = RX*2 + 1                  !no. of values on row of filter
      NW2 = RY*2 + 1                  !no. of values on column of filter
      ICOUNT = NW1 * NW2
      COUNT = FLOAT(ICOUNT)          !no. of pixels in neighbourhood with center
C
      RYNX = RY * NX
      NXMRX = NX - RX
      RX1 = RX + 1
      RXM1 = RX - 1
      OFF = RYNX                        !offset for center row
C
      IF (ALL.EQ.0) THEN
         FACTOR = 1. / (COUNT-1.)
      ELSE
         FACTOR = 1. / COUNT
      ENDIF
C
C  fork according to test criterion (for speed reasons)
      IF (METH(2:2).EQ.'F') THEN
         GOTO 4000
      ELSE
         THRESH = TR(1)
      ENDIF
      IF (THRESH.GT.0.) GOTO 2000
C
C
C --------------------------
C
C  here we handle boxcar filtering without thresholding
C
C --------------------------
C
C
C  set up first row sums
      NOZ = 0                        		!offset for lower left
C
      DO 300 N=1,NW2                        	!we have NW2 row sums
         V(N) = 0.D0
         DO 200 NP=NOZ+1,NOZ+NW1
            V(N) = V(N) + A(NP)
200      CONTINUE
         NOZ = NOZ + NX                        !follow with offset...
300   CONTINUE
C
C                 ***
C  do average calculations for first "center" line
C                 ***
C
C  average in the inner part
      TOP = 1                        !init to lower left corner of matrix
      NOSUM = 1                        !point for array of average values
C
      IF (ALL.EQ.0) THEN
         DO 700 INDX=OFF+RX1,OFF+NXMRX
            DIFX = TOP                  !moving offset along line
            SUM(NOSUM)= 0.D0 
C
C  sum over neighbourhood + update the row sums on the fly
C  [ if RX = 0, there should be no out-of-range problem,
C    just spill over into the next lines  (RX = 0 and RY = 0  is not allowed!) ]
C
            DO 500 NRL=1,NW2
               SUM(NOSUM) = SUM(NOSUM) + V(NRL)
               V(NRL) = V(NRL) - A(DIFX) + A(DIFX+NW1)
               DIFX = DIFX + NX
500         CONTINUE
C
C  replace pixel in question by average value
            TEST = (SUM(NOSUM) - A(INDX)) * FACTOR
            B(BINDX) = TEST 
            BINDX = BINDX + 1
            TOP = TOP + 1
            NOSUM = NOSUM + 1
700      CONTINUE
C
      ELSE
         DO 900 INDX=OFF+RX1,OFF+NXMRX
            DIFX = TOP                  !moving offset along line
            SUM(NOSUM) = 0.D0 
C
            DO 800 NRL=1,NW2
               SUM(NOSUM) = SUM(NOSUM) + V(NRL)
               V(NRL) = V(NRL) - A(DIFX) + A(DIFX+NW1)
               DIFX = DIFX + NX
800         CONTINUE
C
            TEST = SUM(NOSUM) * FACTOR
            B(BINDX) = TEST
            BINDX = BINDX + 1
            TOP = TOP + 1
            NOSUM = NOSUM + 1
900      CONTINUE
      ENDIF
C
      OFF = OFF + NX                        !update offset of center line
C
C                 ***
C  now handle average calculations for all following lines
C                 ***
C
      NO = 0                        	!pointer to elements which are out...
      DO 1600 NL=RY+2,NY-RY
C
C  copy first row sums + set up the diff-row sum
         VNEW = 0.D0
         NN = NO - NOZ
         DO 1200 NP=NOZ+1,NOZ+NW1
            VNEW = VNEW + A(NP) - A(NP+NN)
1200     CONTINUE
         TOP = NOZ + 1                  !init to lower left corner of matrix
         NOZ = NOZ + NX                  !follow with offset...
C
C  average in the inner part
         NOSUM = 1
         DIFX = NO + 1
C
         IF (ALL.EQ.0) THEN
            DO 1300 INDX=OFF+RX1,OFF+NXMRX-1
C
C  sum over neighbourhood + update the diff-row sum on the fly
               SUM(NOSUM) = SUM(NOSUM) + VNEW
               VNEW = VNEW - A(TOP) + A(TOP+NW1)
     +                     + A(DIFX) - A(DIFX+NW1)
               TEST = (SUM(NOSUM) - A(INDX)) * FACTOR
C
C  replace pixel in question by average value
               B(BINDX) = TEST 
               BINDX = BINDX + 1
               TOP = TOP + 1
               DIFX = DIFX + 1
               NOSUM = NOSUM + 1
1300        CONTINUE
C
C  handle last pixel on line specially to avoid out_of_bounds for (TOP+NW1)
            SUM(NOSUM) = SUM(NOSUM) + VNEW
            TEST = (SUM(NOSUM) - A(OFF+NXMRX)) * FACTOR
         ELSE
            DO 1400 INDX=OFF+RX1,OFF+NXMRX-1
C
               SUM(NOSUM) = SUM(NOSUM) + VNEW
               VNEW = VNEW - A(TOP) + A(TOP+NW1)
     +                     + A(DIFX) - A(DIFX+NW1)
               TEST = SUM(NOSUM) * FACTOR
               B(BINDX) = TEST
C
C  replace pixel in question by average value
               BINDX = BINDX + 1
               TOP = TOP + 1
               DIFX = DIFX + 1
               NOSUM = NOSUM + 1
1400        CONTINUE
C
C  handle last pixel on line specially to avoid out_of_bounds for (TOP+NW1)
            SUM(NOSUM) = SUM(NOSUM) + VNEW
            TEST = SUM(NOSUM) * FACTOR
         ENDIF
C
         B(BINDX) = TEST                       !for last pixel...
         BINDX = BINDX + 1
         OFF = OFF + NX                        !update offset of center line
         NO = NO + NX
C
1600  CONTINUE
C
C finally update the cuts
      CALL CHCUTS(B,BINDX,CUTS)
      RETURN
C
C  check, if absolute or relative thresholding
C
2000  IF (METH(2:2).NE.'A') GOTO 3000
C
C
C
C --------------------------
C
C  here we handle boxcar filtering with absolute thresholding
C
C --------------------------
C
C
C
C  set up first row sums
      NOZ = 0                        		!offset for lower left
C
      DO 2100 N=1,NW2                        	!we have NW2 row sums
         V(N) = 0.D0
         DO 2060 NP=NOZ+1,NOZ+NW1
            V(N) = V(N) + A(NP)
2060     CONTINUE
         NOZ = NOZ + NX                        !follow with offset...
2100  CONTINUE
C
C                 ***
C  do average calculations for first "center" line
C                 ***
C
C  average in the inner part
      TOP = 1                        !init to lower left corner of matrix
      NOSUM = 1                        !point for array of average values
C
      IF (ALL.EQ.0) THEN
         DO 2200 INDX=OFF+RX1,OFF+NXMRX
            DIFX = TOP                  !moving offset along line
            W = A(INDX)
            SUM(NOSUM)= 0.D0 
C
            DO 2180 NRL=1,NW2
               SUM(NOSUM) = SUM(NOSUM) + V(NRL)
               V(NRL) = V(NRL) - A(DIFX) + A(DIFX+NW1)
               DIFX = DIFX + NX
2180        CONTINUE
C
C  replace pixel in question by average value
            TEST = (SUM(NOSUM) - W) * FACTOR
            IF (ABS(TEST-W).LE.THRESH) TEST = W
            B(BINDX) = TEST
            BINDX = BINDX + 1
            TOP = TOP + 1
            NOSUM = NOSUM + 1
2200     CONTINUE
C
      ELSE
         DO 2300 INDX=OFF+RX1,OFF+NXMRX
            DIFX = TOP                  !moving offset along line
            W = A(INDX)
            SUM(NOSUM) = 0.0
C
            DO 2250 NRL=1,NW2
               SUM(NOSUM) = SUM(NOSUM) + V(NRL)
               V(NRL) = V(NRL) - A(DIFX) + A(DIFX+NW1)
               DIFX = DIFX + NX
2250        CONTINUE
C
            TEST = SUM(NOSUM) * FACTOR
            IF (ABS(TEST-W).LE.THRESH) TEST = W
            B(BINDX) = TEST
            BINDX = BINDX + 1
            TOP = TOP + 1
            NOSUM = NOSUM + 1
2300     CONTINUE
      ENDIF
C
      OFF = OFF + NX                        !update offset of center line
C
C                 ***
C  now handle average calculations for all following lines
C                 ***
C
      NO = 0                        	!offset for elements which are out...
      DO 2900 NL=RY+2,NY-RY
C
C  copy first row sums + set up the diff-row sum
         VNEW = 0.D0
         NN = NO - NOZ
         DO 2450 NP=NOZ+1,NOZ+NW1
            VNEW = VNEW + A(NP) - A(NP+NN)
2450     CONTINUE
         TOP = NOZ + 1                  !init to lower left corner of matrix
         NOZ = NOZ + NX                  !follow with offset...
C
C  average in the inner part
         NOSUM = 1
         DIFX = NO + 1
C
         IF (ALL.EQ.0) THEN
            DO 2600 INDX=OFF+RX1,OFF+NXMRX-1
               W = A(INDX)
               SUM(NOSUM) = SUM(NOSUM) + VNEW
               VNEW = VNEW - A(TOP) + A(TOP+NW1)
     +                     + A(DIFX) - A(DIFX+NW1)
               TEST = (SUM(NOSUM) - W) * FACTOR
C
C  replace pixel in question by average value
               IF (ABS(TEST-W).LE.THRESH) TEST = W
               B(BINDX) = TEST
               BINDX = BINDX + 1
               TOP = TOP + 1
               NOSUM = NOSUM + 1
               DIFX = DIFX + 1
2600        CONTINUE
C
C  handle last pixel on line specially to avoid out_of_bounds for (TOP+NW1)
            W = A(OFF+NXMRX)
            SUM(NOSUM) = SUM(NOSUM) + VNEW
            TEST = (SUM(NOSUM) - W) * FACTOR
            IF (ABS(TEST-W).LE.THRESH) TEST = W
         ELSE
            DO 2700 INDX=OFF+RX1,OFF+NXMRX-1
               W = A(INDX)
               SUM(NOSUM) = SUM(NOSUM) + VNEW
               VNEW = VNEW - A(TOP) + A(TOP+NW1)
     +                     + A(DIFX) - A(DIFX+NW1)
               TEST = SUM(NOSUM) * FACTOR
C
C  replace pixel in question by average value
               IF (ABS(TEST-W).LE.THRESH) TEST = W
               B(BINDX) = TEST
               BINDX = BINDX + 1
               TOP = TOP + 1
               NOSUM = NOSUM + 1
               DIFX = DIFX + 1
2700        CONTINUE
C
C  handle last pixel on line specially to avoid out_of_bounds for (TOP+NW1)
            W = A(OFF+NXMRX)
            SUM(NOSUM) = SUM(NOSUM) + VNEW
            TEST = SUM(NOSUM) * FACTOR
            IF (ABS(TEST-W).LE.THRESH) TEST = W
         ENDIF
C
         B(BINDX) = TEST                       !for last pixel...
         BINDX = BINDX + 1
         OFF = OFF + NX                        !update offset of center line
         NO = NO + NX
C
2900  CONTINUE
C
C finally update the cuts
      CALL CHCUTS(B,BINDX,CUTS)
      RETURN
C
C
C --------------------------
C
C  here we handle boxcar filtering with relative thresholding
C
C --------------------------
C
C
C  set up first row sums
3000  NOZ = 0                        		!offset for lower left
C
      DO 3100 N=1,NW2                        	!we have NW2 row sums
         V(N) = 0.D0
         DO 3060 NP=NOZ+1,NOZ+NW1
            V(N) = V(N) + A(NP)
3060     CONTINUE
         NOZ = NOZ + NX                        !follow with offset...
3100  CONTINUE
C
C                 ***
C  do average calculations for first "center" line
C                 ***
C
C  average in the inner part
      TOP = 1                        !init to lower left corner of matrix
      NOSUM = 1                        !point for array of average values
C
      IF (ALL.EQ.0) THEN
         DO 3300 INDX=OFF+RX1,OFF+NXMRX
            W = A(INDX)
            DIFX = TOP                  !moving offset along line
            SUM(NOSUM)= 0.D0          
C
            DO 3200 NRL=1,NW2
               SUM(NOSUM) = SUM(NOSUM) + V(NRL)
               V(NRL) = V(NRL) - A(DIFX) + A(DIFX+NW1)
               DIFX = DIFX + NX
3200        CONTINUE
C
C  replace pixel in question by average value
            TEST = (SUM(NOSUM) - W) * FACTOR
            IF (ABS(TEST-W).LE.ABS(THRESH*W)) TEST = W
            B(BINDX) = TEST
            BINDX = BINDX + 1
            TOP = TOP + 1
            NOSUM = NOSUM + 1
3300     CONTINUE
C
      ELSE
         DO 3500 INDX=OFF+RX1,OFF+NXMRX
            W = A(INDX)
            DIFX = TOP                  !moving offset along line
            SUM(NOSUM) = 0.D0     
C
            DO 3400 NRL=1,NW2
               SUM(NOSUM) = SUM(NOSUM) + V(NRL)
               V(NRL) = V(NRL) - A(DIFX) + A(DIFX+NW1)
               DIFX = DIFX + NX
3400        CONTINUE
C
            TEST = SUM(NOSUM) * FACTOR
            IF (ABS(TEST-W).LE.ABS(THRESH*W)) TEST = W
            B(BINDX) = TEST
            BINDX = BINDX + 1
            TOP = TOP + 1
            NOSUM = NOSUM + 1
3500     CONTINUE
      ENDIF
C
      OFF = OFF + NX                        !update offset of center line
C
C                 ***
C  now handle average calculations for all following lines
C                 ***
C
      NO = 0                        	!offset for elements which are out...
      DO 3900 NL=RY+2,NY-RY
C
C  copy first row sums + set up the diff-row sum
         VNEW = 0.D0
         NN = NO - NOZ
         DO 3650 NP=NOZ+1,NOZ+NW1
            VNEW = VNEW + A(NP) - A(NP+NN)
3650     CONTINUE
         TOP = NOZ + 1                  !init to lower left corner of matrix
         NOZ = NOZ + NX                  !follow with offset...
C
C  average in the inner part
         NOSUM = 1
         DIFX = NO + 1
C
         IF (ALL.EQ.0) THEN
            DO 3700 INDX=OFF+RX1,OFF+NXMRX-1
               W = A(INDX)
               SUM(NOSUM) = SUM(NOSUM) + VNEW
               VNEW = VNEW - A(TOP) + A(TOP+NW1)
     +                     + A(DIFX) - A(DIFX+NW1)
               TEST = (SUM(NOSUM) - W) * FACTOR
C
C  replace pixel in question by average value
               IF (ABS(TEST-W).LE.ABS(THRESH*W)) TEST = W
               B(BINDX) = TEST
               BINDX = BINDX + 1
               TOP = TOP + 1
               NOSUM = NOSUM + 1
               DIFX = DIFX + 1
3700        CONTINUE
C
C  handle last pixel on line specially to avoid out_of_bounds for (TOP+NW1)
            W = A(OFF+NXMRX)
            SUM(NOSUM) = SUM(NOSUM) + VNEW
            TEST = (SUM(NOSUM) - W) * FACTOR
            IF (ABS(TEST-W).LE.ABS(THRESH*W)) TEST = W
         ELSE
            DO 3800 INDX=OFF+RX1,OFF+NXMRX-1
               W = A(INDX)
               SUM(NOSUM) = SUM(NOSUM) + VNEW
               VNEW = VNEW - A(TOP) + A(TOP+NW1)
     +                     + A(DIFX) - A(DIFX+NW1)
               TEST = SUM(NOSUM) * FACTOR
C
C  replace pixel in question by average value
               IF (ABS(TEST-W).LE.ABS(THRESH*W)) TEST = W
               B(BINDX) = TEST
               BINDX = BINDX + 1
               TOP = TOP + 1
               NOSUM = NOSUM + 1
               DIFX = DIFX + 1
3800        CONTINUE
C
C  handle last pixel on line specially to avoid out_of_bounds for (TOP+NW1)
            W = A(OFF+NXMRX)
            SUM(NOSUM) = SUM(NOSUM) + VNEW
            TEST = SUM(NOSUM) * FACTOR
            IF (ABS(TEST-W).LE.ABS(THRESH*W)) TEST = W
         ENDIF
C
         B(BINDX) = TEST                       !for last pixel...
         BINDX = BINDX + 1
         OFF = OFF + NX                        !update offset of center line
         NO = NO + NX
C
3900  CONTINUE
C
C finally update the cuts
      CALL CHCUTS(B,BINDX,CUTS)
      RETURN
C
C
C --------------------------
C
C  here we handle boxcar filtering with fixed thresholding interval
C
C --------------------------
C
C
C  set up first row sums
4000  NOZ = 0                        		!offset for lower left
C
      DO 4100 N=1,NW2                        	!we have NW2 row sums
         V(N) = 0.D0
         DO 4060 NP=NOZ+1,NOZ+NW1
            V(N) = V(N) + A(NP)
4060     CONTINUE
         NOZ = NOZ + NX                        !follow with offset...
4100  CONTINUE
C
C                 ***
C  do average calculations for first "center" line
C                 ***
C
C  average in the inner part
      TOP = 1                        !init to lower left corner of matrix
      NOSUM = 1                        !point for array of average values
C
      IF (ALL.EQ.0) THEN
         DO 4300 INDX=OFF+RX1,OFF+NXMRX
            TEST = A(INDX)
            DIFX = TOP                  !moving offset along line
            SUM(NOSUM)= 0.D0   
C
            DO 4200 NRL=1,NW2
               SUM(NOSUM) = SUM(NOSUM) + V(NRL)
               V(NRL) = V(NRL) - A(DIFX) + A(DIFX+NW1)
               DIFX = DIFX + NX
4200        CONTINUE
C
C  replace pixel in question by average value
            IF ( (TEST.GE.TR(1)) .AND. (TEST.LE.TR(2)) ) 
     +         TEST = (SUM(NOSUM) - TEST) * FACTOR
            B(BINDX) = TEST
            BINDX = BINDX + 1
            TOP = TOP + 1
            NOSUM = NOSUM + 1
4300     CONTINUE
C
      ELSE
         DO 4500 INDX=OFF+RX1,OFF+NXMRX
            TEST = A(INDX)
            DIFX = TOP                  !moving offset along line
            SUM(NOSUM) = 0.D0  
C
            DO 4400 NRL=1,NW2
               SUM(NOSUM) = SUM(NOSUM) + V(NRL)
               V(NRL) = V(NRL) - A(DIFX) + A(DIFX+NW1)
               DIFX = DIFX + NX
4400        CONTINUE
C
            IF ( (TEST.GE.TR(1)) .AND. (TEST.LE.TR(2)) )
     +         TEST = SUM(NOSUM) * FACTOR
            B(BINDX) = TEST
            BINDX = BINDX + 1
            TOP = TOP + 1
            NOSUM = NOSUM + 1
4500     CONTINUE
      ENDIF
C
      OFF = OFF + NX                        !update offset of center line
C
C                 ***
C  now handle average calculations for all following lines
C                 ***
C
      NO = 0                        	!offset for elements which are out...
      DO 4900 NL=RY+2,NY-RY
C
C  copy first row sums + set up the diff-row sum
         VNEW = 0.D0
         NN = NO - NOZ
         DO 4650 NP=NOZ+1,NOZ+NW1
            VNEW = VNEW + A(NP) - A(NP+NN)
4650     CONTINUE
         TOP = NOZ + 1                  !init to lower left corner of matrix
         NOZ = NOZ + NX                  !follow with offset...
C
C  average in the inner part
         NOSUM = 1
         DIFX = NO + 1
C
         IF (ALL.EQ.0) THEN
            DO 4700 INDX=OFF+RX1,OFF+NXMRX-1
               TEST = A(INDX)
               SUM(NOSUM) = SUM(NOSUM) + VNEW
               VNEW = VNEW - A(TOP) + A(TOP+NW1)
     +                     + A(DIFX) - A(DIFX+NW1)
C
C  replace pixel in question by average value
               IF ( (TEST.GE.TR(1)) .AND. (TEST.LE.TR(2)) )
     +            TEST = (SUM(NOSUM) - TEST) * FACTOR
               B(BINDX) = TEST
               BINDX = BINDX + 1
               TOP = TOP + 1
               NOSUM = NOSUM + 1
               DIFX = DIFX + 1
4700        CONTINUE
C
C  handle last pixel on line specially to avoid out_of_bounds for (TOP+NW1)
            TEST = A(OFF+NXMRX)
            SUM(NOSUM) = SUM(NOSUM) + VNEW
            IF ( (TEST.GE.TR(1)) .AND. (TEST.LE.TR(2)) )
     +         TEST = (SUM(NOSUM) - TEST) * FACTOR
         ELSE
            DO 4800 INDX=OFF+RX1,OFF+NXMRX-1
               TEST = A(INDX)
               SUM(NOSUM) = SUM(NOSUM) + VNEW
               VNEW = VNEW - A(TOP) + A(TOP+NW1)
     +                     + A(DIFX) - A(DIFX+NW1)
C
C  replace pixel in question by average value
               IF ( (TEST.GE.TR(1)) .AND. (TEST.LE.TR(2)) )
     +            TEST = SUM(NOSUM) * FACTOR
               B(BINDX) = TEST
               BINDX = BINDX + 1
               TOP = TOP + 1
               NOSUM = NOSUM + 1
               DIFX = DIFX + 1
4800        CONTINUE
C
C  handle last pixel on line specially to avoid out_of_bounds for (TOP+NW1)
            TEST = A(OFF+NXMRX)
            SUM(NOSUM) = SUM(NOSUM) + VNEW
            IF ( (TEST.GE.TR(1)) .AND. (TEST.LE.TR(2)) )
     +         TEST = SUM(NOSUM) * FACTOR
         ENDIF
C
         B(BINDX) = TEST                       !for last pixel...
         BINDX = BINDX + 1
         OFF = OFF + NX                        !update offset of center line
         NO = NO + NX
C
4900  CONTINUE
C 
C finally update the cuts
      CALL CHCUTS(B,BINDX,CUTS)
      RETURN
C 
      END

      SUBROUTINE BVER2A(METH,A,B,NX,NY,RADIUS,TR,ALL,V,SUM,CUTS)
C
      IMPLICIT NONE
C
      INTEGER      NX,NY,ALL
      INTEGER      DIFX,ICOUNT
      INTEGER      INDX,BINDX
      INTEGER      N,NN,NL,NO,NOZ
      INTEGER      NP,NRL,NXMRX
      INTEGER      OFF,RX,RX1,RXM1,RY,RYNX,TOP
      INTEGER      RADIUS(2),NW1,NW2,NOSUM
C
      REAL         A(*),B(*),TR(*),CUTS(*)
      REAL         W,TEST,COUNT,FACTOR,THRESH
C
      DOUBLE PRECISION   V(*),SUM(*)
      DOUBLE PRECISION   VNEW
C
      CHARACTER*(*)      METH
C
C  initialize
C
      BINDX = 1
      RX = RADIUS(1)
      RY = RADIUS(2)
      NW1 = RX*2 + 1                  !no. of values on row of filter
      NW2 = RY*2 + 1                  !no. of values on column of filter
      ICOUNT = NW1 * NW2
      COUNT = FLOAT(ICOUNT)          !no. of pixels in neighbourhood with center
C
      RYNX = RY * NX
      NXMRX = NX - RX
      RX1 = RX + 1
      RXM1 = RX - 1
      OFF = RYNX                        !offset for center row
C
      IF (ALL.EQ.0) THEN
         FACTOR = 1. / (COUNT-1.)
      ELSE
         FACTOR = 1. / COUNT
      ENDIF
C
C  fork according to test criterion (for speed reasons)
      IF (METH(2:2).EQ.'F') THEN
         GOTO 4000
      ELSE
         THRESH = TR(1)
      ENDIF
      IF (THRESH.GT.0.) GOTO 2000
C
C
C --------------------------
C
C  here we handle boxcar filtering without thresholding
C
C --------------------------
C
C
C  set up first row sums
      NOZ = 0                        		!offset for lower left
C
      DO 300 N=1,NW2                        	!we have NW2 row sums
         V(N) = 0.D0
         DO 200 NP=NOZ+1,NOZ+NW1
            V(N) = V(N) + A(NP)
200      CONTINUE
         NOZ = NOZ + NX                        !follow with offset...
300   CONTINUE
C
C                 ***
C  do average calculations for first "center" line
C                 ***
C
C  average in the inner part
      TOP = 1                        !init to lower left corner of matrix
      NOSUM = 1                        !point for array of average values
C
      IF (ALL.EQ.0) THEN
         DO 700 INDX=OFF+RX1,OFF+NXMRX
            DIFX = TOP                  !moving offset along line
            SUM(NOSUM)= 0.D0  
C
C  sum over neighbourhood + update the row sums on the fly
C  [ if RX = 0, there should be no out-of-range problem,
C    just spill over into the next lines  (RX = 0 and RY = 0  is not allowed!) ]
C
            DO 500 NRL=1,NW2
               SUM(NOSUM) = SUM(NOSUM) + V(NRL)
               V(NRL) = V(NRL) - A(DIFX) + A(DIFX+NW1)
               DIFX = DIFX + NX
500         CONTINUE
C
C  replace pixel in question by average value
            TEST = (SUM(NOSUM) - A(INDX)) * FACTOR
            B(BINDX) = TEST 
            BINDX = BINDX + 1
            TOP = TOP + 1
            NOSUM = NOSUM + 1
700      CONTINUE
C
      ELSE
         DO 900 INDX=OFF+RX1,OFF+NXMRX
            DIFX = TOP                  !moving offset along line
            SUM(NOSUM) = 0.D0  
C
            DO 800 NRL=1,NW2
               SUM(NOSUM) = SUM(NOSUM) + V(NRL)
               V(NRL) = V(NRL) - A(DIFX) + A(DIFX+NW1)
               DIFX = DIFX + NX
800         CONTINUE
C
            TEST = SUM(NOSUM) * FACTOR
            B(BINDX) = TEST
            BINDX = BINDX + 1
            TOP = TOP + 1
            NOSUM = NOSUM + 1
900      CONTINUE
      ENDIF
C
      OFF = OFF + NX                        !update offset of center line
C
C                 ***
C  now handle average calculations for all following lines
C                 ***
C
      NO = 0                        	!pointer to elements which are out...
      DO 1600 NL=RY+2,NY-RY
C
C  copy first row sums + set up the diff-row sum
         VNEW = 0.D0
         NN = NO - NOZ
         DO 1200 NP=NOZ+1,NOZ+NW1
            VNEW = VNEW + A(NP) - A(NP+NN)
1200     CONTINUE
         TOP = NOZ + 1                  !init to lower left corner of matrix
         NOZ = NOZ + NX                  !follow with offset...
C
C  average in the inner part
         NOSUM = 1
         DIFX = NO + 1
C
         IF (ALL.EQ.0) THEN
            DO 1300 INDX=OFF+RX1,OFF+NXMRX
C
C  sum over neighbourhood + update the diff-row sum on the fly
               SUM(NOSUM) = SUM(NOSUM) + VNEW
               VNEW = VNEW - A(TOP) + A(TOP+NW1)
     +                     + A(DIFX) - A(DIFX+NW1)
               TEST = (SUM(NOSUM) - A(INDX)) * FACTOR
C
C  replace pixel in question by average value
               B(BINDX) = TEST 
               BINDX = BINDX + 1
               TOP = TOP + 1
               DIFX = DIFX + 1
               NOSUM = NOSUM + 1
1300        CONTINUE
C
         ELSE
            DO 1400 INDX=OFF+RX1,OFF+NXMRX
C
               SUM(NOSUM) = SUM(NOSUM) + VNEW
               VNEW = VNEW - A(TOP) + A(TOP+NW1)
     +                     + A(DIFX) - A(DIFX+NW1)
               TEST = SUM(NOSUM) * FACTOR
               B(BINDX) = TEST
               BINDX = BINDX + 1
C
C  replace pixel in question by average value
               TOP = TOP + 1
               DIFX = DIFX + 1
               NOSUM = NOSUM + 1
1400        CONTINUE
         ENDIF
C
         OFF = OFF + NX                        !update offset of center line
         NO = NO + NX
C
1600  CONTINUE
C
C finally update the cuts
      CALL CHCUTS(B,BINDX,CUTS)
      RETURN
C
C
C  check, if absolute or relative thresholding
C
2000  IF (METH(2:2).NE.'A') GOTO 3000
C
C
C
C --------------------------
C
C  here we handle boxcar filtering with absolute thresholding
C
C --------------------------
C
C
C
C  set up first row sums
      NOZ = 0                        		!offset for lower left
C
      DO 2100 N=1,NW2                        	!we have NW2 row sums
         V(N) = 0.D0
         DO 2060 NP=NOZ+1,NOZ+NW1
            V(N) = V(N) + A(NP)
2060     CONTINUE
         NOZ = NOZ + NX                        !follow with offset...
2100  CONTINUE
C
C                 ***
C  do average calculations for first "center" line
C                 ***
C
C  average in the inner part
      TOP = 1                        !init to lower left corner of matrix
      NOSUM = 1                        !point for array of average values
C
      IF (ALL.EQ.0) THEN
         DO 2200 INDX=OFF+RX1,OFF+NXMRX
            DIFX = TOP                  !moving offset along line
            W = A(INDX)
            SUM(NOSUM)= 0.D0
C
            DO 2180 NRL=1,NW2
               SUM(NOSUM) = SUM(NOSUM) + V(NRL)
               V(NRL) = V(NRL) - A(DIFX) + A(DIFX+NW1)
               DIFX = DIFX + NX
2180        CONTINUE
C
C  replace pixel in question by average value
            TEST = (SUM(NOSUM) - W) * FACTOR
            IF (ABS(TEST-W).LE.THRESH) TEST = W
            B(BINDX) = TEST
            BINDX = BINDX + 1
            TOP = TOP + 1
            NOSUM = NOSUM + 1
2200     CONTINUE
C
      ELSE
         DO 2300 INDX=OFF+RX1,OFF+NXMRX
            DIFX = TOP                  !moving offset along line
            W = A(INDX)
            SUM(NOSUM) = 0.D0
C
            DO 2250 NRL=1,NW2
               SUM(NOSUM) = SUM(NOSUM) + V(NRL)
               V(NRL) = V(NRL) - A(DIFX) + A(DIFX+NW1)
               DIFX = DIFX + NX
2250        CONTINUE
C
            TEST = SUM(NOSUM) * FACTOR
            IF (ABS(TEST-W).LE.THRESH) TEST = W
            B(BINDX) = TEST
            BINDX = BINDX + 1
            TOP = TOP + 1
            NOSUM = NOSUM + 1
2300     CONTINUE
      ENDIF
C
      OFF = OFF + NX                        !update offset of center line
C
C                 ***
C  now handle average calculations for all following lines
C                 ***
C
      NO = 0                        	!offset for elements which are out...
      DO 2900 NL=RY+2,NY-RY
C
C  copy first row sums + set up the diff-row sum
         VNEW = 0.D0
         NN = NO - NOZ
         DO 2450 NP=NOZ+1,NOZ+NW1
            VNEW = VNEW + A(NP) - A(NP+NN)
2450     CONTINUE
         TOP = NOZ + 1                  !init to lower left corner of matrix
         NOZ = NOZ + NX                  !follow with offset...
C
C  average in the inner part
         NOSUM = 1
         DIFX = NO + 1
C
         IF (ALL.EQ.0) THEN
            DO 2600 INDX=OFF+RX1,OFF+NXMRX
               W = A(INDX)
               SUM(NOSUM) = SUM(NOSUM) + VNEW
               VNEW = VNEW - A(TOP) + A(TOP+NW1)
     +                     + A(DIFX) - A(DIFX+NW1)
               TEST = (SUM(NOSUM) - W) * FACTOR
C
C  replace pixel in question by average value
               IF (ABS(TEST-W).LE.THRESH) TEST = W
               B(BINDX) = TEST
               BINDX = BINDX + 1
               TOP = TOP + 1
               NOSUM = NOSUM + 1
               DIFX = DIFX + 1
2600        CONTINUE
C
         ELSE
            DO 2700 INDX=OFF+RX1,OFF+NXMRX
               W = A(INDX)
               SUM(NOSUM) = SUM(NOSUM) + VNEW
               VNEW = VNEW - A(TOP) + A(TOP+NW1)
     +                     + A(DIFX) - A(DIFX+NW1)
               TEST = SUM(NOSUM) * FACTOR
C
C  replace pixel in question by average value
               IF (ABS(TEST-W).LE.THRESH) TEST = W
               B(BINDX) = TEST
               TOP = TOP + 1
               NOSUM = NOSUM + 1
               DIFX = DIFX + 1
2700        CONTINUE
         ENDIF
C
         OFF = OFF + NX                        !update offset of center line
         NO = NO + NX
C
2900  CONTINUE
C
C finally update the cuts
      CALL CHCUTS(B,BINDX,CUTS)
      RETURN
C
C
C --------------------------
C
C  here we handle boxcar filtering with relative thresholding
C
C --------------------------
C
C
C  set up first row sums
3000  NOZ = 0                        		!offset for lower left
C
      DO 3100 N=1,NW2                        	!we have NW2 row sums
         V(N) = 0.D0
         DO 3060 NP=NOZ+1,NOZ+NW1
            V(N) = V(N) + A(NP)
3060     CONTINUE
         NOZ = NOZ + NX                        !follow with offset...
3100  CONTINUE
C
C                 ***
C  do average calculations for first "center" line
C                 ***
C
C  average in the inner part
      TOP = 1                        !init to lower left corner of matrix
      NOSUM = 1                        !point for array of average values
C
      IF (ALL.EQ.0) THEN
         DO 3300 INDX=OFF+RX1,OFF+NXMRX
            W = A(INDX)
            DIFX = TOP                  !moving offset along line
            SUM(NOSUM)= 0.D0       
C
            DO 3200 NRL=1,NW2
               SUM(NOSUM) = SUM(NOSUM) + V(NRL)
               V(NRL) = V(NRL) - A(DIFX) + A(DIFX+NW1)
               DIFX = DIFX + NX
3200        CONTINUE
C
C  replace pixel in question by average value
            TEST = (SUM(NOSUM) - W) * FACTOR
            IF (ABS(TEST-W).LE.ABS(THRESH*TEST)) TEST = W
            B(BINDX) = TEST
            BINDX = BINDX + 1
            TOP = TOP + 1
            NOSUM = NOSUM + 1
3300     CONTINUE
C
      ELSE
         DO 3500 INDX=OFF+RX1,OFF+NXMRX
            W = A(INDX)
            DIFX = TOP                  !moving offset along line
            SUM(NOSUM) = 0.D0
C
            DO 3400 NRL=1,NW2
               SUM(NOSUM) = SUM(NOSUM) + V(NRL)
               V(NRL) = V(NRL) - A(DIFX) + A(DIFX+NW1)
               DIFX = DIFX + NX
3400        CONTINUE
C
            TEST = SUM(NOSUM) * FACTOR
            IF (ABS(TEST-W).LE.ABS(THRESH*TEST)) TEST = W
            B(BINDX) = TEST
            BINDX = BINDX + 1
            TOP = TOP + 1
            NOSUM = NOSUM + 1
3500     CONTINUE
      ENDIF
C
      OFF = OFF + NX                        !update offset of center line
C
C                 ***
C  now handle average calculations for all following lines
C                 ***
C
      NO = 0                        	!offset for elements which are out...
      DO 3900 NL=RY+2,NY-RY
C
C  copy first row sums + set up the diff-row sum
         VNEW = 0.D0
         NN = NO - NOZ
         DO 3650 NP=NOZ+1,NOZ+NW1
            VNEW = VNEW + A(NP) - A(NP+NN)
3650     CONTINUE
         TOP = NOZ + 1                  !init to lower left corner of matrix
         NOZ = NOZ + NX                  !follow with offset...
C
C  average in the inner part
         NOSUM = 1
         DIFX = NO + 1
C
         IF (ALL.EQ.0) THEN
            DO 3700 INDX=OFF+RX1,OFF+NXMRX
               W = A(INDX)
               SUM(NOSUM) = SUM(NOSUM) + VNEW
               VNEW = VNEW - A(TOP) + A(TOP+NW1)
     +                     + A(DIFX) - A(DIFX+NW1)
               TEST = (SUM(NOSUM) - W) * FACTOR
C
C  replace pixel in question by average value
               IF (ABS(TEST-W).LE.ABS(THRESH*TEST)) TEST = W
               B(BINDX) = TEST
               BINDX = BINDX + 1
               TOP = TOP + 1
               NOSUM = NOSUM + 1
               DIFX = DIFX + 1
3700        CONTINUE
C
         ELSE
            DO 3800 INDX=OFF+RX1,OFF+NXMRX
               W = A(INDX)
               SUM(NOSUM) = SUM(NOSUM) + VNEW
               VNEW = VNEW - A(TOP) + A(TOP+NW1)
     +                     + A(DIFX) - A(DIFX+NW1)
               TEST = SUM(NOSUM) * FACTOR
C
C  replace pixel in question by average value
               IF (ABS(TEST-W).LE.ABS(THRESH*TEST)) TEST = W
               B(BINDX) = TEST
               BINDX = BINDX + 1
               TOP = TOP + 1
               NOSUM = NOSUM + 1
               DIFX = DIFX + 1
3800        CONTINUE
         ENDIF
C
         OFF = OFF + NX                        !update offset of center line
         NO = NO + NX
C
3900  CONTINUE
C
C finally update the cuts
      CALL CHCUTS(B,BINDX,CUTS)
      RETURN
C
C
C
C --------------------------
C
C  here we handle boxcar filtering with fixed thresholding interval
C
C --------------------------
C
C
C  set up first row sums
4000  NOZ = 0                        		!offset for lower left
C
      DO 4100 N=1,NW2                        	!we have NW2 row sums
         V(N) = 0.D0
         DO 4060 NP=NOZ+1,NOZ+NW1
            V(N) = V(N) + A(NP)
4060     CONTINUE
         NOZ = NOZ + NX                        !follow with offset...
4100  CONTINUE
C
C                 ***
C  do average calculations for first "center" line
C                 ***
C
C  average in the inner part
      TOP = 1                        !init to lower left corner of matrix
      NOSUM = 1                        !point for array of average values
C
      IF (ALL.EQ.0) THEN
         DO 4300 INDX=OFF+RX1,OFF+NXMRX
            TEST = A(INDX)
            DIFX = TOP                  !moving offset along line
            SUM(NOSUM)= 0.D0  
C
            DO 4200 NRL=1,NW2
               SUM(NOSUM) = SUM(NOSUM) + V(NRL)
               V(NRL) = V(NRL) - A(DIFX) + A(DIFX+NW1)
               DIFX = DIFX + NX
4200        CONTINUE
C
C  replace pixel in question by average value
            IF ( (TEST.GE.TR(1)) .AND. (TEST.LE.TR(2)) ) 
     +         TEST = (SUM(NOSUM) - TEST) * FACTOR
            B(BINDX) = TEST
            TOP = TOP + 1
            NOSUM = NOSUM + 1
4300     CONTINUE
C
      ELSE
         DO 4500 INDX=OFF+RX1,OFF+NXMRX
            TEST = A(INDX)
            DIFX = TOP                  !moving offset along line
            SUM(NOSUM) = 0.D0
C
            DO 4400 NRL=1,NW2
               SUM(NOSUM) = SUM(NOSUM) + V(NRL)
               V(NRL) = V(NRL) - A(DIFX) + A(DIFX+NW1)
               DIFX = DIFX + NX
4400        CONTINUE
C
            IF ( (TEST.GE.TR(1)) .AND. (TEST.LE.TR(2)) )
     +         TEST = SUM(NOSUM) * FACTOR
            B(BINDX) = TEST
            BINDX = BINDX + 1
            TOP = TOP + 1
            NOSUM = NOSUM + 1
4500     CONTINUE
      ENDIF
C
      OFF = OFF + NX                        !update offset of center line
C
C                 ***
C  now handle average calculations for all following lines
C                 ***
C
      NO = 0                        	!offset for elements which are out...
      DO 4900 NL=RY+2,NY-RY
C
C  copy first row sums + set up the diff-row sum
         VNEW = 0.D0
         NN = NO - NOZ
         DO 4650 NP=NOZ+1,NOZ+NW1
            VNEW = VNEW + A(NP) - A(NP+NN)
4650     CONTINUE
         TOP = NOZ + 1                  !init to lower left corner of matrix
         NOZ = NOZ + NX                  !follow with offset...
C
C  average in the inner part
         NOSUM = 1
         DIFX = NO + 1
C
         IF (ALL.EQ.0) THEN
            DO 4700 INDX=OFF+RX1,OFF+NXMRX
               TEST = A(INDX)
               SUM(NOSUM) = SUM(NOSUM) + VNEW
               VNEW = VNEW - A(TOP) + A(TOP+NW1)
     +                     + A(DIFX) - A(DIFX+NW1)
C
C  replace pixel in question by average value
               IF ( (TEST.GE.TR(1)) .AND. (TEST.LE.TR(2)) )
     +            TEST = (SUM(NOSUM) - TEST) * FACTOR
               B(BINDX) = TEST
               BINDX = BINDX + 1
               TOP = TOP + 1
               NOSUM = NOSUM + 1
               DIFX = DIFX + 1
4700        CONTINUE
C
         ELSE
            DO 4800 INDX=OFF+RX1,OFF+NXMRX
               TEST = A(INDX)
               SUM(NOSUM) = SUM(NOSUM) + VNEW
               VNEW = VNEW - A(TOP) + A(TOP+NW1)
     +                     + A(DIFX) - A(DIFX+NW1)
C
C  replace pixel in question by average value
               IF ( (TEST.GE.TR(1)) .AND. (TEST.LE.TR(2)) )
     +            TEST= SUM(NOSUM) * FACTOR
               B(BINDX) = TEST
               BINDX = BINDX + 1
               TOP = TOP + 1
               NOSUM = NOSUM + 1
               DIFX = DIFX + 1
4800        CONTINUE
         ENDIF
C
         OFF = OFF + NX                        !update offset of center line
         NO = NO + NX
C
4900  CONTINUE
C
C finally update the cuts
      CALL CHCUTS(B,BINDX,CUTS)
      RETURN
C
      END
