/************************************************************************** * * * Copyright (c) 1996 Michael Richmond and Richard Treffers * * * * This software may be copied and distributed for * * educational, research and not for profit services * * provided that this copyright and statement are * * included in all such copies. * * * **************************************************************************/ /* computes histogram HIST filename [BIN=] [BOX=] [PIX=] [MIN=] [MAX=] [QUIET] [help] writes filename.HIS 9/22/92 - usage, exit(0) -rrt 9/22/94 - add circle - combined with histlib 1/23/96 - added "help" option. MWR 4/18/97 - changed some "long" vars to "double" in makehist in fact, probably should check over all code in this file to find and remove assumptions about integers and maximum possible values. MWR */ #include "pcvista.h" #include "fits.h" #include #include #include #include #include int main(int,char **); static void lookhist(int); static void writehist(char *, int); static int makehist(FITS_HANDLE); static char *warning(); /* this holds # of pixels with values in each bin */ static long hist[NHIST]; /* this holds the index into "hist" for each possible pixel value */ static int hist_index[(MAX_DATA_VAL - MIN_DATA_VAL) + 1]; static long hist_under; /* histogram underflow number */ static long hist_over; /* overflow counter */ static long hist_pix; static int binsize = 1; static int hist_min = 0, hist_max = MAX_DATA_VAL; static int circle_flag=0; static int cr,cc,rad; static int quiet; static int nskip=1; static int nr, nc, sr, sc, nrow, ncol; double total_counts; #define USAGE \ "usage: hist filename [box=] [circle=] [bin=] [min=] [max=] [pix=] [quiet] [help]" int main(argc, argv) int argc; char *argv[]; { FITS_HANDLE fh; int i, boxno, circno; int topbin; char *gotit; if (argc == 1) { error(-1, USAGE); } if (strcmp(argv[1], "help") == 0) { printf("%s\n", USAGE); exit(0); } fh = fits_open(argv[1], "r", &nrow, &ncol); nr = nrow; nc = ncol; for (i = 2; i < argc; i++) { if (keyword("help", argv[i])) { printf("%s\n", USAGE); exit(0); } if ((gotit = find("bin", argv[i])) != NULL) { binsize = (int)(evaluate(gotit) + 0.5); continue; } if ((gotit = find("box", argv[i])) != NULL) { boxno = (int)(evaluate(gotit) + 0.5); if (getbox(boxno, &sr, &sc, &nr, &nc)) error(-1, "can't find box"); check_box(boxno, nrow, ncol); continue; } if ((gotit = find("circle", argv[i])) != NULL) { circno = (int)(evaluate(gotit) + SMALL); if (getcircle(circno, &cr, &cc, &rad)) error(-1, "mn: circle not found"); check_circle(circno, nrow, ncol); circle_flag=1; continue; } if ((gotit = find("min", argv[i])) != NULL) { hist_min = (unsigned)(evaluate(gotit) + 0.5); continue; } if ((gotit = find("max", argv[i])) != NULL) { hist_max = (unsigned)(evaluate(gotit) + 0.5); continue; } if ((gotit = find("pix",argv[i])) != NULL) { nskip = (int)(evaluate(gotit) + SMALL); continue; } if (keyword("quiet", argv[i])) { quiet = 1; continue; } fprintf(stderr,"hist: unknown argument '%s'\n", argv[i]); exit(1); } /* give a warning message if the bin size may be too small */ if ((hist_max - hist_min)/binsize > NHIST) { fprintf(stderr, "hist warning: any binsize less than %d may be too small\n", ((hist_max - hist_min)/NHIST) + 1); } topbin = makehist(fh); lookhist(topbin); writehist(argv[1], topbin); exit(0); } static void lookhist(topbin) int topbin; { int i, ncountmax, median, n1; long hmax, sum, half, s1; ncountmax = 0; /* determine mode -- most popular bin */ total_counts = 0.0; for (i = 0, hmax = 0; i <= topbin; i++) { total_counts += hist[i]; if (hist[i] > hmax) { hmax = hist[i]; ncountmax = i*binsize + hist_min; } } half = total_counts/2; /* * now determine median, by walking up bins * until we reach half of all pixels */ for (sum = 0, median = 0; sum < half; median++) sum += hist[median]; n1 = hist_min + binsize*(median - 2); s1 = sum - hist[median - 1]; if ((sum - s1) > 0) { median = n1 + (int)((double)binsize*((half - s1)/(sum - s1))); } else { median = 0; } if (!quiet) { printf("Median %10d adus %s\n", median,warning()); printf("Mode %10d adus %s\n", ncountmax,warning()); } } /* write the histogram into a disk file */ /* skip any set of more than 2 consecutive entries for 0 pixels */ static void writehist(name, topbin) char *name; int topbin; { int i, val; FILE *fout; char fname[NBUF]; long sum = 0; strcpy(fname, name); exten(fname, ".his"); if ((fout = fopen(fname, "w")) == NULL) error(-1, "can't open output file"); if (total_counts == 0.0) total_counts = 1.0; for (i = 0; i <= topbin; i++) { if ((i > 0) && (i < topbin)) { if ((hist[i - 1] == 0) && (hist[i] == 0) && (hist[i + 1] == 0)) { continue; } } val = i*binsize + hist_min; sum += hist[i]; if (fprintf(fout, "%8d %16ld %5.3f\n", val, hist[i], sum/total_counts) == 0) { fprintf(stderr, "hist: error writing histogram to disk file\n"); break; } } fclose(fout); } /* fill up the array 'hist' with a histogram of data values for the given file. return the number of the highest bin used, up to a maximum of NHIST - 1. */ static int makehist(fh) FITS_HANDLE fh; { int index; int number, dmin = MAX_DATA_VAL, dmax = MIN_DATA_VAL; static int16 data[NMAX]; int i, row, cmin, cmax, rmin, rmax, bin, topbin, pixval; double total, sumsq, pixnum; double mean, diff, stdev; total = 0; sumsq = 0; topbin = 1; pixnum = 0; hist_under=hist_over=0; /* placate compiler */ rmax = -1; rmin = -1; cmax = -1; cmin = -1; /* initialize the histogram lookup table */ for (pixval = MIN_DATA_VAL; pixval <= MAX_DATA_VAL; pixval++) { index = pixval - MIN_DATA_VAL; if ((bin = (pixval - hist_min)/binsize) < 0) { hist_index[index] = 0; } else if (bin >= NHIST) { hist_index[index] = NHIST - 1; /* we can just set all the rest to (NHIST - 1), too */ for ( ; index < NHIST; index++) { hist_index[index] = NHIST - 1; } } else { hist_index[index] = bin; } } if(circle_flag){ sr=cr-rad; nr=2*rad; } for (row = sr; row < (sr + nr); row++) { if(circle_flag){ int chord; chord=find_chord(row-cr,rad); nc=2*chord; sc=cc-chord; } fits_get_data(fh, row, sc, data, nc); for (i = 0; i < nc; i += nskip) { number = data[i]; if (number > dmax) { dmax = number; cmax = i + sc; rmax = row; } if (number < dmin) { dmin = number; cmin = i + sc; rmin = row; } total += number; sumsq += number*number; pixnum++; /* don't use numbers less than the minimum, or more than the maximum, in making the histogram */ if ((number < hist_min)){ hist_under++; continue; } if( (number > hist_max)){ hist_over++; continue; } if ((index = hist_index[number - MIN_DATA_VAL]) > topbin) { topbin = index; } hist[index]++; } } /* figure out what the top bin used was */ if (topbin >= NHIST) topbin = NHIST - 1; mean = (double) total/pixnum; hist_pix = pixnum; diff = (sumsq - pixnum*mean*mean); if (diff > 0.0) stdev = sqrt(diff/(pixnum - 1.0)); else stdev = 0.0; if (!quiet) { printf("Mean adus %8.2f\n", mean); printf("Total adus %10f\n", total); printf("Max. adus %6d at row= %3d col=%3d \n", dmax, rmax, cmax); printf("Min. adus %6d at row= %3d col=%3d\n", dmin, rmin, cmin); printf("Total pixels %10f\n", pixnum); printf("Overflow pixels %10ld\n", hist_over); printf("Underflow pixels %10ld\n", hist_under); printf("Standard deviation = %10.2f\n", stdev); } return(topbin); } /* returns warning if there is overflow and/or undeflow */ static char *warning() { if(hist_under) return " *** warning underflow"; if(hist_over) return " *** warning overflow"; return ""; }