C++********************************************************************* C C HIST.F FMINT FOR HIST RANGE JUN 2000 ArDean Leith C NaN OK JAN 2001 ArDean Leith C OLD MODE BUG JAN 2002 ArDean Leith C FNUMEL OVERFLOW BUG JUL 2002 ArDean Leith C INTEGER OUTPUT OCT 2002 ArDean Leith C C ********************************************************************** C=* FROM: SPIDER - MODULAR IMAGE PROCESSING SYSTEM. AUTHOR: J.FRANK * C=* Copyright (C) 1985-2005 Health Research Inc. * C=* * C=* HEALTH RESEARCH INCORPORATED (HRI), * C=* ONE UNIVERSITY PLACE, RENSSELAER, NY 12144-3455. * C=* * C=* Email: spider@wadsworth.org * 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 the * C=* 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 GNU * C=* General Public License for more details. * C=* * C=* You should have received a copy of the GNU General Public License * C=* along with this program; if not, write to the * C=* Free Software Foundation, Inc., * C=* 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. * C=* * C ********************************************************************** C C HIST(LUN,LUNMA,LUNDOC,NSAM,NROW,NSLICE,IDUST,HMIN,HMAX) C C PURPOSE: COMPUTE 128 PLACE HISTOGRAM FROM IMAGE RECORDS, C DISPLAY HISTOGRAM ON LINE PRINTER, TERMINAL, OR C IN DOC. FILE AND, OPTIONALLY, REMOVE DATA THAT ARE C OUT OF A SPECIFIED STATISTICAL RANGE. C C PARAMETERS: LUN IO UNIT NUMBER OF IMAGE FILE C LUNDMA IO UNIT NUMBER FOR MASK FILE C LUNDOC IO UNIT NUMBER FOR DOCUMENT FILE C NSAM,NROW DIMENSIONS OF IMAGE C NSLICE DIMENSIONS OF IMAGE C HMIN,HMAX HIST. MIN & MAX (RET.) C HSIG,HMODE HIST. S.D. & MODE (USED BY DUST) C C23456789 123456789 123456789 123456789 123456789 123456789 123456789 12 C--********************************************************************* SUBROUTINE HIST(LUN,LUNMA,LUNDOC,NSAM,NROW,NSLICE, & HMIN,HMAX,HSIG,HMODE) INCLUDE 'CMBLOCK.INC' DIMENSION FREQ(128) DIMENSION PLIST(4) DOUBLE PRECISION HAV, HAV2, DTOP, FNUMEL CHARACTER *1 ANS,NULL LOGICAL MASK NULL = CHAR(0) IF (IFORM .LE. 0) THEN C CAN NOT DO HISTOGRAM ON FOURIER FILES CALL ERRT(2,'HISTM',NE) RETURN ENDIF C MAKE SURE STATISTICS ARE CURRENT IF (IMAMI .NE. 1) CALL NORM3(LUN,NSAM,NROW,NSLICE,FMAX,FMIN,AV) HMAX = FMAX HMIN = FMIN C NBINS MUST BE <= 1000!!! NBINS = 128 NDEV = NDAT IDEV = 1 IF (FCHAR(4:4) .EQ. 'R' .OR. FCHAR(5:5) .EQ. 'R') THEN C ONE OF THE LETTERS IS 'R' ASK USER FOR HISTOGRAM RANGE CALL RDPRM2(HMIN,HMAX,NOT_USED,'HISTOGRAM RANGE MIN, MAX') ELSEIF (FCHAR(4:4) .EQ. 'M') THEN C ASK USER FOR HISTOGRAM SINK CALL RDPRMC(ANS,NLET,.TRUE., & 'OUTPUT TO RESULTS FILE, DOC. FILE, OR TERMINAL? (R/D/T)', & NULL,IRTFLG) ENDIF IF (FCHAR(4:4) .EQ. 'T' .OR. FCHAR(5:5) .EQ. 'T' .OR. & (FCHAR(4:4) .EQ. 'M' .AND. ANS .EQ. 'T')) THEN C OUTPUT TO TERMINAL, NOT RESULTS FILE NBINS = 70 NDEV = NOUT IDEV = 0 ENDIF C ZERO THE HISTOGRAM FREQUENCIES DO K = 1,NBINS FREQ(K) = 0.0 ENDDO C MPI 3-MAR-80 ADDED TO AVOID "BEAT" PHENOMENON (REMOVED AUG 98) C APPEARED TO EXPAND RANGE OF HIST. MIN/MAX FOR < 0...2 MASK = .FALSE. NREC = NSLICE * NROW IF (FCHAR(4:4) .EQ. 'M') THEN C HISTOGRAM UNDER MASK ONLY MASK = .TRUE. C DETERMINE RANGE UNDER MASK CALL HISTMINMAX(LUN,LUNMA,NSAM,NREC,HMIN,HMAX) ENDIF HDIFF = HMAX - HMIN FF = (NBINS - 1.0) / HDIFF BSIZ = 1.0 / FF FNUMEL = 0.0 HAV = 0.0 HAV2 = 0.0 CALL HISTIMAGE(LUN,LUNMA,NSAM,NREC,FREQ, & FNUMEL,HAV,HAV2,HMIN,FF,NBINS,MASK) IF (FCHAR(4:4) .EQ. 'D' .OR. FCHAR(5:5) .EQ. 'D' .OR. & (FCHAR(4:4) .EQ. 'M' .AND. ANS .EQ. 'D')) THEN C OUTPUT TO DOC FILE PLIST(2) = HMIN DO IBIN=1,NBINS PLIST(1) = IBIN PLIST(2) = HMIN + (BSIZ *(IBIN - 1)) PLIST(3) = FREQ(IBIN) CALL SAVD(LUNDOC,PLIST,3,IRTFLG) ENDDO CALL SAVDC CLOSE(LUNDOC) ELSE C OUTPUT TO RESULTS FILE OR TERMINAL WRITE(NDEV,*) ' ' C GRAPHX IN GRAPHS USES FMIN & FMAX FOR HISTOGRAM LABELS FMINT = FMIN FMIN = HMIN FMAXT = FMAX FMAX = HMAX CALL GRAPHS(NDEV,FREQ,NBINS,1,IDEV,1.0,IRTFLG) FMIN = FMINT FMAX = FMAXT WRITE(NDEV,*) ' ' IF (IRTFLG .NE. 0) THEN CALL ERRT(100,'HIST',NE) RETURN ENDIF ENDIF C FIND MAXIMUM FREQUENCY OCCURING IN HISTOGRAM (HISMAX) & LOCATION HISMAX = FREQ(1) MAXBIN = 1 DO IBIN = 2,NBINS IF (FREQ(IBIN) .GE. HISMAX) THEN HISMAX = FREQ(IBIN) MAXBIN = IBIN ENDIF ENDDO C CONVERT LOCATION MAXBIN TO CORRESPONDING IMAGE INTENSITY (MODE) IF (MAXBIN .EQ. 1) THEN C SUB-BIN ESTIMATE OF MODE BMODE = 0.5 ELSEIF (MAXBIN .EQ. NBINS) THEN C SUB-BIN ESTIMATE OF MODE C (MAYBE THIS SHOULD AVERAGE ALL BINS AT END WITH SAME NUMBER?) BMODE = FLOAT(NBINS) - 0.5 ELSE C SUB-BIN ESTIMATE OF MODE YM1 = FREQ(MAXBIN-1) YP1 = FREQ(MAXBIN+1) BMODE = FLOAT(MAXBIN-1) + (YM1-YP1)*0.5/ (YM1+YP1-2.*HISMAX) ENDIF HMODE = HMIN + BMODE * BSIZ DTOP = HAV2 - HAV * HAV / FNUMEL IF (DTOP .LT. 0.0) THEN C SQRT OF NEGATIVE NUMBER WRITE(NOUT,*) '*** in HIST, SQRT(',DTOP,') IMPOSSIBLE' CALL ERRT(100,'HIST',NE) RETURN ELSEIF (FNUMEL .EQ. 1) THEN WRITE(NOUT,*) '*** CAN NOT DETERMINE S.D. (ONLY ONE PIXEL) ' CALL ERRT(100,'HIST',NE) RETURN ENDIF HAV = HAV / FNUMEL HSIG = DSQRT(DTOP / (FNUMEL - 1)) FNPIX = NSAM * NROW * NSLICE FNBINS = NBINS WRITE(NDEV,*) ' ' IBIG = HUGE(IBIG) IF (FNPIX .LT. IBIG) THEN NPIX = NSAM * NROW * NSLICE NUMEL = FNUMEL WRITE(NDEV,90) FMIN,FMAX,HMIN,HMAX,NPIX,NUMEL, & NBINS,BSIZ,HAV,HMODE,HSIG 90 FORMAT( & ' FILE RANGE: ',1PG11.4,' ......... ',1PG11.4,/, & ' HISTOGRAM RANGE: ',1PG11.4,' ......... ',1PG11.4,/, & ' TOTAL PIXELS: ',I11, ' PIXELS IN HIST.: ',I11,/, & ' NO. OF BINS: ',I11, ' BIN SIZE: ',1PG11.4,/, & ' HIST. MEAN: ',1PG11.4,' HIST. MODE: ',1PG11.4,/, & ' HIST. S.D.: ',1PG11.4/) ELSE WRITE(NDEV,91) FMIN,FMAX,HMIN,HMAX,FNPIX,FNUMEL, & FNBINS,BSIZ,HAV,HMODE,HSIG 91 FORMAT( & ' FILE RANGE: ',1PG11.4,' ......... ',1PG11.4,/, & ' HISTOGRAM RANGE: ',1PG11.4,' ......... ',1PG11.4,/, & ' TOTAL PIXELS: ',1PG11.4,' PIXELS IN HIST.: ',1PG11.4,/, & ' NO. OF BINS: ',1PG11.4,' BIN SIZE: ',1PG11.4,/, & ' HIST. MEAN: ',1PG11.4,' HIST. MODE: ',1PG11.4,/, & ' HIST. S.D.: ',1PG11.4/) ENDIF WRITE(NDEV,*) ' ' IF (NDEV .NE. NOUT) THEN WRITE(NOUT,*) ' ' IF (FNPIX .LT. IBIG) THEN WRITE(NOUT,90) FMIN,FMAX,HMIN,HMAX,NPIX, & NUMEL,NBINS,BSIZ,HAV,HMODE,HSIG ELSE WRITE(NOUT,91) FMIN,FMAX,HMIN,HMAX,NPIX, & FNUMEL,FNBINS,BSIZ,HAV,HMODE,HSIG ENDIF WRITE(NOUT,*) ' ' ENDIF RETURN END C -------------------- HISTIMAGE -------------------------------- SUBROUTINE HISTIMAGE(LUN,LUNMA,NSAM,NREC,FREQ, & FNUMEL,HAV,HAV2,HMIN,FF,NBINS,MASK) INCLUDE 'CMBLOCK.INC' DIMENSION FREQ(128), BUFMASK(NSAM), REDBUF(NSAM) DOUBLE PRECISION HAV, HAV2, FNUMEL LOGICAL MASK C GET HISTOGRAM FROM IMAGE VALUES IF (MASK) THEN C MASKED, HANDLES NaN for NOT MASK OK DO IREC=1,NREC CALL REDLIN(LUN,REDBUF,NSAM,IREC) CALL REDLIN(LUNMA,BUFMASK,NSAM,IREC) DO ISAM = 1,NSAM IF (BUFMASK(ISAM) .GE. 0.5) THEN C HISTOGRAM THIS PIXEL, MASK HAS POSITIVE VALUE) C FIND BIN NUMBER BVAL = REDBUF(ISAM) JBIN = INT((BVAL - HMIN) * FF) + 1.5 IF (JBIN.GE.1 .AND. JBIN.LE.NBINS) THEN C WITHIN HISTOGRAM RANGE FREQ(JBIN) = FREQ(JBIN) + 1.0 FNUMEL = FNUMEL + 1 HAV = HAV + BVAL HAV2 = HAV2 + DBLE(BVAL) * DBLE(BVAL) ENDIF ENDIF ENDDO ENDDO ELSE C NO MASK DO IREC=1,NREC CALL REDLIN(LUN,REDBUF,NSAM,IREC) DO ISAM = 1,NSAM C HISTOGRAM THIS PIXEL C FIND BIN NUMBER BVAL = REDBUF(ISAM) JBIN = INT((BVAL - HMIN) * FF) + 1.5 IF (JBIN.GE.1 .AND. JBIN.LE.NBINS) THEN C WITHIN HISTOGRAM RANGE FREQ(JBIN) = FREQ(JBIN) + 1.0 FNUMEL = FNUMEL + 1 HAV = HAV + BVAL HAV2 = HAV2 + DBLE(BVAL) * DBLE(BVAL) ENDIF ENDDO ENDDO ENDIF RETURN END C -------------------- HISTMINMAX -------------------------------- SUBROUTINE HISTMINMAX(LUN,LUNMA,NSAM,NREC, & HMIN,HMAX) INCLUDE 'CMBLOCK.INC' DIMENSION BUFMASK(NSAM), REDBUF(NSAM) C DETERMINE RANGE UNDER MASK DO IREC=1, NREC CALL REDLIN(LUN, REDBUF,NSAM,IREC) CALL REDLIN(LUNMA,BUFMASK,NSAM,IREC) DO ISAM = 1,NSAM IF (BUFMASK(ISAM) .GE. 0.5) THEN C PIXEL HAS POSITIVE MASK VALUE BVAL = REDBUF(ISAM) IF (BVAL .LT. HMIN) HMIN = BVAL IF (BVAL .GT. HMAX) HMAX = BVAL ENDIF ENDDO ENDDO RETURN END