C @(#)necnord.for 17.1.1.1 (ESO-DMD) 01/25/02 17:51:30 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 ECHNORD C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ C.COPYRIGHT: Copyright (c) 1987 European Southern Observatory, C all rights reserved C C.VERSION: 1.0 ESO-FORTRAN Conversion, AA 21:22 - 3 DEC 1987 C C.LANGUAGE: F77+ESOext C C.AUTHOR: D.PONZ C C.MODIFICATIONS: M.PERON C PB 910313 Suppress creation of TRACE.bdf C C.IDENTIFICATION C echnord.for C C.KEYWORDS C C echelle, order searching, maxima following C C.PURPOSE C C Search echelle orders on a flat field image. C Output is a table with the following columns : C 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. 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(2) T THRESHOLD IN THE DETECTION C C- IMPLICIT NONE INTEGER MADRID INTEGER ACTVAL,I,ISTAT,IXSTR INTEGER NACOL,NAXIS,NDIM,NORDER,NP INTEGER*8 ITRACE INTEGER NTOT INTEGER*8 PNTRA INTEGER STATUS,WINDOW,IACOL,ID INTEGER NPIX(3),ICOL(7), KUN, KNUL, IMNO C INTEGER*8 IMSK CC INTEGER II INTEGER IJ INTEGER UPPER(200),LOWER(200),INEW REAL RPAR(4),LEVEL,FORD(200),FBACK(200),SLOPE DOUBLE PRECISION START(3),STEP(3) DOUBLE PRECISION ST CHARACTER*80 CA CHARACTER*80 TABLE,TRACE CHARACTER*72 IDENA CHARACTER*80 LINE CHARACTER*64 CUNITA CHARACTER*16 LABCOL(5),UNIT(5) CHARACTER*6 FORM(5) CHARACTER*60 FRMIN INCLUDE 'MID_INCLUDE:ST_DEF.INC' COMMON /VMR/MADRID(1) INCLUDE 'MID_INCLUDE:ST_DAT.INC' 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)/'World Coord'/ DATA UNIT(3)/'World Coord'/ DATA UNIT(4)/'World Coord'/ DATA UNIT(5)/'DN'/ DATA NACOL/5/ DATA TRACE/'TRACE.bdf'/ CA = 'TEST' ST = 1.D0 NDIM = 3 C C ... initialize system and read parameters C CALL STSPRO('ECHNORD') CALL STKRDC('IN_A',1,1,60,ACTVAL,FRMIN,KUN,KNUL,STATUS) CALL CLNFRA(FRMIN,FRMIN,0) CALL STKRDC('IN_B',1,1,60,ACTVAL,TABLE,KUN,KNUL,STATUS) CALL CLNTAB(TABLE,TABLE,0) CALL STKRDR('INPUTR',1,3,ACTVAL,RPAR,KUN,KNUL,STATUS) WINDOW = RPAR(1) LEVEL = RPAR(2) SLOPE = RPAR(3) 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) WINDOW CALL STTPUT(LINE,ISTAT) WRITE (LINE,9010) RPAR(2) 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,IMNO,STATUS) C C ......................... CREATE ORDER MASK FOR DEBUGGING C C CALL STIPUT('MASK',D_R4_FORMAT,F_O_MODE,F_IMA_TYPE, C . 1,NPIX(2),ST,ST,CA,CA,IMSK,II,STATUS) CALL STIPUT(TRACE,D_R4_FORMAT,F_O_MODE,F_IMA_TYPE, . 1,NPIX(2),ST,ST,CA,CA,ITRACE,IJ,STATUS) C C ... generate trace C IXSTR = NPIX(1)/2 CALL EXTR3(MADRID(PNTRA),NPIX(1),NPIX(2),IXSTR, . MADRID(ITRACE)) C C ... find order positions in the middle of the image by simple thresholding C NORDER = 0 NP = NPIX(2) - 2*WINDOW CALL ABSTHR(MADRID(ITRACE),NPIX(2),LEVEL, . NORDER,UPPER,LOWER) CALL ORDBAC(NPIX(2),MADRID(ITRACE), . NORDER,UPPER,LOWER,FORD,FBACK) CALL ORDREJ(WINDOW,NORDER,UPPER,LOWER,FORD,FBACK) C C ... while new orders are found detect previous and next orders C using variable thresholding C 50 CONTINUE CALL LVRTHR(MADRID(ITRACE),NPIX(2), . NORDER,UPPER,LOWER,FORD,FBACK,INEW) IF (INEW.GT.0) GOTO 50 60 CONTINUE CALL RVRTHR(MADRID(ITRACE),NPIX(2), . NORDER,UPPER,LOWER,FORD,FBACK,INEW) IF (INEW.GT.0) GOTO 60 C C ... generate mask for control C C CALL MASKIN(MADRID(IMSK),NPIX(2),NORDER,UPPER,LOWER,FORD,FBACK) C IF (NORDER.EQ.0) THEN CALL STTPUT('Ordr detection failed',STATUS) CALL STTPUT('Use a lower threshold',STATUS) GO TO 20 END IF WRITE (LINE,9030) NORDER CALL STTPUT(LINE,ISTAT) C C ... allocate table C NTOT = (NORDER+1)*NPIX(1)/10 IACOL = 12 CALL TBTINI(TABLE,F_TRANS,F_O_MODE,IACOL,NTOT,ID,STATUS) DO 10 I = 1,NACOL CALL TBCINI(ID,D_R4_FORMAT,1,FORM(I), . UNIT(I),LABCOL(I),ICOL(I),STATUS) 10 CONTINUE C C ... follow orders C CALL FOLLOW(MADRID(PNTRA),NPIX(1),NPIX(2),NORDER,UPPER, . LOWER,FORD,FBACK,ID,SLOPE) C C ... update number of elements C C NTOT = NPOINT CALL STKWRI('OUTPUTI',NORDER,1,1,KUN,STATUS) CALL STKWRI('OUTPUTI',NPIX,2,2,KUN,STATUS) IDENA = 'ORDER POSITION' CALL TBSINI(ID,STATUS) CALL TBTCLO(ID,STATUS) C C ... free data C 20 CONTINUE CALL STFDEL(TRACE,STATUS) CALL STSEPI 9000 FORMAT (' ORDER WIDTH : ',I6) 9010 FORMAT (' THRESHOLD : ',F10.1) 9030 FORMAT (' NUMBER OF DETECTED ORDERS : ',I4) END SUBROUTINE FOLLOW(X,N1,N2,NORDER,UPPER,LOWER,FO,FB,TID,SLOPE) C C DEFINE ORDER POSTIONS BY DETECTING THE UPPER AND LOWER C EDGES C IMPLICIT NONE C INTEGER N1 ! NUMBER OF SAMPLES INTEGER N2 ! NUMBER OF ROWS REAL X(N1,N2) ! INPUT IMAGE INTEGER NORDER ! NUMBER OF DETECTED ORDERS INTEGER UPPER(*) ! ROW NUMBER OF ASCENDING EDGE INTEGER LOWER(*) ! ROW NUMBER OF DESCENDING EDGE REAL FO(*) ! ORDER FLUX IN THE MIDDLE OF THE IMAGE REAL FB(*) ! BACKGROUND IN THE MIDDLE OF THE IMAGE INTEGER TID ! TABLE IDENTIFIER C INTEGER ISTAT, ICOL(5), IXSTR, NP, IX, IL, IU, NC, W,IP INTEGER IORD, WINDOW, ILIM1, ILIM2, IHW, IBACK, K, IBUF INTEGER LOWBND, UPPBND INTEGER ILPOS, IUPOS, IRL, ISTEP, NSAMP, ORDPOS, NF1, NF2 C REAL THRES, T1, T2, BACK, RBACK1, RBACK2, RVAL(5) REAL F1, F2, BUF(100, 4),SLOPE C CHARACTER*80 LINE C LOGICAL NEXT C DATA ICOL/1,2,3,4,5/ C NSAMP = 10 WINDOW = 10 IXSTR = N1/2 ILIM2 = N2 - WINDOW ILIM1 = N1 - WINDOW C LOWBND = 2 UPPBND = N2 - 2 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 NP = 0 DO 40 IORD = 1,NORDER IF (IORD.NE.NORDER) IHW = (UPPER(IORD+1)-LOWER(IORD))*0.5 W = (LOWER(IORD)-UPPER(IORD))*0.5 IF (IORD.EQ.1) THEN BACK = FB(1) ELSE IF (IORD.EQ.NORDER) THEN BACK = FB(NORDER-1) ELSE BACK = 0.5*(FB(IORD-1) + FB(IORD)) ENDIF THRES = BACK + 0.5*(FO(IORD) - BACK) C C ... LEFT/RIGHT PART OF THE IMAGE C NC = 3 IBUF = 0 DO 30 IRL = 1,2 C C ... INITIAL THRESHOLD C IF (IRL.EQ.1) THEN ISTEP = -1 ELSE C ... INCLUDE START POSITION OF THE ORDER NP = NP + 1 RVAL(1) = IORD RVAL(2) = 1. CALL TBRWRR(TID,NP,2,ICOL,RVAL,ISTAT) DO 5 K = IBUF, 1, -1 NP = NP + 1 NC = NC + 1 RVAL(1) = IORD RVAL(2) = BUF(K,1) RVAL(3) = BUF(K,2) IF (BUF(K,3).LT.0.) THEN CALL TBRWRR(TID,NP,3,ICOL,RVAL,ISTAT) ELSE RVAL(4) = BUF(K,3) RVAL(5) = BUF(K,4) CALL TBRWRR(TID,NP,5,ICOL,RVAL,ISTAT) ENDIF 5 CONTINUE C ... INCLUDE CENTRAL POSITION OF THE ORDER NP = NP + 1 RVAL(1) = IORD RVAL(2) = IXSTR RVAL(3) = UPPER(IORD) + (LOWER(IORD)-UPPER(IORD))/2 ORDPOS = RVAL(3) IF (IORD.LT.NORDER) THEN IBACK = LOWER(IORD) + IHW IF (IBACK.LT.ILIM2) THEN RVAL(4) = IBACK RVAL(5) = BACK CALL TBRWRR(TID,NP,5,ICOL,RVAL,ISTAT) ENDIF ELSE CALL TBRWRR(TID,NP,3,ICOL,RVAL,ISTAT) END IF ISTEP = 1 END IF T1 = THRES T2 = THRES IX = IXSTR IU = UPPER(IORD) IL = LOWER(IORD) NEXT = .TRUE. IP = 1 10 CONTINUE IX = IX + ISTEP C C ... FIND ASCENDING EDGE C F1 = 0. NF1 = 0 IUPOS = 0 RBACK1 = 10000. DO 50 K = IU-1-IP*SLOPE, IU+W+IP*SLOPE IF (K.GE.LOWBND) THEN RBACK1 = MIN(RBACK1, X(IX, K)) IF (X(IX,K-1).LT.T1.AND.X(IX,K).GT.T1) IUPOS = K IF (IUPOS.GT.0) THEN F1 = F1 + X(IX,K) NF1 = NF1 + 1 ENDIF ENDIF 50 CONTINUE C C ... FIND DESCENDING EDGE C ILPOS = 0 F2 = 0. NF2 = 0 RBACK2 = 10000. DO 70 K = IL+1+IP*SLOPE, IL-W-IP*SLOPE, -1 IF (K.LE.UPPBND) THEN RBACK2 = MIN(RBACK2, X(IX, K)) IF (X(IX,K+1).LT.T2.AND.X(IX,K).GT.T2) ILPOS = K IF (ILPOS.GT.0) THEN F2 = F2 + X(IX,K) NF2 = NF2 + 1 ENDIF ENDIF 70 CONTINUE IF ((IUPOS.GT.0).AND.(ILPOS.GT.0)) THEN T1 = RBACK1 + 0.5*(F1/NF1 - RBACK1) T2 = RBACK2 + 0.5*(F2/NF2 - RBACK2) C C ... STORE INTO TABLE C IF (MOD(IX,NSAMP).EQ.0) THEN IF (IRL.EQ.1) THEN ! STORE ON BUFFER IF (IORD.LT.NORDER) THEN IBACK = ILPOS + IHW IF (IBUF.GE.200) GOTO 30 IF (IBACK.LT.ILIM2) THEN IBUF = IBUF + 1 BUF(IBUF,1) = IX BUF(IBUF,2) = IUPOS+(ILPOS-IUPOS)/2 BUF(IBUF,3) = IBACK BUF(IBUF,4) = 0.5*(RBACK1 + RBACK2) ENDIF ELSE IBUF = IBUF + 1 BUF(IBUF,1) = IX BUF(IBUF,2) = IUPOS+(ILPOS-IUPOS)/2 BUF(IBUF,3) = -1. BUF(IBUF,4) = -1. END IF ELSE NC = NC + 1 NP = NP + 1 RVAL(1) = IORD RVAL(2) = IX RVAL(3) = IUPOS+(ILPOS-IUPOS)/2 IF (IORD.LT.NORDER) THEN IBACK = ILPOS + IHW IF (IBACK.LT.ILIM2) THEN RVAL(4) = IBACK RVAL(5) = 0.5*(RBACK1 + RBACK2) CALL TBRWRR(TID,NP,5,ICOL,RVAL,ISTAT) ENDIF ELSE CALL TBRWRR(TID,NP,3,ICOL,RVAL,ISTAT) END IF END IF END IF ENDIF IF (ILPOS.GT.0.AND.IUPOS.GT.0) THEN IL = ILPOS IU = IUPOS IP = 1 ELSE IP=IP+1 ENDIF C NEXT = (IU.LT.IL).AND.(IL.LT.ILIM2).AND.(IU.GT.IHW) C . .AND.(IX.GT.WINDOW).AND.(IX.LT.ILIM1) NEXT = (IU.LT.IL).AND.(IL.LT.ILIM2).AND.(IU.GT.IHW).AND. . (IP.LT.11).AND.(IX.GT.WINDOW).AND.(IX.LT.ILIM1) IF (NEXT) GO TO 10 C 20 CONTINUE 30 CONTINUE C ... INCLUDE END POSITION OF THE ORDER NP = NP + 1 RVAL(1) = IORD RVAL(2) = N1 CALL TBRWRR(TID,NP,2,ICOL,RVAL,ISTAT) C WRITE (LINE,9010) IORD,IXSTR,ORDPOS,IHW,NC,T1,T2 WRITE (LINE,9010) IORD,IXSTR,ORDPOS,IHW,NC CALL STTPUT(LINE,ISTAT) 40 CONTINUE WRITE (LINE,9040) CALL STTPUT(LINE,ISTAT) RETURN C 9000 FORMAT (I4) C 9010 FORMAT (1X,I7,2X,I8,2X,I8,2X,I10,2X,I8) 9010 FORMAT (1X,I7,2X,I8,2X,I8,2X,I10,2X,I8,2F10.1) 9020 FORMAT (' SEQ.NO. X CENTER Y CENTER INTERORDER TEMPLATE') 9030 FORMAT (' ------- -------- -------- ---------- --------') 9040 FORMAT (' -------------------------------------------------') END SUBROUTINE ABSTHR(X,N,THRES,NORDER,UPPER,LOWER) C C FIND EDGES OF THE ORDERS BY SIMPLE THRESHOLDING C IMPLICIT NONE C INTEGER N ! NUMBER OF SAMPLES REAL X(N) ! INPUT IMAGE TRACE REAL THRES ! DETECTING THRESHOLD INTEGER NORDER ! NUMBER OF DETECTED ORDERS INTEGER UPPER(*) ! ROW NUMBER OF ASCENDING EDGE INTEGER LOWER(*) ! ROW NUMBER OF DESCENDING EDGE C INTEGER I, NO LOGICAL UP C C ... ITERATION C UP = .TRUE. NO = 0 C DO 10 I = 1, N IF (UP) THEN IF (X(I).GT.THRES .AND. X(I-1).LT.THRES) THEN C ... ASCENDING EDGE UP = .FALSE. NO = NO + 1 UPPER(NO) = I ENDIF ELSE IF (X(I).LT.THRES .AND. X(I-1).GT.THRES) THEN C ... DESCENDING EDGE UP = .TRUE. LOWER(NO) = I-1 ENDIF ENDIF 10 CONTINUE IF (UP) THEN NORDER = NO ELSE NORDER = NO - 1 ENDIF RETURN END SUBROUTINE MASKIN(X,N,NO,UPPER,LOWER,FO,FB) C C DERIVE A MASK WITH THE ORDER POSITIONS C IMPLICIT NONE C INTEGER N ! NUMBER OF PIXELS REAL X(N) ! MASK INTEGER NO ! INPUT/OUTPUT NUMBER OF ORDERS INTEGER UPPER(*) ! ASCENDING EDGE POSITION INTEGER LOWER(*) ! DESCENDING EDGE POSITION REAL FO(*) REAL FB(*) C INTEGER I, J REAL ON, OFF, THRES, BACK C ON = 10000. OFF = 0. THRES = ON/2. C C ... BUILD THE MASK C BACK = 0. J = 1 DO 10 I = 1, N IF (I .GE. UPPER(J) .AND. I .LE. LOWER(J) ) THEN X(I) = FO(J) ELSE IF (I.GT.LOWER(J)) THEN BACK = FB(J) J = J + 1 IF (J.GT.NO) GOTO 20 ENDIF X(I) = BACK ENDIF 10 CONTINUE RETURN 20 DO 30 J = I, N X(J) = BACK 30 CONTINUE C RETURN END SUBROUTINE EXTR3(X,N1,N2,I,OUT) C C EXTRACE PERPENDICULAR TRACE C INTEGER N1 INTEGER N2 REAL X(N1, N2) INTEGER I REAL OUT(N2) C INTEGER J C DO 10 J = 1, N2 OUT(J) = (X(I-1,J)+X(I,J)+X(I+1,J))/3. 10 CONTINUE RETURN END SUBROUTINE ORDBAC(N,X, . NORDER,UPPER,LOWER,FORD,FBACK) IMPLICIT NONE C INTEGER N ! NUMBER OF SAMPLES REAL X(N) ! INPUT IMAGE INTEGER NORDER ! INPUT/OUTPUT NUMBER OF ORDERS INTEGER UPPER(*) ! ASCENDING EDGE INTEGER LOWER(*) ! DESCENDING EDGE REAL FORD(*) REAL FBACK(*) C INTEGER I, I1, I2, NBACK, IORD, IO1, IO2, NORD REAL ORDER, BACK C C ... COMPUTE FOR EACH ORDER THE AVERAGE FLUX AND THE NEXT BACKGROUND C DO 40 IORD = 1, NORDER IO1 = UPPER(IORD) + 1 !.............. IO2 = LOWER(IORD) - 1 !.............. NORD = IO2-IO1+1 ORDER = 0. DO 20 I = IO1,IO2 ORDER = ORDER + X(I)/NORD 20 CONTINUE FORD(IORD) = ORDER I1 = LOWER(IORD) + 1 IF (IORD.EQ.NORDER) THEN FBACK(IORD) = 0. ELSE I2 = UPPER(IORD+1) - 1 NBACK = I2-I1+1 C BACK = 0. BACK = 10000. !.......... DO 30 I = I1, I2 BACK = MIN(BACK, X(I)) !.......... C BACK = BACK + X(I)/NBACK 30 CONTINUE FBACK(IORD) = BACK ENDIF 40 CONTINUE RETURN END SUBROUTINE ORDREJ(WINDOW,NORDER,UPPER,LOWER,FORD,FBACK) C C REJECT SPURIOUS ORDERS IMPLICIT NONE C INTEGER WINDOW ! ORDER WIDTH INTEGER NORDER ! INPUT/OUTPUT NUMBER OF ORDERS INTEGER UPPER(*) ! ASCENDING EDGE INTEGER LOWER(*) ! DESCENDING EDGE REAL FORD(*) REAL FBACK(*) C INTEGER I, J, IW, STATUS, KK CHARACTER*70 LINE C I = 1 KK = NORDER 10 CONTINUE IW = LOWER(I)-UPPER(I) + 1 IF (ABS(IW-WINDOW).GT.2)THEN WRITE(LINE,1000)I,IW CALL STTPUT(LINE,STATUS) IF (I.NE.KK) THEN DO 20 J = I, KK-1 LOWER(J) = LOWER(J+1) UPPER(J) = UPPER(J+1) FORD(J) = FORD(J+1) FBACK(J) = FBACK(J+1) 20 CONTINUE ENDIF LOWER(KK) = 0 UPPER(KK) = 0 FORD(KK) = 0. FBACK(KK) = 0. KK = KK-1 I = I-1 ENDIF I = I+1 IF (I.LE.KK) GO TO 10 NORDER = KK RETURN 1000 FORMAT(' Reject order ',I4, ' computed width = ',I4) END SUBROUTINE LVRTHR(X,N,NORDER,UPPER,LOWER,FORD,FBACK,INEW) C C FIND EDGES OF THE ORDERS BY RELATIVE THRESHOLDING C IMPLICIT NONE C INTEGER N ! NUMBER OF SAMPLES REAL X(N) ! INPUT IMAGE TRACE INTEGER NORDER ! NUMBER OF DETECTED ORDERS INTEGER UPPER(*) ! ROW NUMBER OF ASCENDING EDGE INTEGER LOWER(*) ! ROW NUMBER OF DESCENDING EDGE REAL FORD(*) REAL FBACK(*) INTEGER INEW C INTEGER I, IUP, ILOW, N1, I1, I2 LOGICAL UP REAL THRES, F C INEW = 0 THRES = FBACK(1) + (FORD(1) - FBACK(1))/2. C C ... ITERATION C UP = .TRUE. C I1 = UPPER(1) - 2 DO 10 I = I1, 1, -1 IF (UP) THEN IF (X(I).GT.THRES .AND. X(I+1).LT.THRES) THEN C ... ASCENDING EDGE UP = .FALSE. IUP = I ENDIF ELSE IF (X(I).LT.THRES .AND. X(I+1).GT.THRES) THEN C ... DESCENDING EDGE UP = .TRUE. ILOW = I-1 INEW = 1 GO TO 20 ENDIF ENDIF 10 CONTINUE RETURN 20 CONTINUE C C ... INCLUDE NEW ORDER C NORDER = NORDER + 1 DO 30 I = NORDER, 2, -1 UPPER(I) = UPPER(I-1) LOWER(I) = LOWER(I-1) FORD(I) = FORD(I-1) FBACK(I) = FBACK(I-1) 30 CONTINUE UPPER(1) = ILOW LOWER(1) = IUP C C ... COMPUTE NEW BACKGROUND AND ORDER FLUX C F = 0. C I1 = UPPER(1) + 1 !......... I2 = LOWER(1) - 1 !......... N1 = I2 - I1 + 1 DO 40 I = I1, I2 F = F + X(I)/N1 40 CONTINUE FORD(1) = F F = 10000. I1 = UPPER(2)-1 I2 = IUP + 1 N1 = I1-I2 + 1 DO 50 I = I2, I1 F = MIN(F,X(I)) 50 CONTINUE FBACK(1) = F RETURN END SUBROUTINE RVRTHR(X,N,NORDER,UPPER,LOWER,FORD,FBACK,INEW) C C FIND EDGES OF THE ORDERS BY RELATIVE THRESHOLDING C IMPLICIT NONE C INTEGER N ! NUMBER OF SAMPLES REAL X(N) ! INPUT IMAGE TRACE INTEGER NORDER ! NUMBER OF DETECTED ORDERS INTEGER UPPER(*) ! ROW NUMBER OF ASCENDING EDGE INTEGER LOWER(*) ! ROW NUMBER OF DESCENDING EDGE REAL FORD(*) REAL FBACK(*) INTEGER INEW C INTEGER I, IUP, ILOW, N1, I1, I2 LOGICAL UP REAL THRES, F C ................................. NOTE THAT FBACK(NORDER) IS 0. INEW = 0 THRES = FBACK(NORDER-1)+(FORD(NORDER)-FBACK(NORDER-1))/2. C C ... ITERATION C UP = .TRUE. I1 = LOWER(NORDER) + 2 C DO 10 I = I1, N IF (UP) THEN IF (X(I).GT.THRES .AND. X(I-1).LT.THRES) THEN C ... ASCENDING EDGE UP = .FALSE. IUP = I ENDIF ELSE IF (X(I).LT.THRES .AND. X(I-1).GT.THRES) THEN C ... DESCENDING EDGE UP = .TRUE. ILOW = I-1 INEW = 1 GO TO 20 ENDIF ENDIF 10 CONTINUE RETURN 20 CONTINUE C C ... INCLUDE NEW ORDER C NORDER = NORDER + 1 UPPER(NORDER) = IUP LOWER(NORDER) = ILOW C C ... COMPUTE NEW BACKGROUND AND ORDER FLUX C F = 0. I1 = IUP + 1 I2 = ILOW -1 N1 = I2-I1+1 DO 40 I = I1, I2 F = F + X(I)/N1 40 CONTINUE FORD(NORDER) = F FBACK(NORDER) = 0. F = 10000. I2 = UPPER(NORDER)-1 I1 = LOWER(NORDER-1) + 1 N1 = I2-I1 + 1 DO 50 I = I1, I2 F = MIN(F,X(I)) 50 CONTINUE FBACK(NORDER-1) = F RETURN END