C @(#)search.for 17.1.1.1 (ES0-DMD) 01/25/02 17:15:45 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 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 C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ C.COPYRIGHT: Copyright (c) 1987 European Southern Observatory, C all rights reserved C.IDENTIFICATION: SEARCH.FOR C.PURPOSE: Searches image frame for objects and write found C objects into array ACAT. C C.ALGORITHM: Described in additional documentation. C C.LANGUAGE: ESO-FOR C.AUTHOR: A. Kruszewski C.KEYWORDS: GALAXIES, IMAGES, SEARCH, STARS C.ENVIRONMENT: Portable MIDAS C.COMMENTS: Subroutines SEARCH and SEAREG used by program INVSEARCH C.VERSION: 1.0 JUL 1981 Creation ESO-Garching C.VERSION: 1.1 16 AUG 1983 Modification ESO-Garching C.VERSION: 2.0 6 SEP 1983 Modified by Ch. Ounnas ESO-Garching C.VERSION: 22 MAY 1984 Modified by Ch. Ounnas ESO-Garching C.VERSION: 2.1 19 FEB 1987 Modified for FX Obs. de Geneve C.VERSION: 3.0 JUN 1987 Indirect pixel addresing and SEAREG C added ESO-Garching C.VERSION: 4.0 18 OCT 1988 Modified for portability Warsaw U. Obs. C.VERSION: 2 MAY 1989 Few minor changes ESO-Garching C----------------------------------------------------------------------- SUBROUTINE SEARCH(IMF, A, JAPY, NX, NY, & ACAT, IDET, IDA, IARR, RARR, MM) C IMPLICIT NONE INCLUDE 'MID_REL_INCL:INVENT.INC/NOLIST' C INTEGER IMF ! IN: Image file REAL A(1) ! LOC: Image data buffer INTEGER JAPY(1) ! LOC: Pointers to image lines INTEGER NX ! IN: Image X-dimension INTEGER NY ! IN: Image Y-dimension REAL ACAT(5,MAXCNT) ! OUT: Catalog of objects INTEGER IDET(1) ! LOC: Detection mapping array INTEGER IDA(1) ! LOC: Multiple detections limits INTEGER IARR(32) ! IN: Integer INVENTORY keywords REAL RARR(64) ! IN: Real INVENTORY keywords INTEGER MM ! OUT: Number of detected objects C INTEGER IBUF(4) INTEGER IHED, ILIM, ISTAT INTEGER INET(0:MAXNET) INTEGER ITMP INTEGER IUSD(4), IXYU(4) INTEGER JEND, JMSG, JS, JSRCHD INTEGER JSTART, K, M INTEGER NC, NH1, NHED, NL INTEGER MCAT(4,MAXDET) C INTEGER MMB INTEGER NHSEG, NLPB, NVSEG INTEGER NXS, NXU, NYS, NYU C REAL BCAT(2,MAXDET), CRMD(3), DLIM, SKYNET(2,0:MAXNET) C CHARACTER*80 TEXT C C *** Recall keywords. C NC = MAXCNT NL = MAXDET CRMD(1) = RARR(39) CRMD(2) = RARR(40) CRMD(3) = RARR(41) IHED = IARR(8) NHED = IARR(23) C C *** Define used part of frame IXYU. C IXYU(1) = MAX( 1 , IARR(12)-IHED ) IXYU(2) = MAX( 1 , IARR(13)-IHED ) IXYU(3) = MIN( NX , IARR(14)+IHED ) IXYU(4) = MIN( NY , IARR(15)+IHED ) C C *** Part of frame in buffer is limited by IBUF. C c IBUF(1) = IXYU(1) c IBUF(3) = IXYU(3) C C *** Presently buffer is identical with the frame. C IBUF(1) = 1 IBUF(2) = 1 IBUF(3) = NX IBUF(4) = NY C C *** Analysed part of frame shall be divided onto vertical C *** segments. Borders of actual segment are set by IUSD. C IUSD(1) = IXYU(1) IUSD(3) = IXYU(3) cC cC *** Find number of lines in buffer. cC c NPPL = IXYU(3) - IXYU(1) + 1 c NLPB = MIN ( NYBUF , NIBUF/NPPL ) NLPB = NY C C *** Searched part of frame is defined by IARR(12)-IARR(15) C NXU = IXYU(3) - IXYU(1) + 1 NYU = IXYU(4) - IXYU(2) + 1 NXS = IARR(14) - IARR(12) + 1 NYS = IARR(15) - IARR(13) + 1 C C *** Set array for combining detections. C DLIM = RARR(43) ILIM = INT( DLIM ) CALL LMTDET( ILIM , DLIM , IDA ) C C *** Initialize detection array. C ITMP = (ILIM+1) * NXS DO 40 K = 1 , ITMP IDET(K) = 0 40 CONTINUE C C *** Calculate number of vertical segments. C NVSEG = INT( RARR(48) * FLOAT(NYS) / FLOAT(2*NHED+1) ) + 1 cC cC *** Make segments small enough to fit buffer. cC c IF ( (NYS/NVSEG+1+2*IHED)*NXU .GT. NIBUF .OR. c & NYS/NVSEG+1+2*IHED .GT. NYBUF ) THEN c NLPB = NIBUF / NXU c IF ( NLPB-2*IHED .GT. IHED .AND. NLPB .LE. c & NYBUF ) THEN c NVSEG = NYS / (NLPB-2*IHED) + 1 c ELSE c CALL STTPUT('Check connection and include files',ISTAT) c RETURN c ENDIF c ENDIF C C *** Divide searched part of a frame into horizontal segments. C NHSEG = MIN ( INT( RARR(48)*FLOAT(NXS) / FLOAT(2*NHED+1) ) + 1 , & MAXNET ) INET(0) = IARR(12) INET(NHSEG) = IARR(14) NH1 = NHSEG-1 DO 10 K = 1 , NH1 INET(K) = INET(0) + ( K*NXS ) / NHSEG 10 CONTINUE C C *** Initialize counters. C M = 0 MM = 0 C C *** Start loop over vertical segments C c IBUF(2) = 0 c JBS = 0 c IBUF(4) = 0 c JBE = 0 JSRCHD = IARR(13) - 1 CALL STTPUT( 'Search started' , ISTAT ) DO 20 JS = 1 , NVSEG C C *** Find vertical limits of searched area. C JSTART = JSRCHD + 1 IF ( JS .LT. NVSEG ) THEN JEND = IARR(13) + ( JS*NYS ) / NVSEG ELSE JEND = IARR(15) ENDIF cC cC *** Find vertical limits of used area. cC c IEXT = ( NLPB - JEND + JSTART ) / 2 c IEXT = MAX( NHED , IEXT ) c JUS = MAX( JSTART-IEXT , IXYU(2) ) c JUE = MIN( JEND+IEXT , IXYU(4) , JUS+NLPB-1 ) C C *** Take care of image buffers. C c IF ( JBE .LT. JUE ) THEN c IUSD(2) = JUS c IUSD(4) = JUE IF ( JS .EQ. 1 ) THEN CALL FILBUF( IMF , A , JAPY , NX , IXYU , IUSD , IBUF ) ENDIF c JBS = IBUF(2) c JBE = IBUF(4) c ENDIF C C *** Update sky background net. C CALL SBGNET ( A , JAPY , IBUF , IXYU , JS , & JSTART , JEND , INET , SKYNET , NHSEG , & CRMD , NHED ) C C *** Search JS-th vertical segment. C CALL SEAREG ( A , JAPY , IBUF , IXYU , JSTART , & JEND , NHSEG , INET , SKYNET , ACAT , & BCAT , MCAT , IDET , IDA , IARR , & RARR , M , MM ) C C *** Mark last searched line. C JSRCHD = JEND JMSG = ( 100 * (JSRCHD-IARR(13)+1) ) / NYS C C *** Write message. C WRITE(TEXT,'(I4,A,I8,A)') & JMSG,'% of frame searched ',MM,' objects detected' CALL STTPUT( TEXT , ISTAT ) 20 CONTINUE cC cC ****** Write down objects left in ACAT. cC c MMB = MOD( MM-1 , NC ) + 1 c IF ( MMB .GT. 0 ) THEN c CALL CTLG( ITF , ACAT , NC , START , STEP , MM , MMB ) c END IF C RETURN C END C C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ C C.PURPOSE: Searches horizontal image segment C C.ALGORITHM: Described in additional documentation C C.REMARKS: C----------------------------------------------------------------------- SUBROUTINE SEAREG(A, JAPY, IBUF, IXYU, JSTART, & JEND, NHSEG, INET, SKYNET, ACAT, & BCAT, MCAT, IDET, IDA, IARR, & RARR, M, MM) C IMPLICIT NONE INCLUDE 'MID_REL_INCL:INVENT.INC/NOLIST' C REAL A(1) ! IN: Image buffer INTEGER JAPY(1) ! IN: Pointers to image lines INTEGER IBUF(4) ! IN: Limits of image buffer INTEGER IXYU(4) ! IN: Limits of investigated area INTEGER JSTART ! IN: First searched line INTEGER JEND ! IN: Last searched line INTEGER NHSEG ! IN: Number of horizontal segments INTEGER INET(0:NHSEG) ! IN: Horizontal segments limits REAL SKYNET(2,0:NHSEG) ! IN: Net of sky values REAL ACAT(5,MAXCNT) ! MOD: Array with objects data REAL BCAT(2,MAXDET) ! MOD: Real detection data INTEGER MCAT(4,MAXDET) ! MOD: Integer detection data INTEGER IDET(1) ! MOD: Rolling array with detections INTEGER IDA(1) ! MOD: Definition of nearby pixels INTEGER IARR(32) ! IN: Integer INVENTORY keywords REAL RARR(64) ! IN: Real INVENTORY keywords INTEGER M ! MOD: Number of actual detections INTEGER MM ! MOD: Number of actual objects C INTEGER I , IARG , IB INTEGER IEND , IFLD , IFULL , IHED INTEGER IL , ILIM INTEGER IOF , IOFF INTEGER IP , IPLIM(MAXNET) INTEGER IS , ISDET , ISTART, ISTAT INTEGER IUNIT INTEGER J , JBS , JBE , JXY(4) INTEGER JOF , JOFF C INTEGER KREC, KCOR INTEGER MINCR , MK INTEGER NC , NL , NXS C REAL AVER C REAL AR REAL BGRD , BLIM , BLIME , BLIMS REAL CRMD(3) REAL FLTR , FRCTI1 , FRCTI2 , FRCTJ1 , FRCTJ2 REAL HCUT , LCUT REAL PLIM , PLIML , PRLIM(MAXNET) REAL RTRSH REAL SB2 , SKYMIN , SKYMAX REAL TRSH C REAL TRLM C LOGICAL DETECT , OBJECT C C NC = MAXCNT NL = MAXDET JBS = IBUF(2) JBE = IBUF(4) C C ****** Recall keywords. Sky variations are additive, C ****** limiting threshold RTRSH=RARR(3) is in C ****** measurement units: IFLD=1. Sky variations are C ****** multiplicative, limiting threshold RTRSH=RARR(3) C ****** is expressed in units of local sky: IFLD=0. C IHED = IARR(8) IFLD = IARR(9) IUNIT = IARR(10) NXS = IARR(14) - IARR(12) + 1 MINCR = IARR(22) LCUT = RARR(1) HCUT = RARR(2) RTRSH = RARR(3) FLTR = RARR(4) CRMD(1) = RARR(39) CRMD(2) = RARR(40) CRMD(3) = RARR(41) ILIM = INT( RARR(43) ) C C *** Local sky determination is less accurate when ISDET=0, C *** more accurate when ISDET=1, and most accurate when ISDET=2. C ISDET = IARR(16) C C *** Calculate offsets. C JOFF = JBS - 1 IOFF = IBUF(1) - 1 IOF = IARR(12) - 1 C C *** Find preliminary detection limits. C DO 5 IS = 1 , NHSEG SKYMIN = MIN( SKYNET(1,IS-1) , SKYNET(2,IS-1) , & SKYNET(1,IS) , SKYNET(2,IS) ) SKYMAX = MAX( SKYNET(1,IS-1) , SKYNET(2,IS-1) , & SKYNET(1,IS) , SKYNET(2,IS) ) IF ( IFLD .EQ. 1 ) THEN PRLIM(IS) = SKYMIN + RTRSH TRSH = RTRSH ELSE PRLIM(IS) = SKYMIN * (1.0+RTRSH) TRSH = PRLIM(IS) - SKYMIN ENDIF C C *** Mark segments with disprepant sky determinations. C IF ( SKYMAX-SKYMIN .GT. TRSH ) THEN IPLIM(IS) = 1 ELSE IPLIM(IS) = 0 ENDIF 5 CONTINUE C C *** Search by lines. C DO 10 J = JSTART , JEND JOF = JAPY(J-JOFF) FRCTJ1 = FLOAT(JEND-J) / FLOAT(JEND-JSTART) FRCTJ2 = 1.0 - FRCTJ1 C C *** Search by horizontal segments. C DO 20 IS = 1 , NHSEG IP = IPLIM(IS) IF ( IS .GT. 1 ) THEN ISTART = INET(IS-1) + 1 ELSE ISTART = INET(IS-1) ENDIF IEND = INET(IS) C C *** Find preliminary detection limit. C PLIM = PRLIM(IS) PLIML = MAX( LCUT , PLIM - 4.0*TRSH ) IF ( ISDET .EQ. 2 ) THEN PLIM = PLIM - 0.3 * TRSH ENDIF IB = 0 C C *** Search by pixels. C DO 30 I = ISTART , IEND DETECT = .FALSE. OBJECT = .FALSE. IARG = JOF + I AVER = A(IARG) C C *** Pixel lower than preliminary C *** limit is no longer considered. C IF ( AVER .LT. PLIM ) THEN c IF ( AVER .LT. PLIML ) THEN c MA(IARG) = 0 c ENDIF GOTO 30 ENDIF C C *** Pixels higher than HCUT get special treatment. C IF ( AVER .GE. HCUT ) THEN c MA(IARG) = 2 CALL SATOBJ( A , JAPY , JOFF , I , J , & HCUT , AVER ) IF ( ISDET .EQ. 2 ) THEN BGRD = PLIM - 0.7 * TRSH ELSE BGRD = PLIM - TRSH ENDIF CALL FLTRBP( A , JAPY , IBUF , I , J , & BGRD , FLTR , AVER ) IF ( AVER .GT. 0.9*HCUT ) THEN GOTO 55 ELSE GOTO 56 ENDIF ENDIF C C*********************************************************************** C C ****** This part of code depends on C ****** the detection criterium used. C CALL SRHOBJ( A , JAPY , JOFF , I , J , & DETECT , AVER ) IF ( .NOT. DETECT ) GOTO 30 C C ****** Pixel is higher than 8 neighbours. C C*********************************************************************** C 56 CONTINUE IF ( IB .EQ. 0 ) THEN C C ****** Find background on borders. C BLIMS = SKYNET(1,IS-1)*FRCTJ1 + & SKYNET(2,IS-1)*FRCTJ2 + TRSH BLIME = SKYNET(1,IS)*FRCTJ1 + & SKYNET(2,IS)*FRCTJ2 + TRSH IB = 1 ENDIF C C ****** Calculate interpolated background. C FRCTI1 = FLOAT(I-ISTART) / FLOAT(IEND-ISTART) FRCTI2 = 1.0 - FRCTI1 BLIM = BLIMS*FRCTI1 + BLIME*FRCTI2 IF ( ISDET .EQ. 2 ) THEN BLIM = BLIM - 0.3 * TRSH ENDIF IF ( AVER .GT. BLIM ) THEN IF ( ISDET .EQ. 2 ) THEN BGRD = BLIM - 0.7 * TRSH ELSE BGRD = BLIM - TRSH ENDIF GOTO 40 ENDIF C C ****** Calculate approximate background. C JXY(1) = MAX( IBUF(1) , I-IHED ) JXY(2) = MAX( IBUF(2) , J-IHED ) JXY(3) = MIN( IBUF(3) , I+IHED ) JXY(4) = MIN( IBUF(4) , J+IHED ) CALL APRBGR ( A , JAPY , JOFF , JXY , SB2 ) IF ( AVER .LT. SB2+TRSH ) THEN GOTO 30 ELSE BGRD = SB2 ENDIF 40 CONTINUE C C ****** The pixel has passed approximate C ****** criteria. Now the final criterium C ****** will be applied if necessary. C IF ( IP .EQ. 1 .OR. ISDET .EQ. 1 ) THEN IFULL = 0 CALL SKYMOD( A , JAPY , IBUF , I , J , CRMD , & IHED , IFULL , BGRD ) IF ( AVER .LE. BGRD+TRSH ) GOTO 30 IF ( A(JAPY(J-JOFF)+I) .GT. & BGRD + (FLTR-1.0)*TRSH ) THEN CALL FLTRBP( A , JAPY , IBUF , I , J , & BGRD , FLTR , AVER ) ENDIF IF ( AVER .LE. BGRD+TRSH ) GOTO 30 ENDIF C 50 CONTINUE C C *** One more check. C c TRLM = BGRD + TRSH c CALL RADDET( A , JAPY , IBUF , I , J , c & MINCR , TRLM , AVER , TRSH , AR ) c IF ( AR .LT. 0.0 .AND. AVER .LT. 0.9*HCUT ) GOTO 30 55 CONTINUE C C ****** This is a detection. C M = M + 1 IF (M.GE.NL) THEN CALL STTPUT + ('*** FATAL: Internal buffer overflow; ',ISTAT) CALL STTPUT + (' Restrict our search to smaller subframe', + ISTAT) CALL STTPUT + (' or modify parameter setup fro detection', + ISTAT) CALL STSEPI ENDIF C C ****** Save arrays MCAT and BCAT C ****** if they are full of new data. C c IF ( MOD( M , NL ) .EQ. 1 .AND. M .GT. NL ) THEN c KCOR = M - NL - 1 c DO 60 K = 1 , NL c KREC = K + KCOR c WRITE ( ISF , REC=KREC ) MCAT(1,K) , c & MCAT(2,K) , MCAT(3,K) , MCAT(4,K) , c & BCAT(1,K) , BCAT(2,K) c 60 CONTINUE c ENDIF C C *** Record data for new detection. C MK = MOD( M-1 , NL ) + 1 MCAT(1,MK) = I MCAT(2,MK) = J MCAT(3,MK) = 0 MCAT(4,MK) = 0 BCAT(1,MK) = BGRD BCAT(2,MK) = AVER IP = I - IOF C C *** Update linked list of detections. C CALL UPDTLL( MCAT , NL , IDET , & NXS , ILIM , IDA , IP , M ) 30 CONTINUE 20 CONTINUE C C *** Join multiple detections and update array IDET. C CALL JOINMD( A , JAPY , IBUF , IXYU , ACAT , & NC , BCAT , MCAT , NL , IDET , & NXS , ILIM , IARR , RARR , M , & MM ) IF ( J .EQ. IARR(15) ) THEN DO 70 IL = 1 , ILIM CALL JOINMD( A , JAPY , IBUF , IXYU , ACAT , & NC , BCAT , MCAT , NL , IDET , & NXS , ILIM , IARR , RARR , M , & MM ) 70 CONTINUE ENDIF 10 CONTINUE C RETURN C END