C @(#)histo.for 17.1.1.1 (ESO-DMD) 01/25/02 17:18:59 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