C @(#)lnsearch.for 17.1.1.1 (ES0-DMD) 01/25/02 17:54:00 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 @(#)lnsearch.for 17.1.1.1 (ESO) 01/25/02 17:54:00 C +++++++++++++++++++++++++++++++++++++++++++++++++++++++ C .COPYRIGHT (C) 1993 European Southern Observatory C .IDENT lnsearch.for C .AUTHOR: Daniel PONZ C C .KEYWORDS Spectroscopy, Long-Slit C C.PURPOSE: C SEARCH FOR EMISSION/ABSORPTION LINES IN A 1D- OR 2D-IMAGE AND C CALCULATE PARAMETERS : C - LINE CENTER C - PEAK VALUE C C.ALGORITHM: C C LINES ARE DETECTED BY USING A MEDIAN + THRESHOLDING ALGORITHM C ABSORPTION LINES ARE TREATED AS NEGATIVE EMISSION C CENTER OF THE LINE IS COMPUTED AS C (1)- CENTER OF GRAVITY OF THE 2 HIGHEST PIXELS RELATIVE TO THE C THIRD ONE C (2)- CENTER OF THE GAUSSIAN FIT TO THE POINTS ABOVE THE BACKGROUND C (3)- PIXEL CORRESPONDING TO THE MAXIMUM C C.INPUT/OUTPUT: C THE FOLLOWING KEYWORDS SHOULD BE DEFINED : C P1 : IMAGE TO BE SEARCHED C INPUTR : SEARCH PARAMETERS C WINDOW, THRESHOLD, NAVER C P3 : TABLE FILE C P4 : SEARCH METHOD C GRAVITY - DEFAULT C GAUSSIAN C MAXIMUM C P5 : EMISSION/ABSORPTION C C .VERSION 1.0 Package Creation 17-MAR-1993 C ------------------------------------------------------- C PROGRAM LNSEARCH IMPLICIT NONE INTEGER MD, NMAX PARAMETER (MD=2) PARAMETER (NMAX=1000) INTEGER MADRID INTEGER I,IAV,IC,IMETH,IMODE,IS,ISTAT, ALGO INTEGER NAX,NOBJ,NAVER,ISIZE INTEGER*8 IPM,IPNTA INTEGER KUN,KNUL,TID,IMNO INTEGER NPIX(MD), STEP C CHARACTER*2 METH,IEA CHARACTER*80 INA CHARACTER*80 FOUTA CHARACTER*48 UNIT CHARACTER*72 INST CHARACTER*80 LINE,LINEM,LINEE C REAL RBUF(2,NMAX) REAL STRA(MD),STPA(MD),PARM(5) DOUBLE PRECISION DSTRA(2), DSTPA(2), DPAR(5) C CHARACTER*16 LABEL(5),TUNIT(5) CHARACTER*6 FORM(5) C INCLUDE 'MID_INCLUDE:ST_DEF.INC' COMMON /VMR/MADRID(1) INCLUDE 'MID_INCLUDE:ST_DAT.INC' C DATA LABEL(1)/'X'/,TUNIT(1)/'PIXEL'/,FORM(1)/'F10.2'/ DATA LABEL(2)/'Y'/,TUNIT(2)/'PIXEL'/,FORM(2)/'F10.2'/ DATA LABEL(3)/'PEAK'/,TUNIT(3)/'DN'/,FORM(3)/'E12.3'/ C C ... INITIATE MIDAS C CALL STSPRO('LNSEARCH') C CALL STKRDC('IN_A',1,1,80,IAV,INA,KUN,KNUL,IS) CALL STKRDC('OUT_A',1,1,80,IAV,FOUTA,KUN,KNUL,IS) CALL STKRDC('IN_B',1,1,2,IAV,METH,KUN,KNUL,IS) CALL STKRDC('OUT_B',1,1,1,IAV,IEA,KUN,KNUL,IS) CALL STKRDI('INPUTI',1,1,IAV,ALGO,KUN,KNUL,IS) CALL STKRDR('INPUTR',1,4,IAV,PARM,KUN,KNUL,IS) CALL STKRDD('INPUTD',1,4,IAV,DPAR,KUN,KNUL,IS) CALL FORUPC(METH,METH) CALL FORUPC(IEA,IEA) C CALL STIGET(INA,D_R4_FORMAT,F_I_MODE,F_IMA_TYPE, . MD,NAX,NPIX,DSTRA,DSTPA,INST,UNIT,IPNTA,IMNO,IS) C NPIX(1) = MAX(1,NPIX(1)) NPIX(2) = MAX(1,NPIX(2)) NAX = 2 C STRA(1) = DSTRA(1) C STRA(2) = DSTRA(2) C STPA(1) = DSTPA(1) C STPA(2) = DSTPA(2) STRA(1) = DPAR(1) STRA(2) = DPAR(2) STPA(1) = DPAR(3) STPA(2) = DPAR(4) IF (NPIX(1).EQ.1) THEN NPIX(1) = NPIX(2) NPIX(2) = 1 END IF C C ... READ PARAMETERS C C NMAX = 800 C IF (METH(1:2).EQ.'GA') THEN IMETH = -1 LINEM = ' centering method : gaussian fit' ELSE IF (METH(1:2).EQ.'GR') THEN IMETH = 0 LINEM = ' centering method : center of gravity' ELSE IMETH = 1 LINEM = ' centering method : position of max(min)' END IF END IF IF (IEA(1:1).EQ.'E') THEN IMODE = 1 LINEE = ' search for emission lines' ELSE IMODE = 0 LINEE = ' search for absorption lines' END IF C C ... CREATE THE TABLE FILE PLUS COLUMNS C CALL TBTINI(FOUTA,F_TRANS,F_O_MODE,15,NMAX,TID,IS) CALL STDWRD(TID,'PIXEL',DSTPA,1,1,KUN,IS) DO 10 I = 1,3 CALL TBCINI(TID,D_R4_FORMAT,1,FORM(I), . TUNIT(I),LABEL(I),IC,IS) 10 CONTINUE CALL STTPUT(' search lines ',ISTAT) CALL STTPUT(' ------------',ISTAT) CALL STTPUT(' input image : '//INA,ISTAT) CALL STTPUT(' output table : '//FOUTA,ISTAT) CALL STTPUT(' input parameters ',ISTAT) WRITE (LINE,9000) NINT(PARM(1)) CALL STTPUT(LINE,ISTAT) WRITE (LINE,9010) PARM(2) CALL STTPUT(LINE,ISTAT) NAVER = PARM(3) STEP = PARM(4) WRITE (LINE,9030) NAVER CALL STTPUT(LINE,ISTAT) WRITE (LINE,9040) STEP CALL STTPUT(LINE,ISTAT) CALL STTPUT(LINEM,ISTAT) CALL STTPUT(LINEE,ISTAT) C C ... DO THE SEARCH C ISIZE = NPIX(1)*4 CALL TDMGET(ISIZE,IPM,ISTAT) CALL SRCOD3(MADRID(IPNTA),NAX,NPIX(1),NPIX(2),STRA, + STPA,NMAX,PARM,NOBJ,TID,RBUF,IMETH,IMODE,NAVER,STEP, + MADRID(IPM),ALGO) CALL TDMFRE(ISIZE,IPM,ISTAT) C WRITE (LINE,9020) NOBJ CALL STTPUT(LINE,ISTAT) C C ... FINISHED C CALL TBSINI(TID,IS) CALL TBTCLO(TID,IS) CALL STSEPI 9000 FORMAT (' search window : ',I3,' pixels') 9010 FORMAT (' detection threshold : ',F10.2,' DN') 9030 FORMAT (' average on : ',I3,' scan-lines') 9040 FORMAT (' step of : ',I3,' scan-lines') 9020 FORMAT (' no. of detections: ',I8) END C-------------------------------------------------------------------- SUBROUTINE FNDMAX(X,N,R,NMAX,NACT,START,STEP) C C FIND MAXIMUM VALUE IN CASE OF HIGH THRESHOLD C IMPLICIT NONE INTEGER N, NMAX, NACT REAL X(N),R(2,NMAX) REAL START, STEP C INTEGER I, IP REAL XMAX C XMAX = 1.E-30 DO 10 I = 3,N - 2 IF (X(I).GT.XMAX) THEN XMAX = X(I) IP = I END IF 10 CONTINUE R(1,1) = START + (IP-1)*STEP R(2,1) = XMAX NACT = 1 RETURN END C-------------------------------------------------------------------- SUBROUTINE SRCOD3(F,NAF,NPF1,NPF2,STRA,STPA,NMAX,PARM,NOBJ, + TID,RBUF,IMETH,IMODE,NAVER,STEP,RVAL,ALGO) C C C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ C C.IDENTIFICATION C SUBROUTINE SRCOD3 VERSION 1 23 NOV 1982 C J.D.PONZ C C.KEYWORDS: C SEARCH C C.PURPOSE: C SEARCH SPECTRA FOR EMISSION LINES C C.ALGORITHM: C USE FINDM2 C C.INPUT/OUTPUT C F original image Input C NAF dimension of F Input C NPF number of pixels in F (x,y,z...) Input C STRA starting point in F Input C STPA step Input C PARM floating-point parameters of filter Input C NOBJ number of objects found Output C OBJ charateristics of objects Output C C------------------------------------------------------------- C C IMPLICIT NONE INTEGER NPF1, NPF2, NAF, NMAX, NOBJ, IMETH, IMODE, ALGO INTEGER NAVER, STEP REAL STRA(2),STPA(2),PARM(2),OBJ(6) REAL RBUF(2,NMAX), F(NPF1,NPF2), RVAL(NPF1) REAL XAUX(100) INTEGER IAUX(100), I, NSEQ, IS, NACT, TID INTEGER J, J1, J2, K INTEGER ICOL(5), NAUX, NT, NW, ISTAT, IAC, NW1, IORD REAL T CHARACTER LINE*80 C DATA ICOL/1,2,3,4,5/ C NAUX = 99 C NT = NPF1 NW = PARM(1) T = PARM(2) NOBJ = 0 C C search features C IF (NPF2.GT.1) THEN CALL STTPUT(' seq.no. detected lines ',ISTAT) CALL STTPUT(' ------- -------------- ',ISTAT) END IF IAC = 3 NW1 = 2*NW + 1 NSEQ = 0 DO 20 IORD = 1, NPF2, STEP DO 1 J = 1, NPF1 RVAL(J) = 0. 1 CONTINUE J1 = MAX(IORD-NAVER,1) J2 = MIN(J1 + NAVER,NPF2) DO 3 J = J1, J2 DO 2 K = 1, NPF1 RVAL(K) = RVAL(K) + F(K,J)/(NAVER*2+1) 2 CONTINUE 3 CONTINUE IF (ALGO .EQ. 0) THEN CALL FNDMDN(RVAL,NT,NW1,T,NMAX,RBUF,NACT,XAUX,IAUX, + NAUX, IMETH,IMODE,STRA(1),STPA(1)) ELSE CALL FINDIT(RVAL,NT,NW1,T,NMAX,RBUF,NACT,XAUX,IAUX, + NAUX, IMETH,IMODE,STRA(1),STPA(1)) ENDIF IF (NACT.LE.0) CALL FNDMAX(RVAL,NT,RBUF,NMAX,NACT, + STRA(1),STPA(1)) IF (NPF2.GT.1) THEN WRITE (LINE,9000) IORD,NACT CALL STTPUT(LINE,ISTAT) END IF NOBJ = NOBJ + NACT DO 10 I = 1,NACT OBJ(1) = RBUF(1,I) OBJ(2) = (IORD-1)*STPA(2) + STRA(2) OBJ(3) = RBUF(2,I) NSEQ = NSEQ + 1 CALL TBRWRR(TID,NSEQ,IAC,ICOL,OBJ,IS) 10 CONTINUE 20 CONTINUE IF (NPF2.GT.1) CALL STTPUT(' ----------------------- ',ISTAT) RETURN C 9000 FORMAT (1X,I8,2X,I8) END C----------------------------------------------------------------------- SUBROUTINE FNDMDN(IA,LEN,WIDTH,THRES,NMAX,RBUF,NACT,FA,LINK, + NAUX,IMETH,IMODE,XSTR,XSTP) C C.COPYRIGHT: Copyright (c) 1987 European Southern Observatory, C all rights reserved C C.VERSION: 1.0 ESO-FORTRAN Conversion, AA 23:58 - 3 DEC 1987 C C.LANGUAGE: F77+ESOext C C.AUTHOR: D.PONZ C C.IDENTIFICATION C C SUBROUTINE FNDMDN C C.KEYWORDS C C SEARCH C C.PURPOSE C C SEARCH EMISSION/ABSORPTION FEATURES C C.ALGORITHM C C DATA IN A RUNNING WINDOW IS SORTED TO FORM A LOCAL HISTOGRAM C AN EMISSION(ABSORPTION) FEATURE IS FOUND WHEN THE CENTRAL C VALUE IN THE WINDOW IS GREATER THAN THE MEDIAN VALUE PLUS THE C THRESHOLD. C C.INPUT/OUTPUT C 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 NMAX INTG*4 MAXIMUM NO. OF DETECTED LINES C RBUF REAL*4 RBUF(1,I) - CENTER C RBUF(2,I) - MAXIMUM VALUE C NACT INTG*4 ACTUAL NUMBER OF DETECTED LINES 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 IMETH INTG*4 CENTERING METHOD AS C IMODE INTG*4 0 C XSTR REAL*4 STARTING COORD. C XSTP REAL*4 STEP C C- C IMPLICIT NONE INTEGER NAUX INTEGER LINK(NAUX),MIN,MAX,WIDTH,LEN,MIDPT,FRONT,AVAIL INTEGER RANGE,NEWNOD,I,J,K,NODE,ITMP INTEGER NMAX,NACT,IMODE,IMETH INTEGER ISTART,IFAIL,NPIX,NOBJ,NLOBJ REAL FA(NAUX),IA(LEN),RBUF(2,NMAX),AM REAL THRES,T,XSTR,XSTP,CHICHK,XGO,XCNTR REAL ACOE(4),WORK1(200),WORK2(200),XPREV REAL DELTA LOGICAL IFIRST 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 CHICHK = 0.005 IFIRST = .TRUE. T = ABS(THRES) MIN = 2 MAX = NAUX AVAIL = MIN FA(1) = -1.0 LINK(1) = 0 NOBJ = 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 RANGE = WIDTH/2 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 MIDPT = RANGE + 1 C MIDPT IS THE MIDDLE OF THE WINDOW (IN FA) DO 70 J = WIDTH,LEN,1 NODE = NEWNOD(IA(J),AVAIL,MIN,MAX,FA,LINK,NAUX) CALL INSNOD(NODE,FA,LINK,NAUX) K = K + 1 C C FIND THE MIDPT - - IT IS THE MEDIAN VALUE C ITMP = LINK(1) DO 20 I = 2,MIDPT,1 ITMP = LINK(ITMP) 20 CONTINUE C C GET THE LINK WHICH POINTS TO THE VALUE OF THE MEDIAN C IF (ABS(FA(ITMP)-IA(K)).GT.T) THEN IF (IFIRST) THEN IFIRST = .FALSE. NLOBJ = 1 ISTART = K - 1 NPIX = 3 ELSE NPIX = NPIX + 1 END IF ELSE IF ( .NOT. IFIRST) THEN IFIRST = .TRUE. C C FIND LINE CENTER C XGO = XSTR + (ISTART-1)*XSTP IF (IMETH) 30,40,50 C C GAUSSIAN C 30 CALL SGAUS(IA(ISTART),WORK1,WORK2,IMODE,NPIX,IFAIL, + XGO,XSTP,XCNTR,CHICHK,1,AM,ACOE) GO TO 60 C C GRAVITY C 40 CALL GRAVT(IA(ISTART),NPIX,IMODE,IFAIL,XGO,XSTP,XCNTR, + AM) IF (IFAIL.EQ.0) GO TO 60 C C MAXIMUM C 50 CALL CNTRH(IA(ISTART),NPIX,IMODE,IFAIL,XGO,XSTP,XCNTR, + AM) 60 IF (IFAIL.EQ.0) THEN NOBJ = NOBJ + 1 RBUF(1,NOBJ) = XCNTR RBUF(2,NOBJ) = AM IF (NOBJ.EQ.NMAX) GO TO 80 END IF END IF END IF CALL DELNOD(FRONT,FA,LINK,NAUX) FRONT = FRONT + 1 IF (FRONT.GT.MAX) FRONT = MIN C C DELETE OLDEST AND ADVANCE PTR, FRONT ALSO WRAPS C 70 CONTINUE 80 CONTINUE C C CHECK DOUBLETS C XPREV = RBUF(1,1) I = 1 85 CONTINUE IF ( .NOT. (I.LT.NOBJ)) GO TO 110 DELTA = 2.*XSTP IF (ABS(XPREV-RBUF(1,I+1)).LE.DELTA) THEN DO 90 J = I,NOBJ - 2 RBUF(1,J) = RBUF(1,J+2) RBUF(2,J) = RBUF(2,J+2) 90 CONTINUE NOBJ = NOBJ - 2 ELSE I = I + 1 END IF XPREV = RBUF(1,I) GOTO 85 110 CONTINUE NACT = NOBJ RETURN END C----------------------------------------------------------------------- SUBROUTINE FINDIT(IA,LEN,WIDTH,THRES,NMAX,RBUF,NACT,FA,LINK, + NAUX,IMETH,IMODE,XSTR,XSTP) C+ C.COPYRIGHT: Copyright (c) 1993 European Southern Observatory, C all rights reserved C C.AUTHOR: P. BALLESTER C C.IDENTIFICATION C C SUBROUTINE FINDIT C C.INPUT/OUTPUT C 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 NMAX INTG*4 MAXIMUM NO. OF DETECTED LINES C RBUF REAL*4 RBUF(1,I) - CENTER C RBUF(2,I) - MAXIMUM VALUE C NACT INTG*4 ACTUAL NUMBER OF DETECTED LINES 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 IMETH INTG*4 CENTERING METHOD AS C IMODE INTG*4 0 C XSTR REAL*4 STARTING COORD. C XSTP REAL*4 STEP C C- IMPLICIT NONE INTEGER NAUX INTEGER LINK(NAUX),MIN,MAX,WIDTH,LEN INTEGER RANGE,I,J,K,NODE INTEGER NMAX,NACT,IMODE,IMETH INTEGER ISTART,IFAIL,NPIX,NOBJ REAL FA(NAUX),IA(LEN),RBUF(2,NMAX),AM REAL THRES,T,XSTR,XSTP,CHICHK,XGO,XCNTR REAL ACOE(4),WORK1(200),WORK2(200),XPREV REAL DELTA LOGICAL IFIRST REAL R(2,1) INTEGER IH,PREVK 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 CHICHK = 0.005 IFIRST = .TRUE. T = ABS(THRES) MIN = 2 MAX = NAUX FA(1) = -1.0 NOBJ = 0 PREVK = 0 C C SET UP FIRST FILTER WINDOW C RANGE = WIDTH - 1 RANGE = WIDTH/2 C C FOR EACH POINT FROM WIDTH/2 TO LEN-WIDTH/2+1 C FIND ITS MEDIAN C DO 70 J = RANGE+1,LEN-RANGE C C GET THE LINK WHICH POINTS TO THE VALUE OF THE MEDIAN C MIN = J-RANGE MAX = J+RANGE DO 66 NODE = MIN,MAX FA(NODE-MIN+1) = IA(NODE) 66 CONTINUE CALL FNDMAX(FA,WIDTH,R,1,IH,1.,1.) IH = R(1,1) K = IH+MIN-1 C C FIND LINE CENTER C IF (K.NE.PREVK.AND.IA(K).GT.T) THEN PREVK = K ISTART = K - RANGE NPIX = WIDTH XGO = XSTR + (ISTART-1)*XSTP IF (IMETH) 30,40,50 C C GAUSSIAN C 30 CALL SGAUS(IA(ISTART),WORK1,WORK2,IMODE,NPIX,IFAIL, + XGO,XSTP,XCNTR,CHICHK,1,AM,ACOE) GO TO 60 C C GRAVITY C 40 CALL GRAVT(IA(ISTART),NPIX,IMODE,IFAIL,XGO,XSTP,XCNTR, + AM) IF (IFAIL.EQ.0) GO TO 60 C C MAXIMUM C 50 CALL CNTRH(IA(ISTART),NPIX,IMODE,IFAIL,XGO,XSTP,XCNTR, + AM) C 60 IF (IFAIL.EQ.0) THEN NOBJ = NOBJ + 1 RBUF(1,NOBJ) = XCNTR RBUF(2,NOBJ) = AM IF (NOBJ.EQ.NMAX) GO TO 80 END IF ENDIF 70 CONTINUE 80 CONTINUE C C CHECK DOUBLETS C XPREV = RBUF(1,1) I = 1 85 CONTINUE IF ( .NOT. (I.LT.NOBJ)) GO TO 110 DELTA = WIDTH*XSTP IF (ABS(XPREV-RBUF(1,I+1)).LE.DELTA) THEN DO 90 J = I,NOBJ - 1 RBUF(1,J) = RBUF(1,J+1) RBUF(2,J) = RBUF(2,J+1) 90 CONTINUE NOBJ = NOBJ - 1 ELSE I = I + 1 END IF XPREV = RBUF(1,I) GOTO 85 110 CONTINUE NACT = NOBJ RETURN END