C @(#)histo.for	19.1 (ESO-DMD) 02/25/03 13:30: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 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 HISTO(ARRAY,NAXIS,NPIX,PIXELS,F,CUTS,NSLOT,SLOT)
C
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C.COPYRIGHT: Copyright (c) 1988 European Southern Observatory,
C                                         all rights reserved
C
C.LANGUAGE: F77+ESOext
C
C.AUTHOR: K.Banse
C
C.IDENTIFICATION
C  subroutine HISTO		version 2.02	830615
C             HSTVLS			2.00841130
C             HISTEQ			1.90	830331
C  K. Banse                        ESO - Garching
C
C.KEYWORDS
C  histogram
C
C.PURPOSE
C  1) calculate histogram of given array
C  2) equalize histogram
C
C.ALGORITHM
C  straigth forward + from Gonzalez,Wintz: "Digital Image Processing"
C
C.INPUT/OUTPUT
C  call as HISTO(ARRAY,NAXIS,NPIX,PIXELS,F,CUTS,NSLOT,SLOT)
C
C 1)
C  input par:
C  ARRAY:      R*4 array      input data
C  NAXIS:      I*4            no. of axes of ARRAY
C  NPIX:      I*4 array      original size of each axis
C  PIXELS:      I*4 array      start + end pixels in each axis
C  F:            R*4            scaling factor
C  CUTS:      R*4 array      low + high cuts of ARRAY
C  NSLOT:      I*4            no. of bins (without excess bins...)
C
C  output par:
C  SLOT:      I*4 array      bins for histogram
C
C 2)
C  input par:
C  BINS:      I*4 array      histogram
C  LEVEL:      I*4            no. of levels in histogram (= dim. of BINS)
C  FSIZE:      R*4            no. of pixels in image
C  KBINS:      I*4 array      equalized histogram
C  ITT:            I*2 array      ITT table to realize the equalized histogram on the DeAnza
C
C------------------------------------------------------------------------
C 
      IMPLICIT NONE
C
      INTEGER      NAXIS,NPIX(3),PIXELS(2,2),NSLOT,SLOT(1)
      INTEGER      LOWX,LOWY,HIX,HIY,NX,NY,N,IOFF,X,NOPIX
C 
      REAL         ARRAY(1)
      REAL         F,CUTS(2)
C
C  clear slots
      DO 100 N=1,NSLOT+2
         SLOT(N) = 0
100   CONTINUE
C
C  determine subarea
      LOWX = PIXELS(1,1)
      HIX = PIXELS(1,2)
      IF (NAXIS.GE.2) THEN
         LOWY = PIXELS(2,1)
         HIY = PIXELS(2,2)
      ELSE
         LOWY = 1
         HIY = 1
      ENDIF
C
C  main loop over all pixels in given area
      NOPIX = NPIX(3)
      IOFF = (LOWY - 1) * NOPIX
C 
      DO 1000 NY=LOWY,HIY
         DO 600 NX=LOWX,HIX
            N = IOFF + NX
            IF (ARRAY(N).GT.CUTS(2)) THEN
               X = NSLOT                        !high excess bin
               GOTO 500
            ENDIF
            IF (ARRAY(N).LT.CUTS(1)) THEN
               X = 1                        	!low excess bin
               GOTO 500
            ENDIF
C
            X = INT(F*(ARRAY(N) - CUTS(1))) + 2
500         SLOT(X) = SLOT(X) + 1
            IOFF = IOFF + NOPIX
600      CONTINUE
1000  CONTINUE
C
C  that's it folks...
      RETURN
      END

      SUBROUTINE
     +HSTVLS(ARRAY,NAXIS,NPIX,SUBLO,SUBHI,CUTS,SLTSIZ,NSLOT,SLOT)
C
      IMPLICIT NONE
C 
      INTEGER      NAXIS,NPIX(3),SUBLO(3),SUBHI(3),NSLOT,SLOT(1)
      INTEGER      LOWX,LOWY,LOWZ,HIX,HIY,HIZ
      INTEGER      N,OFF,YOFFA,YOFF,ZOFF,NX,NY,NZ
      INTEGER      NPX,NPXY,X
C 
      REAL         ARRAY(1)
      REAL         SLTSIZ,CUTS(2),F,R
C
C  clear slots
      DO 100 N=1,NSLOT
         SLOT(N) = 0
