C @(#)mask.for	17.1.1.1 (ESO-DMD) 01/25/02 17:52:45
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 BADPIX
C++++++++++++++++++++++++++++++++++++++++++++++++++++
C.COPYRIGHT: Copyright (c) 1992 European Southern Observatory,
C                                         all rights reserved
C.IDENTIFICATION: mask.for
C.PURPOSE:   Replace bad pixels by nearest good pixel according to pixel mask
C.LANGUAGE:  F77+ESOext
C.KEYWORDS:  IRAC2, infrared
C.AUTHOR:    Gert Finger, ESO-VLT
C.VERSION: 
C 
C 010710	last modif
C 
C---------------------------------------------------------------------
C 
      IMPLICIT NONE
C 
      INTEGER      NAXISA,NPIXA(2),IAV,STAT
      INTEGER      IMNOA,IMNOC,IMNOG
      INTEGER      KNULL,KUNIT(1)
      INTEGER      MADRID(1)
      INTEGER      FNO,IBUF(5)
C 
      INTEGER*8    PNTRA,PNTRC,PNTRG
C 
      CHARACTER*60 FRAMEA,FRAMEC,FRAMEG
      CHARACTER    CUNITA*64,IDENTA*72
C 
      DOUBLE PRECISION STEPA(2),STARTA(2)
C
      INCLUDE 'MID_INCLUDE:ST_DEF.INC'   
      COMMON /VMR/ MADRID
      INCLUDE 'MID_INCLUDE:ST_DAT.INC'   
C 
C 
C 
      CALL STSPRO('BADPIX')
C
C  GET NAME OF INPUT FRAME AND MAP IT 
C
      CALL STKRDC('IN_A',1,1,60,IAV,FRAMEA,KUNIT,KNULL,STAT)
C
C  TEST DATA FORMAT (NUMBER OF BYTES PER PIXEL = IBUF(1)
C
      FNO=2 
      CALL STFINF(FRAMEA,FNO,IBUF,STAT)
C     WRITE (6,*) (IBUF(I),I=1,5)
C
C  DOUBLE PRECISSION
C
      IF (IBUF(1).EQ.8) THEN 
      CALL STIGET(FRAMEA,D_R8_FORMAT,F_I_MODE,F_IMA_TYPE,
     +            2,NAXISA,NPIXA,STARTA,STEPA,IDENTA,
     +            CUNITA,PNTRA,IMNOA,STAT)
      GOTO 1000
      ENDIF 
C
C  SINGLE PRECISSION
C
      IF (IBUF(1).EQ.4) THEN 
      CALL STIGET(FRAMEA,D_R4_FORMAT,F_I_MODE,F_IMA_TYPE,
     +            2,NAXISA,NPIXA,STARTA,STEPA,IDENTA,
     +            CUNITA,PNTRA,IMNOA,STAT)
      GOTO 1000
      ENDIF 
C
1000  CONTINUE
      IF (NAXISA.NE.2)
     +CALL STETER(1,'INPUT FRAME MUST HAVE TWO DIMENSIONS')
C
C  GET NAME OF OUTPUT FRAME CREATE AND MAP IT
C
      CALL STKRDC('OUT_A',1,1,60,IAV,FRAMEC,KUNIT,KNULL,STAT) 
      CALL STIPUT(FRAMEC,D_R8_FORMAT,F_O_MODE,F_IMA_TYPE,
     +            NAXISA,NPIXA,STARTA,STEPA,IDENTA,
     +            CUNITA,PNTRC,IMNOC,STAT)
C
C  GET IMAGE GOODPX AND MAP IT ( IMAGE OF GOOD PIXELS )
C
      FRAMEG='GOODPX'
      CALL STIGET(FRAMEG,D_R8_FORMAT,F_I_MODE,F_IMA_TYPE,
     +            2,NAXISA,NPIXA,STARTA,STEPA,IDENTA,
     +            CUNITA,PNTRG,IMNOG,STAT)
C
C  CALL SUBROUTINE TO DO THE REPLACEMENT OF BAD PIXELS
C
C
C  DOUBLE PRECISSION
C
      IF (IBUF(1).EQ.8) THEN 
      CALL RPLBP(MADRID(PNTRA),MADRID(PNTRC),MADRID(PNTRG),
     +                NPIXA(1),NPIXA(2))
      ENDIF
C
C  SINGLE PRECISSION
C
      IF (IBUF(1).EQ.4) THEN 
      CALL RPLBPS(MADRID(PNTRA),MADRID(PNTRC),MADRID(PNTRG),
     +                NPIXA(1),NPIXA(2))
      ENDIF
C
C  COPY ALL DESCRIPTORS FROM INPUT FRAME TO OUTPUT FRAME
C
      CALL STDCOP(IMNOA,IMNOC,3,' ',STAT) 
      CALL STSEPI
      END

C
C-----------------------------------------------------------------------------
C
      SUBROUTINE RPLBP(A,B,G,NX,NY)
C
C  SUBROUTINE TO REPLACE BAD PIXELS BY NEAREST GOOD PIXEL
C
      IMPLICIT NONE
      INTEGER NX,NY
      INTEGER NL,NC,I
      INTEGER ZNEX(11,11),SNEX(11,11)
      INTEGER M,Z,S
      DOUBLE PRECISION A(1),B(1),G(1) 
      DOUBLE PRECISION AA(256,256),BB(256,256)
      INTEGER GOODPX(256,256)
      INTEGER ZZ,SS,ZTRY,STRY
C     REAL THRMIN,THRMAX
C
C  READ INDEX OF NEXT NEIGHBOURS
C
      M=11 
C     CALL RDNXTN(M,ZNEX,SNEX)
      CALL NXTZS(M,ZNEX,SNEX)
C     DO Z=1,M
C     WRITE (6,*) (ZNEX(Z,S),SNEX(Z,S),S=1,M)
C     END DO
C
C  REWRITE 1D IMAGE A INTO 2D IMAGE AA
C
      I=0
      DO NL=1,NX
         DO NC=1,NY
            I=I+1
            AA(NL,NC)=A(I)
            GOODPX(NL,NC)=INT(G(I))
         END DO 
      END DO  
C
C  DO ALL THE PROCESSING OF DATA HERE
C
C**********************************************************************
C
C  SUBSTITUTE BADPIXEL BY NEXT GOOD PIXEL
C
      DO Z=1,NX
         DO S=1,NY
            IF (GOODPX(Z,S).EQ.0) THEN
              DO ZZ=1,M
                DO SS=1,M
                  ZTRY=Z+ZNEX(ZZ,SS)
                  STRY=S+SNEX(ZZ,SS)
                  IF (ZTRY.LT.1) ZTRY=1
                  IF (ZTRY.GT.NX) ZTRY=NX
                  IF (STRY.LT.1) STRY=1
                  IF (STRY.GT.NY) STRY=NY
                  IF (GOODPX(ZTRY,STRY).EQ.1) GOTO 1
                END DO
              END DO
1             CONTINUE
              BB(Z,S)=AA(ZTRY,STRY)
            ELSE
              BB(Z,S)=AA(Z,S)
            ENDIF
         END DO 
      END DO   
C**********************************************************************
C
C  REWRITE 2D IMAGE BB INTO 1D IMAGE B
C
      I=0
      DO NL=1,NX
         DO NC=1,NY
            I=I+1
            B(I)=BB(NL,NC)
        END DO 
      END DO  
      RETURN
      END
C
C-----------------------------------------------------------------------------
C
      SUBROUTINE RPLBPS(A,B,G,NX,NY)
C
C  SUBROUTINE TO CONVERT REAL A DOUBLE PRECISSION AD
C
      IMPLICIT NONE
      INTEGER NX,NY
      INTEGER NTOT,I
      DOUBLE PRECISION G(1),B(1) 
      REAL A(1) 
      DOUBLE PRECISION AD(65536)
      NTOT=NX*NY 
      DO I=1,NTOT
            AD(I)=A(I)
      END DO  
      CALL RPLBP(AD,B,G,NX,NY)
      RETURN
      END
C
C-----------------------------------------------------------------------------
C
      SUBROUTINE RDNXTN(M,ZNEX,SNEX)
C
C  READ INDEX OF NEIGHBOURING PIXELS ORDERD ACCORDING TO 
C  INCREASING DISTANCE FROM FILE NEXTN
C  GFI 17.6.92
C 
      IMPLICIT NONE
      INTEGER Z,M
      INTEGER ZNEX(M,M),SNEX(M,M)
      OPEN (UNIT=3,FILE='NEXTN',STATUS='UNKNOWN')
C
C     DO Z=1,M
C     READ (3,*) (ZNEX(Z,S),SNEX(Z,S),S=1,M)
C     WRITE (6,*) (ZNEX(Z,S),SNEX(Z,S),S=1,M)
C     END DO
C
      CALL NXTZS(M,ZNEX,SNEX)
C
      DO Z=1,M
C     READ (3,*) (ZNEX(Z,S),SNEX(Z,S),S=1,M)
C     WRITE (6,*) (ZNEX(Z,S),SNEX(Z,S),S=1,M)
      END DO
      RETURN 
      END
C
C-----------------------------------------------------------------------------
C
      SUBROUTINE NXTZS(M,ZNEX,SNEX)
C
C  CREATES INDEX ZNEX,SNEX OF NEIGHBOURING PIXELS ORDERD ACCORDING TO 
C  INCREASING DISTANCE USING HEAPSORT ALGORITHM
C  GFI 18.6.92
C 
      IMPLICIT NONE
      INTEGER  Z,S,M,C,I,NTOT
      INTEGER  ZNEX(M,M),SNEX(M,M)
      REAL     DIST(11,11),DIST1(121)
      INTEGER  ZNEX1(121),SNEX1(121),ORDER(121)
      NTOT=M*M
      DO Z=1,M 
C     WRITE (6,*) (ZNEX(Z,S),SNEX(Z,S),S=1,M)
      END DO
C
C  DETERMINE DISTANCE FROM CENTER PIXEL
C
      C=INT(M/2.0)+1
      DO Z=1,M
        DO S=1,M
          ZNEX(Z,S)=Z-C
          SNEX(Z,S)=S-C
          DIST(Z,S)=SQRT(1.0*ZNEX(Z,S)*ZNEX(Z,S)+
     2                   1.0*SNEX(Z,S)*SNEX(Z,S))
C         WRITE (6,*) ZNEX(Z,S),SNEX(Z,S),DIST(Z,S)        
        END DO
      END DO     
C
C  REWRITE ZNEX,SNEX,DIST TO 1D ARRAYS
C
      I=0
      DO Z=1,M
        DO S=1,M
          I=I+1
          ZNEX1(I)=ZNEX(Z,S)
          SNEX1(I)=SNEX(Z,S)
          DIST1(I)=DIST(Z,S)
          ORDER(I)=I
        END DO
      END DO
C
C  GENERATE INDEX VECTOR ORDER FOR SORTING ZNEX,SNEX ACCORDING TO
C  ASCENDING DISTANCE 
C 
      CALL INDEXX(NTOT,DIST1,ORDER)      
C
C  REWRITE AND ORDER ZNEX,SNEX,DIST TO 2D ARRAYS ACCORDING TO 
C  INCREASING DISTANCE
C
      I=0
      DO Z=1,M
        DO S=1,M
          I=I+1
          ZNEX(Z,S)=ZNEX1(ORDER(I))
          SNEX(Z,S)=SNEX1(ORDER(I))
          DIST(Z,S)=DIST1(ORDER(I))
C         WRITE (6,*) ZNEX(Z,S),SNEX(Z,S),DIST(Z,S)        
        END DO
      END DO
      RETURN 
      END
C
C-----------------------------------------------------------------------------
C
      SUBROUTINE INDEXX(N,ARRIN,INDX)
C
C  CREATES INDEX VECTOR INDX ACCORDING TO INCREASING ARRIN 
C  USING HEAPSORT ALGORITHM ( NUMERICAL RECIPES H. PRESS p 233)
C  GFI 18.6.92
C
      IMPLICIT NONE
      INTEGER  N
      REAL     ARRIN(N),Q
      INTEGER  INDX(N)
      INTEGER  I,J,L,IR,INDXT
C
      DO J=1,N
        INDX(J)=J
      END DO
      L=N/2+1
      IR=N
10    CONTINUE
        IF (L.GT.1) THEN
          L=L-1
          INDXT=INDX(L)
          Q=ARRIN(INDXT)
        ELSE
          INDXT=INDX(IR)
          Q=ARRIN(INDXT)
          INDX(IR)=INDX(1)
          IR=IR-1
          IF (IR.EQ.1) THEN
            INDX(1)=INDXT
            RETURN
          ENDIF
        ENDIF
        I=L
        J=L+L
20      IF (J.LE.IR) THEN
          IF (J.LT.IR) THEN
            IF (ARRIN(INDX(J)).LT.ARRIN(INDX(J+1))) J=J+1
          ENDIF
          IF (Q.LT.ARRIN(INDX(J))) THEN
            INDX(I)=INDX(J)
            I=J
            J=J+J
          ELSE
            J=IR+1
          END IF
        GOTO 20
        ENDIF
        INDX(I)=INDXT
      GOTO 10     
      END