C++********************************************************************* C C HISD.F C ADDED NBINS OCT 2006 ARDEAN LEITH C ********************************************************************** C=* * C=* This file is part of: SPIDER - Modular Image Processing System. * C=* SPIDER System Authors: Joachim Frank & ArDean Leith * C=* Copyright 1985-2010 Health Research Inc., * C=* Riverview Center, 150 Broadway, Suite 560, Menands, NY 12204. * C=* Email: spider@wadsworth.org * C=* * C=* SPIDER 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=* SPIDER 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=* You should have received a copy of the GNU General Public License * C=* along with this program. If not, see * C=* * C ********************************************************************** C C HISD(NDOC) C C PURPOSE: COMPUTE HISTOGRAM FROM ONE COLUMN OF C DOCUMENT FILE, DISPLAY HISTOGRAM ON C LINE PRINTER OR IN DOC FILE. C C NOTES: FCHAR = HD - FULL HISTOGRAM TO RESULTS FILE C NOTES: FCHAR = HD D - FULL HISTOGRAM TO DOC FILE C NOTES: FCHAR = HD RD - RANGE HISTOGRAM TO DOC FILE C NOTES: FCHAR = HD R - RANGE HISTOGRAM TO RESULTS FILE C C23456789 123456789 123456789 123456789 123456789 123456789 123456789 12 C--********************************************************************* SUBROUTINE HISD(NDOC) INCLUDE 'CMBLOCK.INC' INCLUDE 'CMLIMIT.INC' INCLUDE 'F90ALLOC.INC' REAL, DIMENSION(:,:), POINTER :: PBUF CHARACTER(LEN=MAXNAM) :: DOCFIL REAL, ALLOCATABLE, DIMENSION(:) :: FREQ DIMENSION PLIST(4) MAXX = 0 MAXY = 0 CALL GETDOCDAT('INPUT DOCUMENT',.TRUE.,DOCFIL,NDOC, & .TRUE.,MAXX, MAXY,PBUF,IRTFLG) IF (IRTFLG .NE. 0) RETURN CALL RDPRI1S(ICOL,NOT_USED, & 'REGISTER (COLUMN) USED FOR HISTOGRAM',IRTFLG) IF (ICOL .LE. 0 .OR. ICOL .GE. MAXX) THEN CALL ERRT(100,'HISD',NE) GOTO 999 ENDIF NBINS = 128 CALL RDPRI1S(NBINS,NOT_USED,'NUMBER OF BINS',IRTFLG) IF (NBINS .LE. 0) THEN CALL ERRT(100,'MUST HAVE AT LEAST 1 BIN',NE) GOTO 999 ENDIF ALLOCATE(FREQ(NBINS),STAT=IRTFLG) IF (IRTFLG. NE. 0) THEN CALL ERRT(46,'HISD, NBINS',NBINS) GOTO 999 ENDIF C IF ONE OF THE FCHAR LETTERS IS 'R' ASK USER FOR THE RANGE IF (FCHAR(4:4) .EQ. 'R' .OR. FCHAR(5:5) .EQ. 'R') THEN CALL RDPRM2(HMIN,HMAX,NOT_USED,'HISTOGRAM RANGE MIN,MAX') ELSE C FIND MIN & MAX FOR HISTOGRAM BY READING DOC FILE HMAX = -1.0E23 HMIN = -HMAX DO IKEY = 1,MAXY IF (PBUF(1, IKEY) .GT. 0) THEN C KEY EXISTS BVAL = PBUF( ICOL + 1, IKEY) HMAX = AMAX1(HMAX,BVAL) HMIN = AMIN1(HMIN,BVAL) ENDIF ENDDO ENDIF C ZERO THE HISTOGRAM FREQ = 0.0 C CALCULATE HISTOGRAM HERE HDIFF = HMAX - HMIN FF = (NBINS - 1.0) / HDIFF BSIZ = 1.0 / FF FNUMEL = 0 FKEYS = 0.0 Y = 0.0 HAV = 0.0 HAV2 = 0.0 DO IKEY = 1,MAXY IF (PBUF( 1, IKEY) .GT. 0) THEN C KEY EXISTS BVAL = PBUF( ICOL + 1, IKEY) JBIN = INT((BVAL - HMIN) * FF) + 1.5 FKEYS = FKEYS + 1 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 NDEV = NDAT IF (FCHAR(4:4) .EQ. 'D' .OR. FCHAR(5:5) .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(NDOC,PLIST,3,IRTFLG) ENDDO CALL SAVDC CLOSE(NDOC) ELSE C OUTPUT IN RESULTS FILE WRITE(NDEV,*) ' ' CALL GRAPHS(NDEV,FREQ,NBINS,1,1,1.0,IRTFLG) WRITE(NDEV,*) ' ' IF (IRTFLG .NE. 0) THEN CALL ERRT(100,'HISD',NE) GOTO 999 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 .OR. MAXBIN .EQ. NBINS) THEN BMODE = MAXBIN - 1 ELSE YM1 = FREQ(MAXBIN-1) Y1 = FREQ(MAXBIN+1) BMODE = FLOAT(MAXBIN-1) + & (YM1-Y1) *.5 / (YM1+Y1-2.* HISMAX) ENDIF HMODE = HMIN + BMODE * BSIZ DTOP = HAV2 - HAV * HAV / DBLE(FNUMEL) IF (DTOP .LT. 0.0) THEN C SQRT OF NEGATIVE NUMBER WRITE(NOUT,*) '*** HISD, SQRT(',DTOP,') IMPOSSIBLE' CALL ERRT(100,'HIST',NE) GOTO 999 ELSEIF (FNUMEL .EQ. 1) THEN WRITE(NOUT,*) '*** CAN NOT DETERMINE S.D. (ONLY ONE PIXEL) ' CALL ERRT(100,'HISD',NE) GOTO 999 ENDIF HSIG = DSQRT(DTOP / DBLE(FNUMEL - 1)) HAV = HAV / DBLE(FNUMEL) FNBINS = NBINS WRITE(NDEV,*) ' ' WRITE(NDEV,90) HMIN,HMAX,FKEYS,FNUMEL, & FNBINS,BSIZ,HAV,HMODE,HSIG 90 FORMAT(' HISTOGRAM RANGE: ',G11.4,' ......... ',G11.4,/, & ' FILE KEYS: ',G11.4,' HIST. KEYS: ',G11.4,/, & ' NO. OF BINS: ',G11.4,' BIN SIZE: ',G11.4,/, & ' HIST. MEAN: ',G11.4,' HIST. MODE: ',G11.4,/, & ' HIST. S.D.: ',G11.4) WRITE(NDEV,*) ' ' IF (NDEV .NE. NOUT) WRITE(NOUT,*) ' ' IF (NDEV .NE. NOUT) WRITE(NOUT,90) HMIN,HMAX,FKEYS,FNUMEL, & FNBINS,BSIZ,HAV,HMODE,HSIG IF (NDEV .NE. NOUT) WRITE(NOUT,*) ' ' 999 IF (ASSOCIATED(PBUF)) DEALLOCATE(PBUF) IF (ALLOCATED(FREQ)) DEALLOCATE(FREQ) END