C @(#)filtfl.for	19.1 (ESO-DMD) 02/25/03 13:31:22
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 Massachusetts Ave, Cambridge, 
C MA 02139, USA.
C
C Correspondence 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
      PROGRAM FILTFL
C
C++++++++++++++++++++++++++++++++++++++++++++++++++
C
C.IDENTIFICATION
C  program    FILTFL  			version 1.00	   881010
C  D. Baade       			ESO - Garching
C  						1.10	   890518
C
C.KEYWORDS
C  filtering, noise reduction
C
C.PURPOSE
C  (Partly) fill up pixels having a flux below a user-supplied threshold with 
C  the flux above the threshold in neighbouring pixels.
C
C.ALGORITHM
C  Get the names of input + output frames from IN_A + OUT_A (both 2-D). 
C  Search for pixels with flux below threshold. Search box around each 
C  such pixel for pixels with flux above threshold. Search box around each 
C  such high-flux pixel for additional low-flux pixels. Now redistribute 
C  flux from high pixels into central low pixel in such a fashion that 
C  a) all high pixels contribute the same amount to the central low pixel,
C  b) a given high pixel gives the same flux to all its neighboring low 
C     pixels, and 
C  c) no high pixel drops below the threshold. 
C	
C.INPUT/OUTPUT
C  the following keys are used:
C
C  IN_A/C/1/60			input frame
C  OUT_A/C/1/60			output frame
C  INPUTI/I/1/2			x and y dimensions of search box
C  INPUTR/R/1/1			threshold
C
C.VERSIONS
C	
C 010201		last modif
C 
C --------------------------------------------------------------------
C
      IMPLICIT NONE
C      
      INTEGER   STAT,NAXIS,NPIX(2),IAV
      INTEGER   INNO,OUTNO
      INTEGER   RADIUS(2),THRESH
      INTEGER   UNI(1),NULO,MADRID(1)
      INTEGER*8   INPT,OUTPT
C
      CHARACTER INFRM*60,OUTFRM*60,CUNIT*64,IDENT*72
C      
      DOUBLE PRECISION  START(2),STEP(2)
C 
      INCLUDE  'MID_INCLUDE:ST_DEF.INC'
C 
      COMMON   /VMR/  MADRID
C 
      INCLUDE  'MID_INCLUDE:ST_DAT.INC'
C      
C  set up MIDAS environment + enable automatic error abort
      CALL STSPRO('FILTFL')
C
C  get name of input image  + map frame 
      CALL STKRDC('IN_A',1,1,60,IAV,INFRM,UNI,NULO,STAT)
      CALL STIGET(INFRM,D_R4_FORMAT,F_I_MODE,F_IMA_TYPE,
     +            2,NAXIS,NPIX,START,STEP,
     +            IDENT,CUNIT,INPT,INNO,STAT)
      IF (NAXIS.NE.2) 
     +   CALL STETER(1,'Currently only 2-D frames supported!')
C
C  Get name of output frame, map it for writing
      CALL STKRDC('OUT_A',1,1,60,IAV,OUTFRM,UNI,NULO,STAT)
      IF (OUTFRM.EQ.INFRM)             ! in- and output frame must be different
     +   CALL STETER
     +   (2,'Names of input + output frames must be different!')
      CALL STIPUT(OUTFRM,D_R4_FORMAT,F_O_MODE,F_IMA_TYPE,
     +            NAXIS,NPIX,START,
     +            STEP,IDENT,CUNIT,OUTPT,OUTNO,STAT)
C 
C  copy input frame to output
      IAV = NPIX(1) * NPIX(2)
      CALL COPYF(MADRID(INPT),MADRID(OUTPT),IAV)
C
C  Get dimensions of filter box and the filtering threshold
      CALL STKRDI('INPUTI',1,2,IAV,RADIUS,UNI,NULO,STAT)
      CALL STKRDR('INPUTR',1,1,IAV,THRESH,UNI,NULO,STAT)
C      
C  Now comes the real thing:
      CALL FILLPIT(MADRID(INPT),MADRID(OUTPT),NPIX(1),NPIX(2),
     +             RADIUS(1),RADIUS(2),THRESH)
C
C  Update descriptor HISTORY
      CALL DSCUPT(INNO,OUTNO,' ',STAT)
