C @(#)fndedg.for	17.1.1.1 (ES0-DMD) 01/25/02 17:40:34
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
      SUBROUTINE FNDEDG(METHOD,A,NX,NY,THR,E)
C
C++++++++++++++++++++++++++++++++++++++++++++++++++
C
C.IDENTIFICATION
C
C  subroutine FNDEDG		version 1.00	880912
C  K. Banse                  	ESO - Garching
C
C.KEYWORDS
C
C  edge detection
C
C.PURPOSE
C
C  find edges in a frame
C
C.ALGORITHM
C
C  method = Band, see Gonzalez + Wintz: "Digital Image Processing", chapter 7
C         = Sobel, see Pratt: "Digital Image Processing", chapter 17.4.2
C
C.INPUT/OUTPUT
C
C  call as FNDEDG(METHOD,A,NX,NY,THR,E)
C
C  input par:
C  METHOD:	char. exp.	method for edge detection
C				= B, compare neighbours "above" + "behind"
C				= S, use Sobel 3x3 edge operator
C  A:		R*4 array	input data frame
C  NX:		I*4		no. of rows
C  NY:		I*4		no. of columns
C  THR:		R*4		threshold for different bands
C
C  output par:
C  E:		R*4 array	output frame
C
C.VERSIONS
C  1.00		from version 2.02 as of 830506
C  		move to FORTRAN 77 + new St interfaces
C 
C--------------------------------------------------
C
      IMPLICIT NONE
C 
      CHARACTER*(*)   METHOD
C 
      REAL        THR,A(*),E(*)
      REAL        X1,Y1,X2,Y2,XGRAD,YGRAD,C1,C2
C 
      INTEGER     NX,NY
      INTEGER     NIN,NOUT
      INTEGER     OFF,ABOVE,BELOW
C      
C  branch on method
      IF (METHOD(1:1).EQ.'S') GOTO 5000
C
C  handle 1. line specially
      E(1) = 0.0
      C1 = A(1)
C 
      DO 400 NIN=2,NX
         C2 = A(NIN)
         IF (C1.LT.THR) THEN
            IF (C2.GE.THR) THEN
               E(NIN) = 1.0
            ELSE
               E(NIN) = 0.0
            ENDIF
         ELSE
            IF (C2.LT.THR) THEN
               E(NIN) = 1.0
            ELSE
               E(NIN) = 0.0
            ENDIF
         ENDIF
         C1 = C2
400   CONTINUE
C
C  for the other lines look for edges "behind" + "above"
C 
      OFF = NX + 1
      DO 1000 NOUT=2,NY
         E(OFF) = 0.0
         C1 = A(OFF)
C
C  1. inner loop along each row
         DO 600 NIN=OFF+1,NOUT*NX
            C2 = A(NIN)
            IF (C1.LT.THR) THEN
               IF (C2.GE.THR) THEN
                  E(NIN) = 1.0
               ELSE
                  E(NIN) = 0.0
               ENDIF
            ELSE
               IF (C2.LT.THR) THEN
                  E(NIN) = 1.0
               ELSE
                  E(NIN) = 0.0
               ENDIF
            ENDIF
            C1 = C2
600      CONTINUE
C  
C  2. inner loop along two adjacent rows
         DO 800 NIN=OFF,NOUT*NX
            C1 = A(NIN-NX)
            C2 = A(NIN)
            IF (C1.LT.THR) THEN
               IF (C2.GE.THR) E(NIN) = 1.0
            ELSE
               IF (C2.LT.THR) E(NIN) = 1.0
            ENDIF
800      CONTINUE
C
         OFF = OFF + NX
1000  CONTINUE
C
C  that's it folks...
      RETURN
C      
C ------------------------------------------------------
C      
C  use Sobel edge detector: 1  2  1  for y-gradient and 1  0 -1 for x-gradient
C                           0  0  0                     2  0 -2
C                          -1 -2 -1                     1  0 -1
C  grad = abs(grady) + abs(gradx)
C      
C ------------------------------------------------------
C      
C    handle 1. line specially...
5000  E(1) = 0.0
      E(NX) = 0.0
C 
      DO 5400 NIN=2,NX-1
         XGRAD = ABS(A(NIN-1)-A(NIN+1))
         IF (XGRAD.GT.THR) THEN
            E(NIN) = 1.0
         ELSE
            E(NIN) = 0.0
         ENDIF
CCC         E(NIN) = XGRAD
5400  CONTINUE
C      
C  handle last line specially...
      OFF = (NY-1)*NX + 1
      E(OFF) = 0.0
      E(NX*NY) = 0.0
C 
      DO 5600 NIN=OFF+1,OFF+NX-1
         XGRAD = ABS(A(NIN-1)-A(NIN+1))
         IF (XGRAD.GT.THR) THEN
            E(NIN) = 1.0
         ELSE
            E(NIN) = 0.0
         ENDIF
CCC         E(NIN) = XGRAD
5600  CONTINUE
C      
C  now work on the "safe" part of the image
      OFF = NX
      DO 7000 NOUT=2,NY-1
         ABOVE = OFF + NX
         BELOW = OFF - NX
         E(OFF+1) = 0.0
         E(OFF+NX) = 0.0
C      
         DO 6000 NIN = 2,NX-1
            X1 = A(BELOW+NIN-1) + 2.*A(OFF+NIN-1) + A(ABOVE+NIN-1)
            X2 = A(BELOW+NIN+1) + 2.*A(OFF+NIN+1) + A(ABOVE+NIN+1)
            XGRAD = ABS(X1-X2)
            Y1 = A(BELOW+NIN-1) + 2.*A(BELOW+NIN) + A(BELOW+NIN+1)
            Y2 = A(ABOVE+NIN-1) + 2.*A(ABOVE+NIN) + A(ABOVE+NIN+1)
            YGRAD = ABS(Y1-Y2)
            IF (MAX(XGRAD,YGRAD).GT.THR) THEN
               E(NIN+OFF) = 1.0
            ELSE
               E(NIN+OFF) = 0.0
            ENDIF
CCC            E(NIN) = XGRAD+YGRAD
6000     CONTINUE
C  
         OFF = OFF + NX
7000  CONTINUE
C
C  that's it folks...
      RETURN
      END