100   CONTINUE
      F = 1./SLTSIZ
C
C  determine subarea
      LOWX = SUBLO(1)
      HIX = SUBHI(1)
      NPX = NPIX(1)
      IF (NAXIS.GE.2) THEN
         LOWY = SUBLO(2)
         HIY = SUBHI(2)
         NPXY = NPX * NPIX(2)
      ELSE
         LOWY = 1
         HIY = 1
         NPXY = NPX
      ENDIF
      IF (NAXIS.GE.3) THEN
         LOWZ = SUBLO(3)
         HIZ = SUBHI(3)
      ELSE
         LOWZ = 1
         HIZ = 1
      ENDIF
C
      ZOFF = (LOWZ-1) * NPXY
      YOFFA = (LOWY-1) * NPX
C
C  test, if we have excess bins
      IF (CUTS(2).LE.CUTS(1)) GOTO 1000
C
C  main loop over all pixels in given area with excess bins
      DO 800 NZ=LOWZ,HIZ
         YOFF = YOFFA                           !reset Y offset for each z-loop
C
         DO 600 NY=LOWY,HIY
            OFF = ZOFF + YOFF
C
            DO 500 NX=LOWX,HIX
               N = OFF + NX
               IF (ARRAY(N).GT.CUTS(2)) THEN
                  X = NSLOT                     !high excess bin
               ELSE
                  R = ARRAY(N) - CUTS(1)
                  IF (R.LT.0.) THEN
                     X = 1                  	!low excess bin
                  ELSE
                     X = INT(F*R) + 2           !valid bin
                  ENDIF
               ENDIF
               SLOT(X) = SLOT(X) + 1
500         CONTINUE
C
            YOFF = YOFF + NPX
600      CONTINUE
C
         ZOFF = ZOFF + NPXY
800   CONTINUE
C
C  that's it
      RETURN
C
C  main loop over all pixels in given area without excess bins
1000  DO 1800 NZ=LOWZ,HIZ
         YOFF = YOFFA                           !reset Y offset for each z-loop
C
         DO 1600 NY=LOWY,HIY
            OFF = ZOFF + YOFF
C
            DO 1500 NX=LOWX,HIX
               N = OFF + NX
               X = INT(F*(ARRAY(N)-CUTS(1))) + 1
               SLOT(X) = SLOT(X) + 1
1500        CONTINUE
C
            YOFF = YOFF + NPX
1600     CONTINUE
C
         ZOFF = ZOFF + NPXY
1800  CONTINUE
C
C  that's it folks...
      RETURN
      END

      SUBROUTINE HISTEQ(BINS,LEVEL,FSIZE,KBINS,ITT)
C
      IMPLICIT NONE
C 
      INTEGER     BINS(1),LEVEL,KBINS(1),ITT(1)
      INTEGER     IS(512)
      INTEGER     N,NN,NST
C 
      REAL        T(512),S(512),RLEVL,FSIZE
C
      RLEVL = FLOAT((LEVEL-1) * 2)
      T(1) = 1./RLEVL                        !T(n) will be set to midpoints between intervals
      S(1) = FLOAT(BINS(1))/FSIZE            !S(n) is the transfer function
C
      DO 100 N=2,LEVEL
         T(N) = (2*N-1)/RLEVL
         S(N) = S(N-1) + FLOAT(BINS(N))/FSIZE
100   CONTINUE
C
C  assign values of transfer function to nearest slot
      NST = 1                        	!start at slot 1
      DO 1000 N=1,LEVEL
C
         IF (NST.EQ.LEVEL) THEN            !if last slot already reached,
            IS(N) = LEVEL                  !index must be the last possible one...
         ELSE
C
            DO 500 NN=NST,LEVEL             !use T(n), to find nearest slot
               IF (S(N).LT.T(NN)) THEN      !in increasing order
                  IS(N) = NN
                  NST = NN
                  GOTO 1000
               ENDIF
500         CONTINUE
         ENDIF
1000  CONTINUE
C
C  calculate equalized histogram
      DO 2100 N=1,LEVEL
         KBINS(N) = 0
2100  CONTINUE
C
      DO 2500 N=1,LEVEL
         NN = IS(N)
         KBINS(NN) = KBINS(NN) + BINS(N)
         ITT(N) = NN - 1
2500  CONTINUE
C
      RETURN
      END
