C @(#)irsdefbad.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 .for C.AUTHOR E. Oliva, Firenze-Arcetri C.KEYWORDS Spectroscopy, IRSPEC C C.PURPOSE Execute command ... C C.ALGORITHM C C C.INPUT/OUTPUT C C C.VERSION 1.0 Creation 02.09.1992 E. Oliva C C------------------------------------------------------- C PROGRAM DEFBAD IMPLICIT REAL(A-H,O-Z) IMPLICIT INTEGER(I-N) CHARACTER*60 FRAMEI,TABLEO CHARACTER*1 SMODE CHARACTER*64 CUNIT CHARACTER*72 IDENT REAL RTH(2) INTEGER*8 INPNTR,OUPNTR INTEGER NPIX(2) DOUBLE PRECISION START(2),STEP(2) COMMON /VMR/ MADRID(1) INCLUDE 'MID_INCLUDE:ST_DEF.INC' INCLUDE 'MID_INCLUDE:ST_DAT.INC' DATA MAXDIM/2/ IRET=1 CALL STSPRO('DEFBAD') C C GET NAME OF INPUT FRAME AND MAP IT C CHECK ALSO WHETHER IT IS A 2D FRAME..... C CALL STKRDC('framei',1,1,60,IRET,FRAMEI,KUNIT,KNUL,ISTAT) 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 Create dummy image wich will contain badpixels. C CALL STIPUT('middumma',D_R4_FORMAT,F_O_MODE,F_IMA_TYPE, , NAXIS,NPIX,START,STEP, , IDENT,CUNIT,OUPNTR,NOU,ISTAT) C C GET FLAG TO DEFINE the selection mode C CALL STKRDC('smode',1,1,1,IRET,SMODE,KUNIT,KNUL,ISTAT) C C Get threshold value(s) C CALL STKRDR('thres',1,2,IRET,RTH,KUNIT,KNUL,ISTAT) C C GET KEY TO DEFINE the number of subframes C CALL STKRDI('nsub',1,1,IRET,NSUB,KUNIT,KNUL,ISTAT) C C GET THE NAME OF THE OUTPUT TABLE WHICH WILL CONTAIN THE COORDINATES OF THE C BAD PIXELS C CALL STKRDC('tableo',1,1,60,IRET,TABLEO,KUNIT,KNUL,ISTAT) C C OPEN TABLE FOR UPDATE (IOPEN=2) C IOPEN=2 CALL TBTOPN(TABLEO,IOPEN,IDETAB,ISTAT) C C GET NUMBERS CORRESPONDING TO COLUMNS C CALL TBLSER(IDETAB,'X_COORD',NXCOL,ISTAT) CALL TBLSER(IDETAB,'Y_COORD',NYCOL,ISTAT) CALL TBLSER(IDETAB,'AVALUE',NAVAL,ISTAT) CALL TBLSER(IDETAB,'RVALUE',NRVAL,ISTAT) C C CALL ROUTINE WHICH DEFINE BAD_PIXELS AND WRITE THEM IN THE TABLE C CALL FIND(MADRID(INPNTR),MADRID(OUPNTR), , NPIX(1),NPIX(2),NSUB, , IDETAB,NXCOL,NYCOL,NAVAL,NRVAL,SMODE,RTH) C C Close table C CALL TBTCLO(IDETAB,ISTAT) C C RELEASE FILES, UPDATE KEYWORDS AND EXIT C CALL STSEPI END C C C SUBROUTINE FIND (A,B,NX,NY,NSUB, , IDETAB,NXCOL,NYCOL,NAVAL,NRVAL,SMODE,RTH) C IMPLICIT REAL(A-H,O-Z) IMPLICIT INTEGER(I-N) C CHARACTER*1 SMODE,DMODE(2),CMODE(2) DIMENSION A(NX,NY),B(NX,NY),RTH(1) CHARACTER*200 OUTPUT C DATA DMODE/'r','a'/ DATA CMODE/'R','A'/ C C Check INPUT IMAGE (A) C SUM=0.0 SUM2=0.0 DO J=1,NY DO I=1,NX SUM=SUM+A(I,J) SUM2=SUM2+A(I,J)*A(I,J) ENDDO ENDDO STD=SUM2-SUM*SUM/(NX*NY) IF(STD.LE.0) THEN WRITE(OUTPUT,192) 192 FORMAT('Input image is a constant, check it please') CALL STTPUT(OUTPUT,ISTAT) RETURN ENDIF C DO I=1,2 IMODE=I IF(SMODE.EQ.DMODE(I) .OR. SMODE.EQ.CMODE(I)) GO TO 10 ENDDO CALL STETER +(1,'Invalid selection mode, only -r-, -a- are recognized') 10 CONTINUE C C Set badpixel_image (B) to zero C DO J=1,NY DO I=1,NX B(I,J)=0. ENDDO ENDDO NBAD=0 C C Divide image in 16 quadrants and define standard deviation in C each of them. C IX1=1 IX2=NX/NSUB DO IXQ=1,NSUB IY1=1 IY2=NY/NSUB DO IYQ=1,NSUB CALL STD_ITER(A,NX,NY,IX1,IX2,IY1,IY2,AVE,STD) DO IY=IY1,IY2 DO IX=IX1,IX2 VAL=A(IX,IY) RVAL=(VAL-AVE)/STD IF(IMODE.EQ.1) THEN IF(ABS(RVAL).GT.RTH(1)) THEN NBAD=NBAD+1 R4=IX CALL TBRWRR(IDETAB,NBAD,1,NXCOL,R4,ISTAT) R4=IY CALL TBRWRR(IDETAB,NBAD,1,NYCOL,R4,ISTAT) CALL TBRWRR(IDETAB,NBAD,1,NAVAL,VAL,ISTAT) CALL TBRWRR(IDETAB,NBAD,1,NRVAL,RVAL,ISTAT) B(IX,IY)=A(IX,IY) ENDIF ENDIF IF(IMODE.EQ.2) THEN IF(VAL.LT.RTH(1).OR.VAL.GT.RTH(2)) THEN NBAD=NBAD+1 R4=IX CALL TBRWRR(IDETAB,NBAD,1,NXCOL,R4,ISTAT) R4=IY CALL TBRWRR(IDETAB,NBAD,1,NYCOL,R4,ISTAT) CALL TBRWRR(IDETAB,NBAD,1,NAVAL,VAL,ISTAT) CALL TBRWRR(IDETAB,NBAD,1,NRVAL,RVAL,ISTAT) B(IX,IY)=A(IX,IY) ENDIF ENDIF ENDDO ENDDO IY1=IY2+1 IY2=IY1+NY/NSUB IY2=MIN(IY2,NY) ENDDO IX1=IX2+1 IX2=IX1+NX/NSUB IX2=MIN(IX2,NX) ENDDO WRITE(OUTPUT,195) NBAD 195 FORMAT(I5,' bad pixels found and stored in table.') CALL STTPUT(OUTPUT,ISTAT) RETURN END SUBROUTINE STD_ITER(A,NX,NY,IX1,IX2,IY1,IY2,AVE,STD) C C Define the std of a given region of the image (IX1,IX2:IY1,IY2) C by iteratively eliminating the points which are more than 3 C sigma from the average C IMPLICIT REAL(A-H,O-Z) IMPLICIT INTEGER(I-N) DIMENSION A(NX,NY) DIMENSION IXOUT(5000),IYOUT(5000) NOUT=0 IXOUT(1)=0 IYOUT(1)=0 10 CONTINUE SUM=0.0 SUM2=0.0 NAVE=0 DO IY=IY1,IY2 DO IX=IX1,IX2 DO J=1,NOUT IF(IX.EQ.IXOUT(J) .AND. IY.EQ.IYOUT(J)) GO TO 11 ENDDO VAL=A(IX,IY) SUM=SUM+VAL SUM2=SUM2+VAL*VAL NAVE=NAVE+1 11 CONTINUE ENDDO ENDDO IF(NAVE.LE.0) THEN CALL STTPUT('---------------------------------------',ISTAT) CALL STTPUT(' NAVE=0, problems in STD_ITER',ISTAT) CALL STTPUT(' You probably chose too many sub_frames.',ISTAT) CALL STTPUT('---------------------------------------',ISTAT) CALL STETER(1,'Aborted') ENDIF AVE=SUM/NAVE STD=SQRT((SUM2-NAVE*AVE*AVE)/(NAVE-1)) IF(STD.LE.0) CALL STETER(1,'STD<=0, problems in STD_ITER') N=0 DO IY=IY1,IY2 DO IX=IX1,IX2 DIST=ABS((A(IX,IY)-AVE)/STD) IF(DIST.GT.3) THEN N=N+1 IXOUT(N)=IX IYOUT(N)=IY ENDIF ENDDO ENDDO IF(N.GT.NOUT) THEN NOUT=N GO TO 10 ENDIF RETURN END