C @(#)irsbadpix.for	17.1.1.1 (ES0-DMD) 01/25/02 17:53:05
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
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C.COPYRIGHT   (C) 1992 European Southern Observatory
C.IDENT       irsbadpix.for
C.AUTHOR      E. Oliva,  Firenze-Arcetri
C.KEYWORDS    Spectroscopy, IRSPEC
C
C.PURPOSE     Modify values of badpixels in a image (FRAMEI).
C             The coordinates of badpixels is contained inside an input
C             table (REFTAB) in the columns X_COORD,Y_COORD.
C             These coordinates are intended as -- pixel coordinates ---
C             => this program will not work if one define the 'badpixel'
C             table in a 'not-original' frame (i.e. a frame with starts
C             and steps not equl to unity)!
C             The output frame is FRAMEO 
C             This fortran code is called by "irsbadix.prg"
C
C.ALGORITHM   Define the value of each badpixel as the average
C             of either:
C             x: the two neighbouring pixels in X
C             y: the two neighbouring pixels in Y
C             b: the four neighbouring pixels in X+Y
C             In each mode the program uses only the good pixels =>
C             => if the bad-pixels are clustered some of them will
C             not be substituted, in which case a warning message
C             (with pixel coordinates) is printed.
C             The choice between x,y,b is done automatically if the input
C             keyword direct is set to "a" (the default), in this case
C             "b" is used whenever the gradients along X and Y do not
C             different by more than a factor of 2, otherwise take
C             the direction with the smaller gradient.
C             The 'direction option' can also be forced manually to
C             "x","y" or "b.
C
C
C.INPUT/OUTPUT keyword IN_A  --> name of input frame
C              keyword OUT_A --> name of output frame
C              keyword IN_B  --> name of output frame
C              local keyword direct/c/1/1 --> 'direction option', see above
C              local keyword debug/i/1/1 --> if >0 prints many details
C                                            (to be used in case of problems)
C
C
C.VERSION     1.0    Creation     02.09.1992   E. Oliva
C
C-------------------------------------------------------
C
      PROGRAM BADPIX
      IMPLICIT REAL(A-H,O-Z)
      IMPLICIT INTEGER(I-N)

      CHARACTER*60 FRAMEI,FRAMEO,REFTAB
      CHARACTER*1 DIRECT
      CHARACTER*64 CUNIT
      CHARACTER*72 IDENT
      LOGICAL    NULL
      INTEGER*8  INPNTR,OUPNTR
      INTEGER    MMPIX(2)
      INTEGER    NPIX(2),IXBAD(500),IYBAD(500)
      DOUBLE PRECISION START(2),STEP(2)
      REAL CUTS(2)
      COMMON /VMR/ MADRID(1)

      INCLUDE 'MID_INCLUDE:ST_DEF.INC'
      INCLUDE 'MID_INCLUDE:ST_DAT.INC'

      DATA MAXDIM/2/

      IRET=1
      CALL STSPRO('BADPIX')
C
C GET NAME OF INPUT FRAME AND MAP IT
C CHECK ALSO WHETHER IT IS A 2D FRAME.....
C
      CALL STKRDC('IN_A',1,1,60,IRET,FRAMEI,KUNIT,KNUL,ISTAT)
      CALL CLNFRA(FRAMEI,FRAMEI,0)
      CALL STIGET(FRAMEI,D_R4_FORMAT,F_I_MODE,F_IMA_TYPE,MAXDIM,
     ,            NAXIS,NPIX,START,STEP,
     ,            IDENT,CUNIT,INPNTR,NIN,ISTAT)
      IF(NAXIS.NE.2) 
     +     CALL STETER(1,'Input frame must be 2-D')
C
C GET NAME OF OUTPUT FRAME AND MAP IT
C
      CALL STKRDC('OUT_A',1,1,60,IRET,FRAMEO,KUNIT,KNUL,ISTAT)
      CALL CLNFRA(FRAMEO,FRAMEO,0)
      CALL STIPUT(FRAMEO,D_R4_FORMAT,F_O_MODE,F_IMA_TYPE,
     ,            NAXIS,NPIX,START,STEP,
     ,            IDENT,CUNIT,OUPNTR,NOU,ISTAT)
C
C GET FLAG TO DEFINE WHICH NEIGHBOURING PIXELS HAVE TO BE CONSIDERED
C
      CALL STKRDC('direct',1,1,1,IRET,DIRECT,KUNIT,KNUL,ISTAT)
C
C get debug option
C
      CALL STKRDI('debug',1,1,IRET,IDEBUG,KUNIT,KNUL,ISTAT)
C
C OPEN  FOR READING (IOPEN=0)
C THE TABLE CONTAINING THE COORDINATES OF THE
C BAD PIXELS WHICH MUST BE SUBSTITUTED BY THE AVERAGE OF THE
C 2 NEIGHBOURING PIXELS ALONG X (DIRECT=X) Y (DIRECT=Y), 
C THE 4 X+Y PIXELS (DIRECT=B) OR THE AUTOMATIC MODE
C
C
      CALL STKRDC('IN_B',1,1,60,IRET,REFTAB,KUNIT,KNUL,ISTAT)
      CALL CLNTAB(REFTAB,REFTAB,0)
      IOPEN=0
      CALL TBTOPN(REFTAB,IOPEN,IDETAB,ISTAT)
C
C GET NUMBERS CORRESPONDING TO XCOL AND YCOL
C
      CALL TBLSER(IDETAB,'X_COORD',NXCOL,ISTAT)
      CALL TBLSER(IDETAB,'Y_COORD',NYCOL,ISTAT)
C
C GET INFORMATION UPON THE TABLE SIZE (N.B. NROW = # OF BAD PIXELS....)
C
      CALL TBIGET(IDETAB,NCOL,NROW,NSORT,NCALL,NRALL,ISTAT)
      NBAD=NROW
C
C COPY TABLE ELEMENTS INTO VECTORS IXBAD,IYBAD
C
      DO IROW=1,NROW
        CALL TBERDR(IDETAB,IROW,NXCOL,RVAL,NULL,ISTAT)
        IXBAD(IROW)=RVAL
        CALL TBERDR(IDETAB,IROW,NYCOL,RVAL,NULL,ISTAT)
        IYBAD(IROW)=RVAL
      ENDDO
      CALL TBTCLO(IDETAB,ISTAT)
C
C CALL ROUTINE WHICH SUBSTITUTES THE BAD PIXELS LISTED IN THE
C
      CALL SUBSTI(MADRID(INPNTR),MADRID(OUPNTR),
     ,                NPIX(1),NPIX(2),IXBAD,IYBAD,NBAD,DIRECT,IDEBUG)
C
C Define and write cuts of output file
C
      CALL MNMX(MADRID(OUPNTR),NPIX,CUTS,MMPIX)
      CALL STDWRR(NOU,'LHCUTS',CUTS,3,2,KUNIT,ISTAT)
C
C RELEASE FILES, UPDATE KEYWORDS AND EXIT
C
      CALL STSEPI
      END
C
C
C
      SUBROUTINE SUBSTI(A,B,NX,NY,IXBAD,IYBAD,NBAD,DIRECT,IDEBUG)
C
C SUBSTITUTE THE VALUE OF THE BAD PIXELS HAVING (PIXEL) COORDINATES
C (IXBAD(J),IYBAD(J)),J=1,NBAD WITH THE AVERAGE VALUE OF THE 2 NEIGHBOURING
C PIXELS ALONG X (DIRECT=X), ALONG Y (DIRECT=Y) OR OF THE 4 PIXELS ALONG X+Y
C (DIRECT=B).
C In case direct=a makes automatic choice.
C WARNING! THE ROUTINE ASSUMES THAT THE BAD PIXELS ARE MOSTLY ISOLATED.
C IT CAN ALSO DEAL WITH 2-3 CONTIGUOS BAD PIXELS BUT IT MUST NOT BE 
C USED WHEN THE BAD PIXELS ARE CLUSTERED. IN PARTICULAR IT MUST NOT BE USED 
C WHEN MORE THAN 2 BAD PIXELS ARE ALIGNED.
C
      IMPLICIT REAL(A-H,O-Z)
      IMPLICIT INTEGER(I-N)
      CHARACTER*1 DIRECT,DMODE(4),CMODE(4)
      DIMENSION A(NX,NY),B(NX,NY)
      DIMENSION IXBAD(NBAD),IYBAD(NBAD)
      CHARACTER*80 OUTPUT
C
      DATA DMODE/'x','y','b','a'/
      DATA CMODE/'X','Y','B','A'/
C
C COPY INPUT IMAGE (A) INTO OUTPUT IMAGE (B)
C
      TOTVAL=0.0
      NGOOD=0
      NOZERO=0
      DO J=1,NY
        DO I=1,NX
          B(I,J)=A(I,J)
          IFLAG=GOODP(I,J,IXBAD,IYBAD,NBAD,NX,NY)
          IF(IFLAG.GT.0) THEN
            TOTVAL=TOTVAL+A(I,J)
            NGOOD=NGOOD+1
            IF(A(I,J).NE.0) NOZERO=1
          ENDIF
        ENDDO
      ENDDO
      IF(NGOOD.LE.0) THEN
        WRITE(OUTPUT,191)
191     FORMAT(
     + 'No good pixel in image, check your reference table')
         CALL STTPUT(OUTPUT,ISTAT)
        RETURN
      ENDIF
      IF(NOZERO.EQ.0) THEN
        WRITE(OUTPUT,192)
192     FORMAT(
     + 'All good pixel in image have intensities =0, check it please')
        CALL STTPUT(OUTPUT,ISTAT)
        RETURN
      ENDIF
      AVER_PIX=TOTVAL/NGOOD
C
      DO I=1,4
        IMODE=I
        IF(DIRECT.EQ.DMODE(I) .OR. DIRECT.EQ.CMODE(I)) GO TO 10
      ENDDO
      CALL STETER
     +(1,'Invalid direction option, only a x y b are recognized')
10    CONTINUE
      IF(AVER_PIX.EQ.0 .AND. IMODE.EQ.4) THEN
        IMODE=3
        WRITE(OUTPUT,193)
193     FORMAT(
     + 'Average=0, option=b used')
        CALL STTPUT(OUTPUT,ISTAT)
      ENDIF
      DO IBAD=1,NBAD
        I=IXBAD(IBAD)
        J=IYBAD(IBAD)
        CALL SIGXY(A,NX,NY,I,J,IXBAD,IYBAD,NBAD,SIGX,SIGY)
        IF(I.LT.1 .OR. J.LT.1 .OR. I.GT.NX .OR. J.GT.NY) GO TO 20
C
C STORE VALUES OF CONTIGUOS PIXELS AND TAKE CARE OF EDGE EFFECTS
C AND OF CONTIGUOS BAD PIXELS
C
        TOTX=0.
        GX1=GOODP(I-1,J,IXBAD,IYBAD,NBAD,NX,NY)
        TOTX=TOTX+GX1
        X1=0.
        IF(GX1.GT.0.) X1=A(I-1,J)
        GX2=GOODP(I+1,J,IXBAD,IYBAD,NBAD,NX,NY)
        TOTX=TOTX+GX2
        X2=0.
        IF(GX2.GT.0.) X2=A(I+1,J)
        TOTY=0.
        GY1=GOODP(I,J-1,IXBAD,IYBAD,NBAD,NX,NY)
        TOTY=TOTY+GY1
        Y1=0.
        IF(GY1.GT.0.) Y1=A(I,J-1)
        GY2=GOODP(I,J+1,IXBAD,IYBAD,NBAD,NX,NY)
        TOTY=TOTY+GY2
        Y2=0.
        IF(GY2.GT.0.) Y2=A(I,J+1)
        IF(IMODE.LE.3) THEN
          JMODE=IMODE
        ELSE
          IF(TOTX.LE.1 .AND. TOTY.LE.1) THEN
            JMODE=3
            GO TO 11
          ENDIF
          IF(TOTX.LE.1) THEN
            JMODE=2
            GO TO 11
          ENDIF
          IF(TOTY.LE.1) THEN
            JMODE=1
            GO TO 11
          ENDIF
          JMODE=3
          IF(SIGX.GT.(4*SIGY)) JMODE=2
          IF(SIGY.GT.(4*SIGX)) JMODE=1
        ENDIF
11      CONTINUE
        IF(IDEBUG.GT.0) THEN
          WRITE(OUTPUT,119) I,J,DMODE(JMODE)
119       FORMAT('Badpixel (' , I3 , ',' , I3 , ') dir= ', A)
          CALL STTPUT(OUTPUT,ISTAT)
        ENDIF
        IF(JMODE.EQ.1) THEN
          IF(TOTX.EQ.0.) THEN
            WRITE(OUTPUT,101) 
101         FORMAT(
     + '>2 bad pixels along X, (',I2,',',I2,') not modified')
            CALL STTPUT(OUTPUT,ISTAT)
          ELSE
            B(I,J)=(X1+X2)/TOTX
          ENDIF
        ENDIF
        IF(JMODE.EQ.2) THEN
          IF(TOTY.EQ.0.) THEN
            WRITE(OUTPUT,102) I,J
102         FORMAT(
     + '>2 bad pixels along Y, (',I2,',',I2,') not modified')
            CALL STTPUT(OUTPUT,ISTAT)
          ELSE
            B(I,J)=(Y1+Y2)/TOTY
          ENDIF
        ENDIF
        IF(JMODE.EQ.3) THEN
          TOT=TOTX+TOTY
          IF(TOT.EQ.0.) THEN
            WRITE(OUTPUT,103) I,J
103         FORMAT(
     + 'region crowded with bad pixels, (',I2,',',I2,') not modified')
            CALL STTPUT(OUTPUT,ISTAT)
          ELSE
            B(I,J)=(X1+X2+Y1+Y2)/TOT
          ENDIF
        ENDIF
20    CONTINUE
      ENDDO
C
      RETURN
      END
C
C
C
      FUNCTION GOODP(I,J,IBAD,JBAD,NBAD,NX,NY)
C
C CHECK WHETHER PIXEL (I,J) IS A BAD PIXEL...
C   GOODP = 0  BAD PIXEL
C   GOODP = 1  GOOD PIXEL
C N.B. GOODP=0 ALSO WHEN I OR J ARE OUTSIDE THE IMAGE
C
      IMPLICIT REAL(A-H,O-Z)
      IMPLICIT INTEGER(I-N)
      INTEGER IBAD(NBAD),JBAD(NBAD)

      GOODP=1
      IF(I.LT.1 .OR. J.LT.1 .OR. I.GT.NX .OR. J.GT.NY) THEN
        GOODP=0
        RETURN
      ENDIF
      DO K=1,NBAD
        IF(I.EQ.IBAD(K).AND.J.EQ.JBAD(K)) THEN
          GOODP=0
          RETURN
        ENDIF
      ENDDO
      RETURN
      END


      SUBROUTINE SIGXY(A,NX,NY,IX,IY,IXBAD,IYBAD,NBAD,SIGX,SIGY)

      IMPLICIT REAL(A-H,O-Z)
      IMPLICIT INTEGER(I-N)
      REAL A(NX,NY)
      INTEGER IXBAD(NBAD),IYBAD(NBAD)

C
C sigma of the neighbouring 10 pixels along Y --> SIGY
C
      SUM=0.0
      SUM2=0.0
      NTOT=0
      I=IX
      NPT=0
      DO J=IY+1,NY
        IFLAG=GOODP(I,J,IXBAD,IYBAD,NBAD,NX,NY)
        IF(IFLAG.GT.0) THEN
          VAL=A(I,J)
          SUM=SUM+VAL
          SUM2=SUM2+VAL*VAL
          NTOT=NTOT+1
          NPT=NPT+1
          IF(NPT.GE.5) GO TO 10
        ENDIF
      ENDDO
10    CONTINUE
      NPT=0
      DO J=IY-1,1,-1
        IFLAG=GOODP(I,J,IXBAD,IYBAD,NBAD,NX,NY)
        IF(IFLAG.GT.0) THEN
          VAL=A(I,J)
          SUM=SUM+VAL
          SUM2=SUM2+VAL*VAL
          NTOT=NTOT+1
          NPT=NPT+1
          IF(NPT.GE.5) GO TO 11
        ENDIF
      ENDDO
11    CONTINUE
      IF(NTOT.GT.0) THEN
        SIGY=SQRT((SUM2-SUM*SUM/NTOT)/(NTOT-1))
      ELSE
        SIGY=1e20
      ENDIF

C
C sigma of the neighbouring 10 pixels along X --> SIGX
C
      SUM=0.0
      SUM2=0.0
      NTOT=0
      J=IY
      NPT=0
      DO I=IX+1,NX
        IFLAG=GOODP(I,J,IXBAD,IYBAD,NBAD,NX,NY)
        IF(IFLAG.GT.0) THEN
          VAL=A(I,J)
          SUM=SUM+VAL
          SUM2=SUM2+VAL*VAL
          NTOT=NTOT+1
          NPT=NPT+1
          IF(NPT.GE.5) GO TO 20
        ENDIF
      ENDDO
20    CONTINUE
      NPT=0
      DO I=IX-1,1,-1
        IFLAG=GOODP(I,J,IXBAD,IYBAD,NBAD,NX,NY)
        IF(IFLAG.GT.0) THEN
          VAL=A(I,J)
          SUM=SUM+VAL
          SUM2=SUM2+VAL*VAL
          NTOT=NTOT+1
          NPT=NPT+1
          IF(NPT.GE.5) GO TO 21
        ENDIF
      ENDDO
21    CONTINUE
      IF(NTOT.GT.0) THEN
        SIGX=SQRT((SUM2-SUM*SUM/NTOT)/(NTOT-1))
      ELSE
        SIGX=1e20
      ENDIF

      RETURN
      END