C @(#)necavmed.for 17.1.1.1 (ES0-DMD) 01/25/02 17:51:28 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 @(#)necavmed.for 17.1.1.1 (ESO) 01/25/02 17:51:28 C C+ C.COPYRIGHT: Copyright (c) 1991 European Southern Observatory, C all rights reserved C C.VERSION: 1.0 23-JULY-1991 C C.LANGUAGE: F77+ESOext C C.AUTHOR: P.BALLESTER C C.IDENTIFICATION C C.KEYWORDS C C ECHELLE, CASPEC, BLAZE FUNCTION C C.PURPOSE C C COMPARE SUCCESSIVE ORDERS TO GET THE MEDIAN AND AVERAGE BLAZE FUNCTION C C.ALGORITHM C C IN DEVELOPMENT ... C C.INPUT/OUTPUT C C- PROGRAM AVEMED IMPLICIT NONE INTEGER NAXISA,NPIXA(2),IAV,STAT,ACTVAL,MAXORD INTEGER*8 PNTRA,PNTRB INTEGER IMNOA,IMNOB INTEGER KNULL,KUNIT(1) INTEGER MADRID(1) PARAMETER (MAXORD=500) CHARACTER FRAMEA*60,FRAMEB*60 CHARACTER CUNIT*64,IDENTA*72 INTEGER ORDSTA(MAXORD),ORDEND(MAXORD),DELTAY REAL THRES DOUBLE PRECISION STEPA(2),STARTA(2) INCLUDE 'MID_INCLUDE:st_def.inc' COMMON /VMR/ MADRID INCLUDE 'MID_INCLUDE:st_dat.inc' CALL STSPRO('AVEMED') CALL STKRDC('IN_A',1,1,60,IAV,FRAMEA,KUNIT,KNULL,STAT) CALL STIGET(FRAMEA,D_R4_FORMAT,F_I_MODE,F_IMA_TYPE, + 2,NAXISA,NPIXA,STARTA,STEPA,IDENTA,CUNIT, + PNTRA,IMNOA,STAT) IF (NPIXA(2).LT.MAXORD) THEN CALL STDRDI(IMNOA,'ORDSTA',1,NPIXA(2),ACTVAL,ORDSTA, + KUNIT,KNULL,STAT) CALL STDRDI(IMNOA,'ORDEND',1,NPIXA(2),ACTVAL,ORDEND, + KUNIT,KNULL,STAT) ELSE CALL STETER(10,'Buffer overflow in AVEMED.') ENDIF CALL STKRDC('OUT_A',1,1,60,IAV,FRAMEB,KUNIT,KNULL,STAT) CALL STIPUT(FRAMEB,D_R4_FORMAT,F_IO_MODE,F_IMA_TYPE, + NAXISA,NPIXA,STARTA,STEPA,IDENTA,CUNIT, + PNTRB,IMNOB,STAT) CALL STKRDI('INPUTI',1,1,IAV,DELTAY,KUNIT,KNULL,STAT) IF (DELTAY.GT.0) THEN CALL STKRDR('INPUTR',1,1,IAV,THRES,KUNIT,KNULL,STAT) CALL MEDIAN(MADRID(PNTRA),NPIXA(1),NPIXA(2),MADRID(PNTRB), + ORDSTA,ORDEND,DELTAY,THRES) ELSE CALL AVERAGE(MADRID(PNTRA),NPIXA(1),NPIXA(2),MADRID(PNTRB), + ORDSTA,ORDEND) ENDIF CALL STSEPI END C ======================Median Filtering==================== SUBROUTINE MEDIAN(INPFRAM,NX,NY,OUTFRAM,ORDSTA,ORDEND, + DELY,THRES) IMPLICIT NONE INTEGER NX,NY,DELY,ROW,COL INTEGER ORDSTA(NY),ORDEND(NY) REAL INPFRAM(NX,NY),OUTFRAM(NX,NY),THRES,RATIO INTEGER STAMIN,STAMAX,ENDMIN,ENDMAX INTEGER ROWINF,ROWSUP,FLAG,ROWCUR,Y REAL TRANS,BUFFER(-50:50) IF (DELY.GT.50) THEN CALL STETER(9,'Width too big in AVEMED. Limited to 50.') ENDIF C --- Initialisations STAMIN = ORDSTA(1) STAMAX = ORDSTA(1) ENDMIN = ORDEND(1) ENDMAX = ORDEND(1) ROWINF = 0 ROWSUP = 0 DO 5 ROW = 1,NY IF (ORDSTA(ROW).LT.STAMIN) STAMIN = ORDSTA(ROW) IF (ORDSTA(ROW).GT.STAMAX) STAMAX = ORDSTA(ROW) IF (ORDEND(ROW).LT.ENDMIN) ENDMIN = ORDEND(ROW) IF (ORDEND(ROW).GT.ENDMAX) ENDMAX = ORDEND(ROW) 5 CONTINUE DO 10 COL = 1,NX C --- Find out what is ROWINF and ROWSUP for column COL C --- This part is easier to program with a DO WHILE. Consider C language. IF (COL.GE.STAMAX.AND.COL.LE.ENDMIN) THEN ROWINF = 1 ROWSUP = NY ELSE IF (COL.LT.STAMAX) THEN ROWINF = 0 ROWSUP = 0 DO 20 ROW = 1,NY IF (ORDSTA(ROW).LT.COL.AND.ROWINF.EQ.0) ROWINF = ROW 20 CONTINUE DO 25 ROW = NY,1,-1 IF (ORDSTA(ROW).LT.COL.AND.ROWSUP.EQ.0) ROWSUP = ROW 25 CONTINUE ELSE IF (COL.GT.ENDMIN) THEN ROWINF = 0 ROWSUP = 0 DO 30 ROW = 1,NY IF (ORDEND(ROW).GT.COL.AND.ROWINF.EQ.0) ROWINF = ROW 30 CONTINUE DO 35 ROW = NY,1,-1 IF (ORDEND(ROW).GT.COL.AND.ROWSUP.EQ.0) ROWSUP = ROW 35 CONTINUE ENDIF C --- Fill in the buffer with defined values, prolonging the frame by a C --- mirror effect. DO 40 ROW = ROWINF,ROWSUP DO 50 Y = (-1)*DELY , DELY ROWCUR = ROW+Y IF (ROWCUR.LT.ROWINF) + BUFFER(Y) = INPFRAM(COL,2*ROWINF-ROWCUR) IF (ROWCUR.GT.ROWSUP) + BUFFER(Y) = INPFRAM(COL,2*ROWSUP-ROWCUR) IF (ROWCUR.GE.ROWINF.AND.ROWCUR.LE.ROWSUP) + BUFFER(Y) = INPFRAM(COL,ROWCUR) 50 CONTINUE C --- Sort the buffer values and get the median value. C --- Consider using HeapSort algorithm (Numerical Recipes) 60 CONTINUE FLAG = 0 DO 70 Y = (-1)*DELY , DELY-1 IF (BUFFER(Y).GT.BUFFER(Y+1)) THEN TRANS = BUFFER(Y) BUFFER(Y) = BUFFER(Y+1) BUFFER(Y+1) = TRANS FLAG = 1 ENDIF 70 CONTINUE IF (FLAG.NE.0) GOTO 60 RATIO = (INPFRAM(COL,ROW)-BUFFER(0))/(INPFRAM(COL,ROW) + +BUFFER(0)) IF (ABS(RATIO).GT.THRES) THEN OUTFRAM(COL,ROW) = BUFFER(0) ELSE OUTFRAM(COL,ROW) = INPFRAM(COL,ROW) ENDIF 40 CONTINUE 10 CONTINUE RETURN END C ====================== Average ==================== SUBROUTINE AVERAGE(INPFRAM,NX,NY,OUTFRAM,ORDSTA,ORDEND) IMPLICIT NONE INTEGER NX,NY,ROW,COL INTEGER ORDSTA(NY),ORDEND(NY) REAL INPFRAM(NX,NY),OUTFRAM(NX,NY) INTEGER STAMIN,STAMAX,ENDMIN,ENDMAX INTEGER ROWINF,ROWSUP,CNT REAL SUM DO 5 ROW = 1,NY IF (ORDSTA(ROW).LT.STAMIN) STAMIN = ORDSTA(ROW) IF (ORDSTA(ROW).GT.STAMAX) STAMAX = ORDSTA(ROW) IF (ORDEND(ROW).LT.ENDMIN) ENDMIN = ORDEND(ROW) IF (ORDEND(ROW).GT.ENDMAX) ENDMAX = ORDEND(ROW) 5 CONTINUE DO 10 COL = 1,NX C --- Find out what is ROWINF and ROWSUP for column COL C --- This part is easier to program with a DO WHILE. Consider C language. IF (COL.GE.STAMAX.AND.COL.LE.ENDMIN) THEN ROWINF = 1 ROWSUP = NY ELSE IF (COL.LT.STAMAX) THEN ROWINF = 0 ROWSUP = 0 DO 20 ROW = 1,NY IF (ORDSTA(ROW).LT.COL.AND.ROWINF.EQ.0) ROWINF = ROW 20 CONTINUE DO 25 ROW = NY,1,-1 IF (ORDSTA(ROW).LT.COL.AND.ROWSUP.EQ.0) ROWSUP = ROW 25 CONTINUE ELSE IF (COL.GT.ENDMIN) THEN ROWINF = 0 ROWSUP = 0 DO 30 ROW = 1,NY IF (ORDEND(ROW).GT.COL.AND.ROWINF.EQ.0) ROWINF = ROW 30 CONTINUE DO 35 ROW = NY,1,-1 IF (ORDEND(ROW).GT.COL.AND.ROWSUP.EQ.0) ROWSUP = ROW 35 CONTINUE ENDIF C --- Fill in the buffer with defined values, prolonging the frame by a C --- mirror effect. SUM = 0. DO 40 ROW = ROWINF,ROWSUP SUM = SUM + INPFRAM(COL,ROW) 40 CONTINUE CNT = ROWSUP - ROWINF + 1 SUM = SUM/CNT DO 50 ROW = 1,NY OUTFRAM(COL,ROW) = SUM 50 CONTINUE 10 CONTINUE RETURN END