C
C  Free data
      CALL STSEPI
      END

      SUBROUTINE FILLPIT(IN,OUT,NX,NY,RX,RY,TR)
C
      IMPLICIT NONE
C      
      INTEGER    NX,NY,RX,RY,DX,DY,I,J,K,L,M,N,IK,JL,IM,JN
      INTEGER    XSTART,XEND,YSTART,YEND,NHI1,NLO2,N20,NPERC
      INTEGER    NCHECK
C      
      REAL       IN(NX,NY),OUT(NX,NY),TR,SURPLUS,TAKE,GIVE
C
      DX=2*RX+1
      DY=2*RY+1
      XSTART=DX+1
      XEND=NX-DX
      YSTART=DY+1
      YEND=NY-DY
C
C  Prepare for progress messages
C
      N20=NX/5
      IF (N20.LT.10) THEN 
         NCHECK=999999                  ! not worth the effort
      ELSE
         NCHECK=N20
         NPERC=20
      ENDIF
C      
C  First treat all pixels farther than DX (DY) from frame boundaries 
C  in X (Y). 
C
      DO 2000 I=XSTART,XEND
         DO 1800 J=YSTART,YEND
            IF (IN(I,J).LT.TR) THEN
               NHI1=0
               GIVE=1.E+12
C
C  Consider all pixels within NX x NY box centered on the pixel to be filled 
C
                DO 1600 K=1,DX
                   IK=I+K-RX-1
                   DO 1400 L=1,DY
                      JL=J+L-RY-1
                      IF (IN(IK,JL).GT.TR) THEN 
                         NHI1=NHI1+1
                         NLO2=0
C
C  For all pixels with a flux above the threshold, check if there are 
C  other low-value pixels in their respective NX x NY neighborhood boxes:
C
                         DO 1200 M=1,DX      
                            IM=IK+M-RX-1
                            DO 1000 N=1,DY
                               JN=JL+N-RY-1
                               IF (IN(IM,JN).LT.TR) NLO2=NLO2+1
1000                        CONTINUE               ! Counted all secondary 
1200                     CONTINUE                  ! low-flux pixels ...
C
C  Pixel (IK,JL) can donate at most SURPLUS. However, it may have to divide 
C  its donation potential among NLO2 low-flux pixels. Therefore, earmark only 
C  an accordingly smaller fraction for donation to current low pixel (I,J). 
C  Furthermore, no pixel shall give more than it can afford.
C
                         SURPLUS=IN(IK,JL)-TR
                         IF (NLO2.GT.1) SURPLUS=SURPLUS/NLO2
                         IF (GIVE.GT.SURPLUS) GIVE=SURPLUS        ! always give maximimum 
                      ENDIF                  ! ... around this high-flux pixel
1400               CONTINUE
1600            CONTINUE
C
C  If no high-flux pixel has been found in the box, skip to next pixel in 
C  input image.
C
                IF (NHI1.EQ.0) GOTO 1800
C
C  All flux donations available for pixel (I,J) have been summed up, now 
C  work out if they really have to give all they can. 
C
                TAKE=NHI1*GIVE
                IF ((TAKE+IN(I,J)).GT.TR) THEN 
                   TAKE=TR-IN(I,J)
                   GIVE=TAKE/NHI1
                ENDIF
C
C  Do the donation: 
C
                DO 1700 K=1,DX
                   IK=I+K-RX-1
                   DO 1660 L=1,DY
                      JL=J+L-RY-1
                      IF (IN(IK,JL).GT.TR) 
     +                   OUT(IK,JL)=OUT(IK,JL)-GIVE 
1660               CONTINUE                  ! all donations 
1700            CONTINUE                     ! made
                OUT(I,J)=OUT(I,J)+TAKE       ! and received
            ENDIF                            ! this low pixel filled
C
1800     CONTINUE                              ! all inner pixels 
C
C  Issue progress messages
C
         IF (I.GE.NCHECK) CALL PROGRS (20,N20,NPERC,NCHECK)
C
2000  CONTINUE                                     ! dealt with
C
C  Next, consider outer pixels close to frame boundaries. Begin at left:
C
      XSTART=1
      XEND=DX
C 
      DO 6000 I=XSTART,XEND
         DO 4000 J=YSTART,YEND
            IF (IN(I,J).LT.TR) THEN
               NHI1=0
               GIVE=1.E+12
               DO 3600 K=1,DX
                  IK=I+K-RX-1
                  IF (IK.LT.1) GOTO 3600
