C @(#)crpixmask.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
C 
      PROGRAM CREATEPIXMASK
C++++++++++++++++++++++++++++++++++++++++++++++++++++
C.COPYRIGHT: Copyright (c) 1992 European Southern Observatory,
C                                         all rights reserved
C.IDENTIFICATION: CREATEPIXMASK
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:   920617 GXF Created
C.VERSION:   920713 RHW Header added
C---------------------------------------------------------------------
      IMPLICIT NONE
      INTEGER      NAXISA,NPIXA(2),IAV,STAT
      INTEGER      IMNOA,IMNOC,IMNOG
      INTEGER*8    PNTRA,PNTRC,PNTRG
      INTEGER      KNULL,KUNIT(1)
      INTEGER      MADRID(1)
      INTEGER      FNO,IBUF(5)
      CHARACTER*60 FRAMEA,FRAMEC,FRAMEG
      CHARACTER    CUNITA*64,IDENTA*72,DEFAUL*1
      DOUBLE PRECISION STEPA(2),STARTA(2)
      REAL LOWMASK,HIGHMASK,INPUTR(2)
C 
      INCLUDE 'MID_INCLUDE:ST_DEF.INC'   
C 
      COMMON /VMR/ MADRID
C 
      INCLUDE 'MID_INCLUDE:ST_DAT.INC'   
C 
C 
      CALL STSPRO('CREATEPIXMASK')
      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
C  DOUBLE PRECISION
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)
      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  CREATE AND MAP IMAGE GOODPX ( IMAGE OF GOOD PIXELS )
C
      FRAMEG='GOODPX'
      CALL STIPUT(FRAMEG,D_R8_FORMAT,F_O_MODE,F_IMA_TYPE,
     +            NAXISA,NPIXA,STARTA,STEPA,IDENTA,
     +            CUNITA,PNTRG,IMNOG,STAT) 
C
C  GET LOW AND HIGH VALUES TO DEFINE BAD PIXELS
C
      CALL STKRDR('INPUTR',1,2,IAV,INPUTR,KUNIT,KNULL,STAT)
      LOWMASK=INPUTR(1)
      HIGHMASK=INPUTR(2)
C
C  GET FLAG IF BAD PIXELS ARE REPLACED           : R
C              OR BAD PIXELS ARE SET TO HIGHMASK : H
C
      CALL STKRDC('DEFAULT',1,1,1,IAV,DEFAUL,KUNIT,KNULL,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 CRBP(MADRID(PNTRA),MADRID(PNTRC),MADRID(PNTRG),
     +                NPIXA(1),NPIXA(2),LOWMASK,HIGHMASK,DEFAUL)
      ENDIF 
C
C  SINGLE PRECISSION
C
      IF (IBUF(1).EQ.4) THEN 
      CALL CRBPS(MADRID(PNTRA),MADRID(PNTRC),MADRID(PNTRG),
     +                NPIXA(1),NPIXA(2),LOWMASK,HIGHMASK,DEFAUL)
      ENDIF 
C
C  COPY ALL DESCRIPTORS FROM INPUT FRAME TO OUTPUT FRAME
C
      CALL STDCOP(IMNOA,IMNOC,3,' ',STAT)
      CALL STDCOP(IMNOA,IMNOG,3,' ',STAT) 
      CALL STSEPI
      END

C
C-----------------------------------------------------------------------------
C
      SUBROUTINE CRBP(A,B,G,NX,NY,THRMIN,THRMAX,FLAG)
C
C  SUBROUTINE TO CREATE MASK GOODPX (1:GOOD PIXEL 0:BAD PIXEL)
C  AND TO REPLACE BAD PIXELS BY NEAREST GOOD PIXEL
C  GFI 17.6.92
C   
      IMPLICIT NONE
C 
      INTEGER NX,NY
      INTEGER NL,NC,I
      INTEGER ZNEX(11,11),SNEX(11,11)
      INTEGER M,Z,S
      INTEGER GOODPX(256,256),BAD,GOOD
      INTEGER ZZ,SS,ZTRY,STRY
C 
      DOUBLE PRECISION A(1),B(1),G(1) 
      DOUBLE PRECISION AA(256,256),BB(256,256),GG(256,256)
C 
      CHARACTER FLAG*1
C 
      REAL THRMIN,THRMAX
C
C  READ INDEX OF NEXT NEIGHBOURS
C
      M=11
C     THRMIN=3000
C     THRMAX=8000 
C     CALL RDNXTN(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)
         END DO 
      END DO  
C
C  DO ALL THE PROCESSING OF DATA HERE
C
C**********************************************************************
C
C  DETERMINE PIXEL MASK
C
      BAD=0
      GOOD=0
      DO Z=1,NX
         DO S=1,NY
            IF ((AA(Z,S).LE.THRMIN).OR.(AA(Z,S).GE.THRMAX)) THEN
              GOODPX(Z,S)=0
              GG(Z,S)=0.
              BAD=BAD+1 
            ELSE
              GOODPX(Z,S)=1
              GG(Z,S)=1. 
              GOOD=GOOD+1
            ENDIF
          END DO 
      END DO
C
C  SPECIAL FOR ROCKWELL FIRST SCIENCE ARRAY WHICH HAS A BAD 
C  COLUMN IN 2ND QUADRAND     
C
      DO Z=129,256
        DO S=230,231
          GOODPX(Z,S)=0
          GG(Z,S)=0.
        END DO
      END DO
      DO Z=128,134
        DO S=228,235
          GOODPX(Z,S)=0
          GG(Z,S)=0.
        END DO
      END DO
      DO S=230,232
        GOODPX(127,S)=0 
        GG(127,S)=0
      END DO
      GOODPX(134,235)=1
      GG(134,235)=1
      GOODPX(128,235)=1
      GG(128,235)=1    
C
C
      WRITE (6,*) GOOD,' GOOD PIXELS',BAD,' BAD PIXELS'
      WRITE (6,*) BAD*1.0/(GOOD+BAD)*100,' % BAD PIXELS'
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
              IF (FLAG.EQ.'R') THEN 
                BB(Z,S)=AA(ZTRY,STRY)
              ELSE
                BB(Z,S)=THRMAX*4.
              ENDIF
            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)
            G(I)=GG(NL,NC)
         END DO 
      END DO  
      RETURN
      END
C
C-----------------------------------------------------------------------------
C
      SUBROUTINE CRBPS(A,B,G,NX,NY,THRMIN,THRMAX,FLAG)
C
C  SUBROUTINE TO CONVERT ARRAY A FROM SINGLE PRECISSION TO DOUBLE
C  PRECISSION    
   
      IMPLICIT NONE
      INTEGER  NX,NY
      INTEGER  NTOT,I
      DOUBLE PRECISION B(1),G(1) 
      REAL    A(1) 
      CHARACTER FLAG*1
      REAL     THRMIN,THRMAX
      DOUBLE PRECISION AD(65536)
C
      NTOT=NX*NY 
      DO I=1,NTOT
            AD(I)=A(I)
      END DO  
      CALL CRBP(AD,B,G,NX,NY,THRMIN,THRMAX,FLAG)
      RETURN
      END 
C
C-----------------------------------------------------------------------------
C
      SUBROUTINE RDNXTN(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  M
      INTEGER  ZNEX(11,11),SNEX(11,11)
      M=11
C     OPEN (UNIT=3,FILE='NEXTN')
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
      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
C     DO Z=1,M
C     WRITE (6,*) (ZNEX(Z,S),SNEX(Z,S),S=1,M)
C     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
      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