C @(#)necsord.for 17.1.1.1 (ESO-DMD) 01/25/02 17:51:31 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 PROGRAM ECHSRD C C+ C.COPYRIGHT: Copyright (c) 1987 European Southern Observatory, C all rights reserved C C.VERSION: 1.0 ESO-FORTRAN Conversion, AA 23:20 - 3 DEC 1987 C C.LANGUAGE: F77+ESOext C C.AUTHOR: D.PONZ C C.IDENTIFICATION C C program ECHSORD version 1.6 830310 C C.MODIFICATIONS C C 910128 P.Ballester Define and set variables C C.KEYWORDS C C ECHELLE, CASPEC C C.PURPOSE C C SEARCH FOR ORDERS IN AN IMAGE. TYPE ECHELLE - FLAT FIELD C OUTPUT IS A TABLE WITH THE FOLLOWING COLUMNS : C COLUMN NO. LABEL UNIT DESCRIPTION C 1 :ORDER - ORDER NUMBER C 2 :X PIXEL SAMPLE NUMBER C 3 :Y PIXEL 1. MOMENT AS SUM(F*Y)/SUM(F) C 4 :YBKG PIXEL POSITION OF THE BACKGROUND C 5 :BKG DN BACKGROUND LEVEL C C.ALGORITHM C C MAXIMA FOLLOWING ALGORITHM IS USED C C CONSTRAINS : MAXIMUM NUMBER OF ORDERS IN THE IMAGE IS 100 C FOR EACH ORDER 19 POINTS ARE FOUND C C.INPUT/OUTPUT C C COMMAND : C SEARCH/ORDER INPUT W,T,S OUTPUT METH C input keywords are C ECHC(1:8) INPUT - IMAGE NAME C ECHC(21:28) OUTPUT - TABLE NAME C W,T,S - PARAMETERS WHERE : C ECHR(1) W IS THE ESTIMATED WIDTH OF THE ORDER ACCROSS THE C DISPERSION DIRECTION IN USER COORDINATES C ECHR(2) T THRESHOLD IN THE DETECTION C ECHR(4) S SLOPE OF THE ORDER IN THE MIDDLE OF THE C IMAGE TAG(ALFA) C ECHC(17:20) METH - STD STANDARD STAR METHOD (DEFAULT) C COM COMPLEMENT C C 010621 last modif C C- IMPLICIT NONE C INTEGER MADRID INTEGER ACTVAL,I,IADD,IDELTA,IPOINT INTEGER IST,ISTAT,NACOL,NAUX INTEGER NAXIS,NDIM,NORDER,NTOT INTEGER NWIN1,NWIN2,PNTRA,STATUS,TID C INTEGER*8 IXCOL,IYAUX C INTEGER NPIX(3),ICOL(7),IAUX(300),KUN,KNUL,IMNI INTEGER NSTART(300),NDATA(300),NPOINT(300),IHW(300) C REAL RPAR(4),LEVEL,SLOPE REAL ORDPOS(300),XAUX(300) C DOUBLE PRECISION START(3),STEP(3) C CHARACTER*8 FRMIN,TABLE CHARACTER*72 IDENA CHARACTER*80 LINE CHARACTER*64 CUNITA CHARACTER*2 METH CHARACTER*16 LABCOL(5),UNIT(5) CHARACTER*6 FORM(5) C INCLUDE 'MID_INCLUDE:ST_DEF.INC' COMMON /VMR/MADRID(1) INCLUDE 'MID_INCLUDE:ST_DAT.INC' C DATA LABCOL(1)/'ORDER'/,FORM(1)/'I6'/ DATA LABCOL(2)/'X'/,FORM(2)/'F7.1'/ DATA LABCOL(3)/'Y'/,FORM(3)/'F7.1'/ DATA LABCOL(4)/'YBKG'/,FORM(4)/'F7.1'/ DATA LABCOL(5)/'BKG'/,FORM(5)/'E12.3'/ DATA UNIT(1)/'UNITLESS'/ DATA UNIT(2)/'PIXEL'/ DATA UNIT(3)/'PIXEL'/ DATA UNIT(4)/'PIXEL'/ DATA UNIT(5)/'DN'/ DATA NACOL/5/ C NDIM = 3 C C initialize system C CALL STSPRO('ECHSRD') CALL STKRDC('ECHC',1,1,8,ACTVAL,FRMIN,KUN,KNUL,STATUS) CALL STKRDC('ECHC',1,17,4,ACTVAL,METH,KUN,KNUL,STATUS) C CALL SXKGET('ECHC', 'C', 21, 8, ACTVAL, TABLE, STATUS) TABLE = 'ORDER' CALL STKRDR('ECHR',1,4,ACTVAL,RPAR,KUN,KNUL,STATUS) LEVEL = RPAR(2) SLOPE = RPAR(4) IF (METH.EQ.'CO') THEN LABCOL(3) = 'YBKG' LABCOL(4) = 'Y' END IF CALL STTPUT(' ECHELLE DEFINITION',ISTAT) CALL STTPUT(' ------------------',ISTAT) CALL STTPUT(' INPUT IMAGE : '//FRMIN,ISTAT) CALL STTPUT(' OUTPUT TABLE : '//TABLE,ISTAT) CALL STTPUT(' PARAMETERS ',ISTAT) WRITE (LINE,9000) RPAR(1) CALL STTPUT(LINE,ISTAT) WRITE (LINE,9010) RPAR(2) CALL STTPUT(LINE,ISTAT) WRITE (LINE,9020) SLOPE CALL STTPUT(LINE,ISTAT) IF (METH.EQ.'CO') THEN LINE = ' METHOD : COMPLEMENT' ELSE LINE = ' METHOD : STANDARD' END IF CALL STTPUT(LINE,ISTAT) C C read frame C CALL STIGET(FRMIN,D_R4_FORMAT,F_I_MODE,F_IMA_TYPE, . NDIM,NAXIS,NPIX,START,STEP,IDENA, . CUNITA,PNTRA,IMNI,STATUS) NWIN2 = RPAR(1) NWIN1 = NWIN2/3 + 1 C C extract column C CALL TDMGET(4*NPIX(2),IXCOL,IST) CALL TDMGET(4*NPIX(2),IYAUX,IST) IPOINT = NPIX(1)/2 CALL EXTR3(MADRID(PNTRA),NPIX(1),NPIX(2),IPOINT,1,NPIX(2), + MADRID(IXCOL)) IF (METH.EQ.'CO') CALL VCOM(NPIX(2),MADRID(IXCOL),MADRID(IXCOL)) C C find order positions in the middle of the image C NAUX = 200 CALL FINDM1(MADRID(IXCOL),NPIX(2),NWIN1,NWIN2,LEVEL,NSTART,NDATA, + NORDER,MADRID(IYAUX),XAUX,IAUX,NAUX) CALL TDMFRE(4*NPIX(2),IXCOL,IST) CALL TDMFRE(4*NPIX(2),IYAUX,IST) DO 10 I = 1,NORDER ORDPOS(I) = FLOAT(NSTART(I)) + FLOAT(NDATA(I)-1)*0.5 10 CONTINUE WRITE (LINE,9030) NORDER CALL STTPUT(LINE,ISTAT) C C get order positions C IDELTA = (NPIX(1)-IPOINT)/10 IF (NORDER.EQ.1) THEN IHW(1) = MIN0(NPIX(2),2*NDATA(1)) ELSE DO 20 I = 1,NORDER - 1 IHW(I) = (ORDPOS(I+1)-ORDPOS(I)) 20 CONTINUE IHW(NORDER) = IHW(NORDER-1) END IF C C allocate table C NTOT = NORDER*30 CALL TBTINI(TABLE,F_TRANS,F_O_MODE,12,NTOT,TID,STATUS) DO 30 I = 1,NACOL CALL TBCINI(TID,D_R4_FORMAT,1,FORM(I),UNIT(I), . LABCOL(I),ICOL(I),STATUS) 30 CONTINUE CALL TBCMAP(TID,0,IADD,STATUS) C C follow maxima C CALL FOLLOW(MADRID(PNTRA),NPIX(1),NPIX(2),IPOINT,ORDPOS,NORDER, + IDELTA,IHW,SLOPE,NWIN1,NWIN2,LEVEL,MADRID(IADD),NTOT, + NACOL,NPOINT,METH) C C update number of elements C NTOT = 0 DO 40 I = 1,NORDER NTOT = NTOT + NPOINT(I) 40 CONTINUE CALL STKWRI('ECHI',NORDER,4,1,KUN,STATUS) CALL STKWRI('ECHI',NPIX,7,2,KUN,STATUS) CALL TBIPUT(TID,NACOL,NTOT,STATUS) IDENA = 'ORDER POSITION' CALL STDWRC(TID,'IDENT',1,IDENA,1,72,KUN,STATUS) C CALL DSCUPX(TABLEF,TABLEF,' ',STATUS) CALL TBSINI(TID,STATUS) CALL TBTCLO(TID,STATUS) C C free data C CALL STSEPI 9000 FORMAT (' ORDER WIDTH : ',F6.1) 9010 FORMAT (' THRESHOLD : ',F6.1) 9020 FORMAT (' SLOPE : ',F6.3) 9030 FORMAT (' NUMBER OF DETECTED ORDERS : ',I4) END SUBROUTINE VCOM(N,X,Y) C C COMPUTE COMPLEMENT IN VECTOR FORM C IMPLICIT NONE INTEGER N,I REAL X(N),Y(N) DO 10 I = 1,N Y(I) = -X(I) 10 CONTINUE RETURN END SUBROUTINE FOLLOW(X,N1,N2,IXST,ORDPOS,NORDER,IDELTA,IHW,SLOPE, + NWI1,NWI2,THRES,TAB,NTOT,NACOL,NP,METH) IMPLICIT NONE C C C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ C C.IDENTIFICATION C SUBROUTINE FOLLOW VERSION 1 28 FEB 1983 C J.D.PONZ C C.KEYWORDS: C SPECTROSCOPY, CASPEC C C.PURPOSE: C DEFINE ORDER POSTIONS BY FOLLOWING STRUCTURES ABOVE A GIVEN THRESHOLD C C.ALGORITHM: C C ORDER POSITIONS ARE CALCULATED BY THE ROUTINE FINDLT C C.INPUT/OUTPUT C INPUT ARGUMENTS C C X(N1, N2) REAL*4 INPUT ARRAY C N1, N2 INTG*4 DIMENSIONS OF X C IXST INTG*4 STARTING X POSITION IN PIXELS C ORDPOS(NORDER)REAL*4 Y POSITION OF THE ORDERS IN PIXELS C NORDER INTG*4 NUMBER OF ORDERS C IDELTA INTG*4 X INTERVAL IN PIXELS C IHW(NORDER) INTG*4 HALFWIDTH OF THE Y-INTERVAL TO BE SEARCHED C SLOPE REAL INITIAL SLOPE OF THE ORDERS C NWI INTG*4 NUMBER OF PIXELS IN THE WINDOW C THRES REAL THRESHOLD C NTOT INTG*4 ALLOCATED SPACE C C OUTPUT ARGUMENTS C C TAB REAL ARRAY WITH THE LOCATED POSITIONS IN FORMAT C (NTOT,0:NACOL) 1 - ORDER NUMBER C 2 - SAMPLE C 3 - MOMENT ORDER 1 C 4 - position of back. C 5 - BACKGROUND C NPOS INTG NUMBER OF LOCATED POSITIONS C C------------------------------------------------------------- C C INTEGER N1,N2,NLIM,I,IXST,NORDER,IDELTA,NWI1,NWI2 INTEGER NTOT,NACOL,NP,ISTAT,IORD,NPF,IDELT1,NL INTEGER NPIX,IX0,IY01,IY0,KZ6907,NPIX1,NMAX,I1,I2 INTEGER IB1,IB2,NR,KZ6910 REAL SLOPE,THRES,SLOPE1,SLMAX,SLMIN,BACK,AVE,Y,A REAL SLOPE2,F REAL X(N1,N2),ORDPOS(NORDER),XPIX(300) REAL XL(50),YL(50),XAUX(300),YAUX(300) INTEGER IHW(NORDER),IAUX(300),NSTART(10),NDATA(10) REAL TAB(NTOT,0:NACOL),TRNULL DOUBLE PRECISION TDNULL INTEGER TINULL CHARACTER*2 METH CHARACTER*80 LINE C CALL TBMNUL(TINULL,TRNULL,TDNULL) C NLIM = N1 - 2 NP = 0 C C ... DISPLAY RESULTS C CALL STTPUT(' ',ISTAT) WRITE (LINE,9020) CALL STTPUT(LINE,ISTAT) WRITE (LINE,9030) CALL STTPUT(LINE,ISTAT) C C ... ITERATION ON ORDERS C DO 90 IORD = 1,NORDER C C ... INCLUDE START + END POINTS C NPF = 0 NP = NP + 1 TAB(NP,1) = IORD TAB(NP,2) = 1. TAB(NP,3) = TRNULL TAB(NP,4) = TRNULL TAB(NP,5) = TRNULL NP = NP + 1 TAB(NP,1) = IORD TAB(NP,2) = N1 TAB(NP,3) = TRNULL TAB(NP,4) = TRNULL TAB(NP,5) = TRNULL C C ... LEFT PART OF THE IMAGE C IDELT1 = IDELTA NL = 1 SLOPE1 = SLOPE SLMAX = SLOPE SLMIN = SLOPE NPIX = 2*IHW(IORD) IX0 = IXST IY01 = ORDPOS(IORD) - IHW(IORD) + 0.5 IY0 = MAX0(1,IY01) XL(1) = IX0 + IDELT1 YL(1) = ORDPOS(IORD) + IDELT1*SLOPE1 NPIX = MIN0(NPIX,N2-IY0) DO 20 KZ6907 = 1,2147483647 IF ( .NOT. (IX0.GE.2.AND.NPIX.GT.NWI2)) GO TO 30 C C ... EXTRACT COLUMN NPIX1 = NPIX - (IY0-IY01) CALL EXTR3(X,N1,N2,IX0,IY0,NPIX1,XPIX) IF (METH.EQ.'CO') CALL VCOM(NPIX1,XPIX,XPIX) C C C ... FIND MAXIMUM IF (METH.EQ.'ST' .OR. METH.EQ.'CO') THEN CALL FINDM1(XPIX,NPIX1,NWI1,NWI2,THRES,NSTART, + NDATA,NMAX,YAUX,XAUX,IAUX,300) ELSE CALL FINDM5(XPIX,NPIX1,NWI2,THRES,NSTART,NDATA, + NMAX) END IF IF (NMAX.EQ.1) THEN NL = NL + 1 XL(NL) = IX0 YL(NL) = FLOAT(NSTART(1)) + FLOAT(NDATA(1)-1)*0.5 + + IY0 - 1. IDELT1 = IDELTA C C ... COMPUTE MOMENTA NP = NP + 1 I1 = NSTART(1) I2 = NSTART(1) + NDATA(1) - 1 BACK = (XPIX(I1-1)+XPIX(I1-2)+XPIX(I1-3)+XPIX(I2+1)+ + XPIX(I2+2)+XPIX(I2+3))/6. AVE = 0. Y = 0. DO 10 I = I1,I2 F = XPIX(I) - BACK AVE = AVE + F Y = Y + F*I 10 CONTINUE TAB(NP,1) = IORD TAB(NP,2) = IX0 A = Y/AVE NPF = NPF + 1 IF (IORD.NE.1 .OR. IX0.EQ.IXST) THEN TAB(NP,3) = A + IY0 - 1 ELSE TAB(NP,3) = TRNULL END IF IB1 = NINT(A-IHW(IORD)/2.) IB2 = NINT(A+IHW(IORD)/2.) IF (IB2.LE.NPIX .AND. IB2.GE.1) THEN IF (METH.EQ.'CO') THEN TAB(NP,5) = -XPIX(NINT(A)) ELSE TAB(NP,5) = XPIX(IB2) END IF TAB(NP,4) = IB2 + IY0 - 1 ELSE TAB(NP,5) = TRNULL TAB(NP,4) = TRNULL END IF C ENDIF ELSE NP = NP + 1 TAB(NP,1) = IORD TAB(NP,2) = IX0 TAB(NP,3) = TRNULL TAB(NP,4) = TRNULL TAB(NP,5) = TRNULL IDELT1 = IDELT1 + IDELTA END IF IX0 = XL(NL) - IDELT1 IF (NL.GT.2) THEN SLOPE2 = (YL(NL)-YL(NL-1))/ (XL(NL)-XL(NL-1)) C IF (ABS(SLOPE2-SLOPE1).GT.ABS(SLOPE1)*1.0) THEN C TYPE *,SLOPE1,SLOPE2 C NP = NP - 1 C GO TO 10 C ENDIF SLOPE1 = SLOPE2 SLMAX = MAX(SLMAX,SLOPE1) SLMIN = MIN(SLMIN,SLOPE1) END IF IY01 = YL(NL) - SLOPE1*IDELT1 - IHW(IORD) + 0.5 IY0 = MAX0(1,IY01) NPIX = 2*IHW(IORD) NPIX = MIN0(NPIX,N2-IY0) 20 CONTINUE 30 CONTINUE C C ... RIGTH PART OF THE IMAGE C IDELT1 = IDELTA NR = 1 SLOPE1 = SLOPE XL(1) = IXST YL(1) = ORDPOS(IORD) NPIX = 2*IHW(IORD) IX0 = IXST + IDELT1 IY01 = ORDPOS(IORD) + IDELT1*SLOPE1 - IHW(IORD) + 0.5 IY0 = MAX0(1,IY01) NPIX = MIN0(NPIX,N2-IY0) DO 60 KZ6910 = 1,2147483647 IF ( .NOT. (IX0.LT.NLIM.AND.NPIX.GT.NWI2)) GO TO 70 C C ... EXTRACT COLUMN NPIX1 = NPIX - (IY0-IY01) CALL EXTR3(X,N1,N2,IX0,IY0,NPIX1,XPIX) IF (METH.EQ.'CO') CALL VCOM(NPIX1,XPIX,XPIX) C C ... FIND MAXIMUM IF (METH.EQ.'ST' .OR. METH.EQ.'CO') THEN CALL FINDM1(XPIX,NPIX1,NWI1,NWI2,THRES,NSTART, + NDATA,NMAX,YAUX,XAUX,IAUX,300) ELSE CALL FINDM5(XPIX,NPIX1,NWI2,THRES,NSTART,NDATA, + NMAX) END IF IF (NMAX.EQ.1) THEN NR = NR + 1 XL(NR) = IX0 YL(NR) = FLOAT(NSTART(1)) + FLOAT(NDATA(1)-1)*0.5 + + IY0 - 1. IDELT1 = IDELTA C C ... COMPUTE MOMENTA NP = NP + 1 I1 = NSTART(1) I2 = NSTART(1) + NDATA(1) - 1 BACK = (XPIX(I1-1)+XPIX(I1-2)+XPIX(I1-3)+XPIX(I2+1)+ + XPIX(I2+2)+XPIX(I2+3))/6. AVE = 0. Y = 0. DO 50 I = I1,I2 F = XPIX(I) - BACK AVE = AVE + F Y = Y + F*I 50 CONTINUE TAB(NP,1) = IORD TAB(NP,2) = IX0 A = Y/AVE NPF = NPF + 1 IF (IORD.NE.NORDER) THEN TAB(NP,3) = A + IY0 - 1 ELSE TAB(NP,3) = TRNULL END IF IB1 = NINT(A-IHW(IORD)/2.) IB2 = NINT(A+IHW(IORD)/2.) IF (IB2.LE.NPIX .AND. IB2.GE.1) THEN IF (METH.EQ.'CO') THEN TAB(NP,5) = -XPIX(NINT(A)) ELSE TAB(NP,5) = XPIX(IB2) END IF TAB(NP,4) = IB2 + IY0 - 1 ELSE TAB(NP,5) = TRNULL TAB(NP,4) = TRNULL END IF ELSE NP = NP + 1 TAB(NP,1) = IORD TAB(NP,2) = IX0 TAB(NP,3) = TRNULL TAB(NP,4) = TRNULL TAB(NP,5) = TRNULL IDELT1 = IDELT1 + IDELTA END IF IX0 = XL(NR) + IDELT1 IF (NR.GT.2) THEN SLOPE2 = (YL(NR)-YL(NR-1))/ (XL(NR)-XL(NR-1)) C IF (ABS(SLOPE2-SLOPE1).GT.ABS(SLOPE1)*1.0) THEN C TYPE *,SLOPE1,SLOPE2 C NP = NP - 1 C GO TO 20 C ENDIF SLOPE1 = SLOPE2 SLMAX = MAX(SLMAX,SLOPE1) SLMIN = MIN(SLMIN,SLOPE1) END IF IY01 = YL(NR) + SLOPE1*IDELT1 - IHW(IORD) + 0.5 IY0 = MAX0(1,IY01) NPIX = 2*IHW(IORD) NPIX = MIN0(NPIX,N2-IY0) 60 CONTINUE 70 CONTINUE WRITE (LINE,9010) IORD,IXST,NINT(ORDPOS(IORD)),IHW(IORD), + NPF,SLMAX,SLMIN CALL STTPUT(LINE,ISTAT) 90 CONTINUE WRITE (LINE,9040) CALL STTPUT(LINE,ISTAT) RETURN C 9000 FORMAT (I4) 9010 FORMAT (1X,I7,2X,I8,2X,I8,2X,I10,2X,I8,2X,F8.2,2X,F8.2) 9020 FORMAT (' SEQ.NO. X CENTER Y CENTER INTERORDER TEMPLA',' MAX', + '.SLOPE MIN.SLOPE') 9030 FORMAT (' ------- -------- -------- ---------- --------',' -', + '-------- ---------') 9040 FORMAT (' -------------------------------------------------','--', + '------------------') END SUBROUTINE EXTR3(X,N1,N2,IXS,IYS,N,XCOL) C+ C.COPYRIGHT: Copyright (c) 1987 European Southern Observatory, C all rights reserved C C.VERSION: 1.0 ESO-FORTRAN Conversion, AA 23:45 - 3 DEC 1987 C C.LANGUAGE: F77+ESOext C C.AUTHOR: D.PONZ C C.KEYWORDS C C C.PURPOSE C C EXTRACTS A COLUMN IN A 2D ARRAY AS AVERAGE OF 3 CONTIGUOUS COLUMNS C C.ALGORITHM C C OUTPUT COLUMN IS THE AVERAGE OF THREEE CONTIGUOUS COLUMNS C C.INPUT/OUTPUT C C INPUT ARGUMENTS ARE C X(N1,N2) REAL*4 INPUT 2D FRAME C N1,N2 INTG*4 X-DIMENSION C IXS INTG*4 CENTRAL COLUMN NUMBER C IYS INTG*4 STARTING ROW NUMBER C N INTG*4 NUMBER OF PIXELS IN THE COLUMN C C OUTPUT ARGUMENTS ARE C XCOL(N) REAL*4 EXTRACTED COLUMN C C- C C EXTRACT COLUMN AT IPOINT C C IMPLICIT NONE INTEGER N1,N2,N,I1,IXS,IYS,I,J,I2,J1 REAL X(N1,N2),XCOL(N) I1 = IXS - 1 I2 = IXS + 1 I = IXS DO 10 J = 1,N J1 = J + IYS - 1 XCOL(J) = (X(I1,J1)+X(I,J1)+X(I2,J1))/3.0 10 CONTINUE RETURN END C SUBROUTINE FINDM1(X,N,NW1,NW2,THRES,IST,NDAT,NPOS,YAUX,XAUX,IAUX, + NAUX) C+ C.COPYRIGHT: Copyright (c) 1987 European Southern Observatory, C all rights reserved C C.VERSION: 1.0 ESO-FORTRAN Conversion, AA 23:50 - 3 DEC 1987 C C.LANGUAGE: F77+ESOext C C.AUTHOR: D.PONZ C C.KEYWORDS C C CASPEC C C.PURPOSE C C FIND VALUES ABOVE A GIVEN THRESHOLD C C.ALGORITHM C C C.INPUT/OUTPUT C C INPUT ARGUMENTS C X(N) REAL INPUT ARRAY C N INTG ARRAY DIMENSION C NW1 INTG SMALL WINDOW C NW2 INTG LARGE WINDOW C THRES REAL THRESHOLD C C OUTPUT ARGUMENTS C IST INTG STARTING POSITION OF THE REGION C NDAT INTG NUMBER OF PIXELS IN THE REGION C NPOS NPOS NUMBER OF REGIONS ABOVE THRESHOLD C C AUXILIARY SPACE C YAUX(N) REAL C XAUX(NAUX) REAL C IAUX(NAUX) INTG C C- C C IMPLICIT NONE REAL X(N),YAUX(N),XAUX(NAUX),THRES INTEGER IST(1),NDAT(1),IAUX(NAUX),N,I,NFIRST INTEGER NW1,NW2,NPIX,NPOS,NAUX,MWIND,NO LOGICAL FIRST C C ... median filtering using threshold C MWIND = 2*NW2 + 1 CALL MEDNTH(X,N,MWIND,THRES,YAUX,XAUX,IAUX,NAUX) C C ... original - filtered C FIRST = .TRUE. NO = 0 DO 10 I = NW2,N - NW2 + 1 YAUX(I) = X(I) - YAUX(I) IF (YAUX(I).GT.0.) THEN C ... above the threshold IF (FIRST) THEN C ... first value above the threshold FIRST = .FALSE. NFIRST = I NPIX = 1 ELSE NPIX = NPIX + 1 END IF ELSE C ... below the threshold IF ( .NOT. FIRST) THEN C ... first value below the threshold FIRST = .TRUE. IF (NPIX.GE.NW1) THEN NO = NO + 1 IST(NO) = NFIRST NDAT(NO) = NPIX END IF END IF END IF 10 CONTINUE NPOS = NO RETURN END SUBROUTINE FINDM5(X,N,NW,T,IST,NDAT,NACT) C+ C.COPYRIGHT: Copyright (c) 1987 European Southern Observatory, C all rights reserved C C.VERSION: 1.0 ESO-FORTRAN Conversion, AA 23:56 - 3 DEC 1987 C C.LANGUAGE: F77+ESOext C C.AUTHOR: D.PONZ C C SUBROUTINE FINDM5 VERSION 1.0 11 MAR 1983 C J.D.PONZ ESO - GARCHING C C.KEYWORDS C C FIND SPECTRAL FEATURES, SPECTROSCOPY C C.PURPOSE C C FIND FLAT FIELD PROFILES C C.ALGORITHM C C SIMPLE THRESHOLDING METHOD C C.INPUT/OUTPUT C C INPUT ARGUMENTS C X(N) REAL INPUT ARRAY C N INTG ARRAY DIMENSION C NW INTG WINDOW SIZE IN PIXELS C T REAL THRESHOLD C C OUTPUT ARGUMENTS C IST INTG POINTER OF THE STARTING POSITION C NDAT INTG NUMBER OF PIXELS ABOVE THE THRESHOLD C NACT INTG NUMBER OF DETECTED FEATURES C C- C IMPLICIT NONE REAL X(N),T INTEGER IST(1),NDAT(1),N,NW,NACT INTEGER I,NPIX,NO LOGICAL BACK C BACK = .FALSE. NO = 0 DO 10 I = 1,N IF (X(I).GT.T) THEN IF (BACK) THEN BACK = .FALSE. NO = NO + 1 IST(NO) = I NPIX = 1 ELSE NPIX = NPIX + 1 END IF ELSE IF ( .NOT. BACK) THEN BACK = .TRUE. IF (NO.GT.0) NDAT(NO) = NPIX END IF END IF 10 CONTINUE NACT = NO IF (NO.EQ.2) THEN NACT = 1 ELSE IF (NO.EQ.3) THEN NACT = 1 NDAT(1) = NDAT(2) IST(1) = IST(2) END IF END IF RETURN END SUBROUTINE MEDNTH(IA,LEN,WIDTH,THRES,OA,FA,LINK,NAUX) C+ C C.IDENTIFICATION C SUBROUTINE MEDIAN VERSION 1.0 8 MAR 1983 C C.COPYRIGHT: Copyright (c) 1987 European Southern Observatory, C all rights reserved C C.VERSION: 1.0 ESO-FORTRAN Conversion, AA 0:34 - 4 DEC 1987 C C.LANGUAGE: F77+ESOext C C.AUTHOR: D.PONZ C C.KEYWORDS C C FILTERING, MEDIAN FILTERING, SPECTROSCOPY, THRESHOLDING C C.PURPOSE C C PERFORM MEDIAN FILTERING AND THRESHOLDING ON THE INPUT ARRAY C C.ALGORITHM C C DATA IN A RUNNING WINDOW IS SORTED TO FORM A LOCAL HISTOGRAM C THE VALUE IN THE MIDDLE OF THE LOCAL HISTOGRAM IS TAKEN AS C OUTPUT VALUE IF ABS(MEDIAN-VALUE) .GT. ABS(THRES) C C.INPUT/OUTPUT C C INPUT ARGUMENTS : C IA REAL*4 INPUT ARRAY C LEN INTG*4 NUMBER OF INPUT PIXELS C WIDTH INTG*4 WINDOW WIDTH C THRES REAL*4 THRESHOLD C C OUTPUT ARGUMENTS : C OA REAL*4 OUTPUT ARRAY C C WORK SPACE : C FA(NAUX) REAL*4 AUXILIARY SPACE C LINK(NAUX) INTG*4 AUXILIARY SPACE C NAUX INTG*4 SPACE SIZE. NAUX .GT. WIDTH+1 C C- C IMPLICIT NONE REAL FA(NAUX),IA(LEN),OA(LEN),THRES,T INTEGER LINK(NAUX),MIN,MAX,WIDTH,LEN,MIDPT,FRONT,AVAIL INTEGER RANGE,NEWNOD,I,J,K,NODE,ITMP,N,NAUX C C SET UP BEGINNING (MIN) AND END (MAX) OF FILTER ARRAY C ALSO SET UP POINTER TO AVAIL LIST AND LIST HEADER (FA(1),LINK(1) C T = ABS(THRES) MIN = 2 MAX = NAUX AVAIL = MIN FA(1) = -1.0 LINK(1) = 0 C C SET UP FIRST FILTER WINDOW C RANGE = WIDTH - 1 DO 10 I = 1,RANGE,1 NODE = NEWNOD(IA(I),AVAIL,MIN,MAX,FA,LINK,NAUX) CALL INSNOD(NODE,FA,LINK,NAUX) 10 CONTINUE C C COPY BORDERS TO OUTPUT BUFFER (CANNOT PROCESS BORDER POINTS) C RANGE = WIDTH/2 K = LEN - RANGE DO 20 J = 1,RANGE,1 OA(J) = IA(J) K = K + 1 C K POINTS TO HIGH END OF ARRAY, J POINTS TO LOW OA(K) = IA(K) 20 CONTINUE C C FOR EACH POINT FROM WIDTH/2 TO LEN-WIDTH/2+1 C FIND ITS MEDIAN C FRONT = MIN C FRONT IS THE OLDEST POINT IN THE FILTER ARRAY K = RANGE C K IS THE POINT BEING PROCESSED (LOC IN OA) MIDPT = RANGE + 1 C MIDPT IS THE MIDDLE OF THE WINDOW (IN FA) DO 40 J = WIDTH,LEN,1 NODE = NEWNOD(IA(J),AVAIL,MIN,MAX,FA,LINK,NAUX) CALL INSNOD(NODE,FA,LINK,NAUX) K = K + 1 C FIND THE MIDPT - - IT IS THE MEDIAN VALUE ITMP = LINK(1) DO 30 I = 2,MIDPT,1 ITMP = LINK(ITMP) 30 CONTINUE C GET THE LINK WHICH POINTS TO THE VALUE OF THE MEDIAN C MODIFY ONLY IF GREATHER THAN THRESHOLD IF (ABS(FA(ITMP)-IA(K)).GT.T) THEN OA(K) = FA(ITMP) ELSE OA(K) = IA(K) END IF CALL DELNOD(FRONT,FA,LINK,NAUX) FRONT = FRONT + 1 IF (FRONT.GT.MAX) FRONT = MIN C DELETE OLDEST AND ADVANCE PTR, FRONT ALSO WRAPS 40 CONTINUE C C USE SMOOTH VALUE AT THE ENDS C N = WIDTH/2 + 1 DO 50 I = 1,N OA(I) = OA(N) OA(LEN-I+1) = OA(LEN-N) 50 CONTINUE RETURN END