C 
                  DO 3400 L=1,DY
                     JL=J+L-RY-1
                     IF (IN(IK,JL).GT.TR) THEN 
                        NHI1=NHI1+1
                        NLO2=0
                        DO 3200 M=1,DX      
                           IM=IK+M-RX-1
                           IF (IM.LT.1) GOTO 3200
                           DO 3000 N=1,DY
                              JN=JL+N-RY-1
                              IF (IN(IM,JN).LT.TR) NLO2=NLO2+1
3000                       CONTINUE                ! Counted all secondary 
3200                    CONTINUE                   ! low-flux pixels ...
                        SURPLUS=IN(IK,JL)-TR
                        IF (NLO2.GT.1) SURPLUS=SURPLUS/NLO2
                        IF (SURPLUS.LT.GIVE) GIVE=SURPLUS
                     ENDIF                  ! ... around this high-flux pixel
3400              CONTINUE
3600           CONTINUE
C 
               IF (NHI1.EQ.0) GOTO 4000
C 
               TAKE=NHI1*GIVE
               IF ((TAKE+IN(I,J)).GT.TR) THEN 
                  TAKE=TR-IN(I,J)
                  GIVE=TAKE/NHI1
               ENDIF
C 
               DO 3800 K=1,DX
                  IK=I+K-RX-1
                  IF (IK.LT.1) GOTO 3800
                  DO 3660 L=1,DY
                     JL=J+L-RY-1
                     IF (IN(IK,JL).GT.TR) OUT(IK,JL)=OUT(IK,JL)-GIVE
3660              CONTINUE                  ! all donations 
3800           CONTINUE                        !      made
               OUT(I,J)=OUT(I,J)+TAKE                  ! and received 
            ENDIF                        ! this low pixel filled
C
4000     CONTINUE
6000  CONTINUE
C
C  Continue at right edge:
C
      XSTART=NX-DX+1
      XEND=NX
C 
      DO 8000 I=XSTART,XEND
         DO 7800 J=YSTART,YEND
            IF (IN(I,J).LT.TR) THEN
               NHI1=0
               GIVE=1.E+12
               DO 7600 K=1,DX
                  IK=I+K-RX-1
                  IF (IK.GT.NX) GOTO 7600
C 
                  DO 7400 L=1,DY
                     JL=J+L-RY-1
C 
                     IF (IN(IK,JL).GT.TR) THEN 
                        NHI1=NHI1+1
                        NLO2=0
                        DO 7200 M=1,DX      
                           IM=IK+M-RX-1
                           IF (IM.GT.NX) GOTO 7200
C 
                           DO 7000 N=1,DY
                              JN=JL+N-RY-1
                              IF (IN(IM,JN).LT.TR) NLO2=NLO2+1
7000                       CONTINUE            ! Counted all secondary 
7200                    CONTINUE                  ! low-flux pixels ... 
                        SURPLUS=IN(IK,JL)-TR
                        IF (NLO2.GT.1) SURPLUS=SURPLUS/NLO2
                        IF (SURPLUS.LT.GIVE) GIVE=SURPLUS
                     ENDIF                  ! ... around this high-flux pixel
7400              CONTINUE
7600           CONTINUE
               IF (NHI1.EQ.0) GOTO 7800
C 
               TAKE=NHI1*GIVE
               IF ((TAKE+IN(I,J)).GT.TR) THEN 
                  TAKE=TR-IN(I,J)
                  GIVE=TAKE/NHI1
               ENDIF
               DO 7700 K=1,DX
                  IK=I+K-RX-1
                  IF (IK.GT.NX) GOTO 7700
C 
                  DO 7660 L=1,DY
                     JL=J+L-RY-1
                     IF (IN(IK,JL).GT.TR) OUT(IK,JL)=OUT(IK,JL)-GIVE
7660              CONTINUE                  ! all donations 
7700           CONTINUE                        !      made
               OUT(I,J)=OUT(I,J)+TAKE            ! and received
            ENDIF                        ! this low pixel filled
C
7800     CONTINUE
8000  CONTINUE
C
C  Switch to bottom:
C
      XSTART=DX+1
      XEND=NX-DX
      YSTART=1
      YEND=DY
C 
      DO 9000 I=XSTART,XEND
         DO 8900 J=YSTART,YEND
