C @(#)search.for	19.1 (ES0-DMD) 02/25/03 13:25:39
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
