C @(#)aini.for	17.1.1.1 (ES0-DMD) 01/25/02 17:11:12
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 aini(NX,NY,MT,IN,H0,HG,HL,EG,EL,SG,SL,YMT)
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C      initializes the algorithm:
C      EG,SG,EL,SL: exspectation value and standard deviation
C      of gradients and Laplace-terms resp. by noise (out)
C      other parameters as in afido
C
C------------------------------------------------------------------------
       IMPLICIT NONE
C
       LOGICALL
       CHARACTER*1 YMT
       INTEGER     NX,NY,I1,J1,N1X,N1Y,I,J,T
       REAL        IN(NX,NY),H0(NX,NY),HG(NX,NY),HL(NX,NY),
     >             EG,EL,SG,SL,A,C,B,D,HX,HY,H,QG,QL,NP,NPP,
     >             MT(NX,NY)
C
C *** 1.order:smoothing
      DO  J=1,NY-1
        J1=J+1
        DO  I=1,NX-1
          I1=I+1
          H0(I1,J1)=(IN(I,J)+IN(I1,J)+IN(I,J1)+IN(I1,J1))/4.
        ENDDO
      ENDDO
C
C *** set edge lines: up and left and upper left corner pixel
      H0(1,1)=H0(2,2)                        ! corner
C
      DO I=2,NX
         H0(I,1)=H0(I,2)                     ! 1.row
      ENDDO
C
      DO J=2,NY
         H0(1,J)=H0(2,J)                     ! 1.column
      ENDDO
C
      L=.FALSE.
      NP=0.
      EL=NP
      EG=NP
      QL=NP
      QG=NP
C
C *** 2.order: centering
      DO J=1,NY-1
         J1=J+1
         DO I=1,NX-1
            I1=I+1
            A=H0(I,J)
            B=H0(I1,J)
            C=H0(I,J1)
            D=H0(I1,J1)
            H=(A+B+C+D)/4.
            H0(I,J)=H                        ! smooth
            H=IN(I,J)-H
            HL(I,J)=H                        ! Laplace
            IF (YMT.EQ.'N') GOTO 15          ! no mask for statistics
            T=MT(I,J)
            IF (T.NE.0) L=.TRUE.
            IF (L) GOTO 16                   ! pixel not used for statistics
   15       NPP=NP
            NP=NP+1
            NPP=NPP/NP
            H=H-EL
            EL=EL+H/NP                       ! estimate statistics
            QL=QL+H*H*NPP                    ! for Laplace
   16       HX=(A-B+C-D)/4.
            HY=(A+B-C-D)/4.
            H=SQRT(HX*HX+HY*HY)
            HG(I,J)=H                        ! gradient
            IF (YMT.EQ.'N') GOTO 17          ! no mask for statistics
            IF (L) GOTO 20                   ! pixel not used for statistics
   17       H=H-EG
            EG=EG+H/NP                       ! estimate statistics
            QG=QG+H*H*NPP                    ! for gradient
   20       L=.FALSE.
         ENDDO
      ENDDO

C correct statistics
      NP=NP-1
      SL=SQRT(QL/NP)
      SG=SQRT(QG/NP)
C
C *** set edge lines down and right and lower right corner pixel
      N1X=NX-1
      N1Y=NY-1
      H0(NX,NY)=H0(N1X,N1Y)                  ! corner
      HG(NX,NY)=HG(N1X,N1Y)
      HL(NX,NY)=HL(N1X,N1Y)

      DO I=1,N1X                             ! last row
         H0(I,NY)=H0(I,N1Y)
         HG(I,NY)=HG(I,N1Y)
         HL(I,NY)=HL(I,N1Y)
      ENDDO

      DO J=1,N1Y                             ! last column
         H0(NX,J)=H0(N1X,J)
         HG(NX,J)=HG(N1X,J)
         HL(NX,J)=HL(N1X,J)
      ENDDO
C
      RETURN
      END