C @(#)rcosmic.for 13.1.1.1 (ES0-DMD) 06/02/98 18:17:32 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.COPYRIGHT: Copyright (c) 1987 European Southern Observatory, C all rights reserved C.IDENTIFICATION: RCOSMIC C.PURPOSE: Remove cosmic ray events on single ccd exposure and replace C them by interpolation on neighbourhood pixels. C.LANGUAGE: F77+ESOext C.AUTHOR: P.Magain, M.Remy, Institut d'Astrophysique de Liege C.INPUT/OUTPUT: IN_A/C/1/60 = Input frame. C IN_B/C/1/60 = median filter of the input frame C OUT_B/C/1/12 = Output frame. C PARAMS/R/1/1 = Mean value of the background. C PARAMS/R/2/1 = Readout noise in ADU units. C PARAMS/R/3/1 = Inverse gain factor (e-/ADU). C VERSION 1.2 New ST interfaces M.Peron 221191 C VERSION 1.3 remove limitation of size, suggestion from Hans Schwengeler C ---------------------------------------------------------------------------- PROGRAM RCOSMIC INTEGER*4 IAV,I INTEGER*4 STATUS,MADRID,SIZEX,IOMODE INTEGER*4 NAXIS,NPIX(2),IMNI,IMNO,IMNF,IMNC INTEGER*8 PNTRI,PNTRF,PNTRO,PNTRC INTEGER*4 KUN,KNUL CHARACTER*60 IMAGE,OBJET,COSMIC CHARACTER*72 IDENT1,IDENT2,IDENT3 CHARACTER*48 CUNIT DOUBLE PRECISION START(2),STEP(2) REAL*4 SKY,GAIN,RON,NS,RC,PARAM(5),CUTS(2) INCLUDE 'MID_INCLUDE:ST_DEF.INC' COMMON/VMR/MADRID(1) INCLUDE 'MID_INCLUDE:ST_DAT.INC' DATA IDENT1 /' '/ DATA IDENT2 /' '/ DATA IDENT3 /'cosmic ray mask '/ DATA CUNIT /' '/ CALL STSPRO('RCOSMIC') CALL STKRDC('IN_A',1,1,60,IAV,IMAGE,KUN,KNUL,STATUS) CALL STIGET(IMAGE,D_R4_FORMAT,F_I_MODE,F_IMA_TYPE, 1 2,NAXIS,NPIX,START,STEP 1 ,IDENT1,CUNIT,PNTRI,IMNI,STATUS) CALL STKRDR('PARAMS',1,5,IAV,PARAM,KUN,KNUL,STATUS) CALL STIGET('middumma',D_R4_FORMAT,F_I_MODE,F_IMA_TYPE, 1 2,NAXIS,NPIX,START,STEP 1 ,IDENT2,CUNIT,PNTRF,IMNF,STATUS) SKY = PARAM(1) GAIN = PARAM(2) RON = PARAM(3) NS = PARAM(4) RC = PARAM(5) CALL STKRDC('OUTIMA',1,1,60,IAV,OBJET,KUN,KNUL,STATUS) CALL STIPUT(OBJET,D_R4_FORMAT,F_O_MODE,F_IMA_TYPE, 1 NAXIS,NPIX,START,STEP 1 ,IDENT1,CUNIT,PNTRO,IMNO,STATUS) SIZEX = 1 DO I=1,NAXIS SIZEX = SIZEX*NPIX(I) ENDDO CALL STKRDC('COSMIC',1,1,60,IAV,COSMIC,KUN,KNUL,STATUS) IF (COSMIC(1:1).EQ.'+') THEN COSMIC = 'dummy_frame' IOMODE = F_X_MODE ELSE IOMODE = F_O_MODE ENDIF CALL STIPUT(COSMIC,D_I2_FORMAT,IOMODE,F_IMA_TYPE 1 ,NAXIS,NPIX,START,STEP 1 ,IDENT3,CUNIT,PNTRC,IMNC,STATUS) CALL COSROUT(MADRID(PNTRI),MADRID(PNTRC),NPIX(1),NPIX(2), 1 RON,GAIN,NS,SKY,RC 1 ,MADRID(PNTRF),MADRID(PNTRO)) CUTS(1) = 0 CUTS(2) = 1 IF (IOMODE.EQ.F_O_MODE) + CALL STDWRR(IMNC,'LHCUTS',CUTS,1,2,KUN,STATUS) CALL DSCUPT(IMNI,IMNO,' ',STATUS) CALL STSEPI END SUBROUTINE COSROUT(AI,COSMIC,I_IMA,J_IMA,RON,GAIN, 1 NS,SKY,RC,AM,AO) INTEGER I_IMA,J_IMA,NUM INTEGER ORD(10000) INTEGER K,L INTEGER IDUMAX,JDUMAX,I1,I2,J1,II,JJ INTEGER I,J,IMAX,JMAX,IMIN,JMIN INTEGER FIRST(2),NEXT(2) INTEGER*2 COSMIC(I_IMA,J_IMA) REAL*4 VECTEUR(10000),FMAX,ASUM,RC REAL*4 AI(I_IMA,J_IMA),AO(I_IMA,J_IMA),AM(I_IMA,J_IMA) REAL*4 SIGMA,SKY,S1,S2 REAL*4 RON,GAIN,NS,AMEDIAN DO 10 J=1,J_IMA DO 5 I=1,I_IMA AO(I,J)=AI(I,J) COSMIC(I,J)= 0 5 CONTINUE 10 CONTINUE C C La boucle suivante selectionne les pixels qui sont C significativement au dessus de l'image filtree medianement. C C C COSMIC =-1 ----> candidate for cosmic C = 0 ----> not a cosmic C = 1 -----> a cosmic (at the end) C = 2 ----> member of the group C = 3 ----> member of a group which has been examined C = 4 ----> neighbourhood of the group K=1 DO 80 J=2,J_IMA-1 DO 70 I=2,I_IMA-1 SIGMA=SQRT(RON**2+AM(I,J)/GAIN) IF ((AI(I,J)-AM(I,J)).GE.(NS*SIGMA)) THEN COSMIC(I,J) = -1 K = K+1 ENDIF 70 CONTINUE 80 CONTINUE C C ------------------------ C C C C C Ces pixels sont regroupes par ensembles connexes dans la boucle C C C ----------------------------- FIRST(1) = 2 FIRST(2) = 2 100 CALL FINDNEXT(COSMIC,I_IMA,J_IMA,FIRST,NEXT) IF (NEXT(1).EQ.-1) RETURN I = NEXT(1) J = NEXT(2) COSMIC(I,J) = 2 IMIN = I IMAX = I JMIN = J JMAX = J IDUMAX = I JDUMAX = J FMAX = AI(I,J) 110 I1 = 0 I2 = 0 CONTINUE DO 125 L=1,2 DO 115 K=1,2 II = I+K-L JJ = J+K+L-3 IF (COSMIC(II,JJ).EQ.-1) THEN I1 = II J1 = JJ IMIN = MIN(IMIN,II) IMAX = MAX(IMAX,II) JMIN = MIN(JMIN,JJ) JMAX = MAX(JMAX,JJ) COSMIC(II,JJ) = 2 IF (AI(II,JJ).GT.FMAX) THEN FMAX = AI(II,JJ) IDUMAX = II JDUMAX = JJ ENDIF ELSE IF (COSMIC(II,JJ).EQ.0) THEN COSMIC(II,JJ) = 4 ENDIF 115 CONTINUE 125 CONTINUE COSMIC(I,J) = 3 IF (I1.NE.0) THEN I = I1 J = J1 GOTO 110 ENDIF DO 140 L = JMIN,JMAX DO 130 K = IMIN,IMAX IF (COSMIC(K,L).EQ.2) THEN I = K J = L GOTO 110 ENDIF 130 CONTINUE 140 CONTINUE FIRST(1) = NEXT(1)+1 FIRST(2) = NEXT(2) C We start here the real work.... C 1- decide if the pixel's group is a cosmic C 2-replace these values by another one S1 = AI(IDUMAX-1,JDUMAX-1) + AI(IDUMAX+1,JDUMAX-1) + + AI(IDUMAX,JDUMAX-1) +AI(IDUMAX,JDUMAX+1) S2 = AI(IDUMAX-1,JDUMAX+1) + AI(IDUMAX+1,JDUMAX+1) + + AI(IDUMAX-1,JDUMAX) + AI(IDUMAX+1,JDUMAX) ASUM = (S1+S2)/8.-SKY IF ((FMAX-SKY).GT.RC*ASUM) THEN NUM = 1 DO L = JMIN-1,JMAX+1 DO K = IMIN-1,IMAX+1 IF (COSMIC(K,L).EQ.4) THEN VECTEUR(NUM) = AI(K,L) NUM = NUM+1 ENDIF ENDDO ENDDO CALL SORT(NUM-1,VECTEUR,ORD) AMEDIAN = VECTEUR(ORD((NUM-1)/2)) DO L = JMIN-1,JMAX+1 DO K = IMIN-1,IMAX+1 IF (COSMIC(K,L).EQ.3) THEN COSMIC(K,L) = 1 AO(K,L) = AMEDIAN ELSE IF (COSMIC(K,L).EQ.4) THEN COSMIC(K,L) = 0 ENDIF ENDDO ENDDO ELSE DO L = JMIN-1,JMAX+1 DO K = IMIN-1,IMAX+1 IF (COSMIC(K,L).NE.-1) COSMIC(K,L) = 0 ENDDO ENDDO ENDIF IF (NEXT(1).GT.0) GOTO 100 C C C RETURN END SUBROUTINE FINDNEXT(COSMIC,I_IMA,J_IMA,FIRST,NEXT) INTEGER I_IMA,J_IMA,FIRST(2),NEXT(2) INTEGER I,J INTEGER*2 COSMIC(I_IMA,J_IMA) DO J = FIRST(2), J_IMA DO I = 2, I_IMA IF (COSMIC(I,J).EQ.-1) THEN NEXT(1) = I NEXT(2) = J RETURN ENDIF ENDDO ENDDO NEXT(1) = -1 NEXT(2) = -1 RETURN END SUBROUTINE SORT(KMAX,INP,ORD) INTEGER KMAX,IMIN,IMAX,I,J,K,L INTEGER ORD(10000) REAL*4 INP(10000),F DO 4100 J=1,KMAX ORD(J)=J 4100 CONTINUE IF (INP(1).GT.INP(2)) THEN ORD(1)=2 ORD(2)=1 END IF DO 4400 J=3,KMAX F=INP(J) L=ORD(J-1) IF (INP(L).LE.F) GO TO 4400 L=ORD(1) IMIN=1 IF (F.LE.INP(L)) GO TO 4250 IMAX=J-1 4200 I=(IMIN+IMAX)/2 L=ORD(I) IF (INP(L).LT.F) THEN IMIN=I ELSE IMAX=I END IF IF (IMAX.GT.(IMIN+1)) GO TO 4200 IMIN=IMAX 4250 DO 4300 K=J-1,IMIN,-1 ORD(K+1)=ORD(K) 4300 CONTINUE ORD(IMIN)=J 4400 CONTINUE RETURN END