C 
            IF (IN(I,J).LT.TR) THEN
               NHI1=0
               GIVE=1.E+12
C 
               DO 8500 K=1,DX
                  IK=I+K-RX-1
                  DO 8400 L=1,DY
                     JL=J+L-RY-1
                     IF (JL.LT.1) GOTO 8400
C 
                     IF (IN(IK,JL).GT.TR) THEN 
                        NHI1=NHI1+1
                        NLO2=0
                        DO 8300 M=1,DX      
                           IM=IK+M-RX-1
                           DO 8200 N=1,DY
                              JN=JL+N-RY-1
                              IF (JN.LT.1) GOTO 8200
                              IF (IN(IM,JN).LT.TR) NLO2=NLO2+1
8200                       CONTINUE            ! Counted all secondary 
8300                    CONTINUE                  ! low-flux pixels ... 
                        SURPLUS=IN(IK,JL)-TR
                        IF (NLO2.GT.1) SURPLUS=SURPLUS/NLO2
                        IF (SURPLUS.LT.GIVE) GIVE=SURPLUS
                     ENDIF                  ! ... around this high-flux pixel
8400              CONTINUE
C 
8500           CONTINUE
C 
               IF (NHI1.EQ.0) GOTO 8900
               TAKE=NHI1*GIVE
               IF ((TAKE+IN(I,J)).GT.TR) THEN 
                  TAKE=TR-IN(I,J)
                  GIVE=TAKE/NHI1
               ENDIF
C 
               DO 8700 K=1,DX
                  IK=I+K-RX-1
                  DO 8600 L=1,DY
                     JL=J+L-RY-1
                     IF (JL.LT.1) GOTO 8600
                     IF (IN(IK,JL).GT.TR) OUT(IK,JL)=OUT(IK,JL)-GIVE
8600              CONTINUE                       ! all donations 
8700           CONTINUE                          ! made
               OUT(I,J)=OUT(I,J)+TAKE            ! and received 
            ENDIF                                ! this low pixel filled
8900     CONTINUE
C
9000  CONTINUE
C
C  Finish at top:
C
      YSTART=NY-DY+1
      YEND=NY
C 
      DO 10000 I=XSTART,XEND
         DO 9900 J=YSTART,YEND
            IF (IN(I,J).LT.TR) THEN
               NHI1=0
               GIVE=1.E+12
               DO 9500 K=1,DX
                  IK=I+K-RX-1
                  DO 9400 L=1,DY
                     JL=J+L-RY-1
                     IF (JL.GT.NY) GOTO 9400
                     IF (IN(IK,JL).GT.TR) THEN 
                        NHI1=NHI1+1
                        NLO2=0
                        DO 9300 M=1,DX      
                           IM=IK+M-RX-1
                           DO 9200 N=1,DY
                              JN=JL+N-RY-1
                              IF (JN.GT.NY) GOTO 9200
                              IF (IN(IM,JN).LT.TR) NLO2=NLO2+1
9200                       CONTINUE            ! Counted all secondary 
9300                    CONTINUE                  ! low-flux pixels ... 
                        SURPLUS=IN(IK,JL)-TR
                        IF (NLO2.GT.1) SURPLUS=SURPLUS/NLO2
                        IF (SURPLUS.LT.GIVE) GIVE=SURPLUS
                     ENDIF                  ! ... around this high-flux pixel
9400              CONTINUE
9500           CONTINUE
C 
               IF (NHI1.EQ.0) GOTO 9900
               TAKE=NHI1*GIVE
               IF ((TAKE+IN(I,J)).GT.TR) THEN 
                  TAKE=TR-IN(I,J)
                  GIVE=TAKE/NHI1
               ENDIF
C 
               DO 9700 K=1,DX
                  IK=I+K-RX-1
                  DO 9600 L=1,DY
                     JL=J+L-RY-1
                     IF (JL.GT.NY) GOTO 9600
                     IF (IN(IK,JL).GT.TR) OUT(IK,JL)=OUT(IK,JL)-GIVE
9600              CONTINUE                       ! all donations 
9700           CONTINUE                          ! made
               OUT(I,J)=OUT(I,J)+TAKE            ! and received 
            ENDIF                                ! this low pixel filled
9900     CONTINUE
C
10000 CONTINUE
C
      RETURN
      END
