/* pgmhist.c - print a histogram of the values in a PGM image ** ** Copyright (C) 1989 by Jef Poskanzer. ** ** Permission to use, copy, modify, and distribute this software and its ** documentation for any purpose and without fee is hereby granted, provided ** that the above copyright notice appear in all copies and that both that ** copyright notice and this permission notice appear in supporting ** documentation. This software is provided "as is" without express or ** implied warranty. */ #include #include #include "pm_c_util.h" #include "mallocvar.h" #include "shhopt.h" #include "pgm.h" struct cmdline_info { /* All the information the user supplied in the command line, in a form easy for the program to use. */ const char * inputFileName; /* Filename of input files */ unsigned int machine; unsigned int median; unsigned int quartile; unsigned int decile; }; static void parseCommandLine(int argc, const char ** argv, struct cmdline_info * const cmdlineP) { /*---------------------------------------------------------------------------- Note that the file spec array we return is stored in the storage that was passed to us as the argv array. -----------------------------------------------------------------------------*/ optStruct3 opt; /* set by OPTENT3 */ optEntry * option_def; unsigned int option_def_index; MALLOCARRAY_NOFAIL(option_def, 100); option_def_index = 0; /* incremented by OPTENT3 */ OPTENT3(0, "machine", OPT_FLAG, NULL, &cmdlineP->machine, 0); OPTENT3(0, "median", OPT_FLAG, NULL, &cmdlineP->median, 0); OPTENT3(0, "quartile", OPT_FLAG, NULL, &cmdlineP->quartile, 0); OPTENT3(0, "decile", OPT_FLAG, NULL, &cmdlineP->decile, 0); opt.opt_table = option_def; opt.short_allowed = FALSE; /* We have no short (old-fashioned) options */ opt.allowNegNum = FALSE; /* We have no parms that are negative numbers */ pm_optParseOptions3(&argc, (char **)argv, opt, sizeof(opt), 0); /* Uses and sets argc, argv, and some of *cmdlineP and others. */ if (cmdlineP->median + cmdlineP->quartile + cmdlineP->decile > 1) pm_error("You may specify only one of -median, -quartile, " "and -decile"); if (argc-1 == 0) cmdlineP->inputFileName = "-"; else if (argc-1 != 1) pm_error("Program takes zero or one argument (filename). You " "specified %d", argc-1); else cmdlineP->inputFileName = argv[1]; free(option_def); } static void buildHistogram(FILE * const ifP, unsigned int ** const histP, gray * const maxvalP) { gray * grayrow; int rows, cols; int format; unsigned int row; unsigned int i; unsigned int * hist; /* malloc'ed array */ gray maxval; pgm_readpgminit(ifP, &cols, &rows, &maxval, &format); if (UINT_MAX / cols < rows) pm_error("Too many pixels (%u x %u) in image. " "Maximum computable is %u", cols, rows, UINT_MAX); grayrow = pgm_allocrow(cols); MALLOCARRAY(hist, maxval + 1); if (hist == NULL) pm_error("out of memory"); for (i = 0; i <= maxval; ++i) hist[i] = 0; for (row = 0; row < rows; ++row) { unsigned int col; pgm_readpgmrow(ifP, grayrow, cols, maxval, format); for (col = 0; col < cols; ++col) { /* Because total pixels in image is limited: */ assert (hist[grayrow[col]] < INT_MAX); ++hist[grayrow[col]]; } } pgm_freerow(grayrow); *histP = hist; *maxvalP = maxval; } static unsigned int sum(unsigned int const hist[], gray const maxval) { unsigned int sum; unsigned int sampleVal; for (sampleVal = 0, sum = 0; sampleVal <= maxval; ++sampleVal) sum += hist[sampleVal]; return sum; } static void findQuantiles(unsigned int const n, unsigned int const hist[], gray const maxval, gray * const quantile) { /*---------------------------------------------------------------------------- Find the order-n quantiles (e.g. n == 4 means quartiles) of the pixel sample values, given that hist[] is the histogram of them (hist[N] is the number of pixels that have sample value N). 'maxval' is the maxval of the image, so the size of hist[] is 'maxval' + 1. We return the ith quantile as quantile[i]. For example, for quartiles, quantile[3] is the least sample value for which at least 3/4 of the pixels are less than or equal to it. quantile[] must be allocated at least to size 'n'. We return -----------------------------------------------------------------------------*/ unsigned int const totalCt = sum(hist, maxval); unsigned int quantSeq; /* 0 is first quantile, 1 is second quantile, etc. */ gray sampleVal; /* As we increment through all the possible sample values, this is the one we're considering now. */ unsigned int cumCt; /* The number of pixels that have sample value 'sampleVal' or less. */ assert(n > 1 && n <= 100); sampleVal = 0; /* initial value */ cumCt = hist[0]; /* initial value */ for (quantSeq = 1; quantSeq <= n; ++quantSeq) { unsigned long int const q = totalCt/n; unsigned long int const r = totalCt%n; unsigned long int const quantCt = q*quantSeq + (r*quantSeq +n -1)/n; /* This is how many pixels are (ignoring quantization) in the quantile. E.g. for the 3rd quartile, it is 3/4 of the pixels in the image. This is equivalent to (float) totalCt * quantSeq / n rounded upwards. We use the int version in spite of complexities for preventing overflow for slight innacuracies in floating point arithmetic causes problems when used as loop counter and array index. The maximum value for n is currently 10; in the future this may be brought up to 100. */ assert(quantCt <= totalCt); /* at sampleVal == maxval, cumCt == totalCt, so because quantCt <= 'totalCt', 'sampleVal' cannot go above maxval. */ while (cumCt < quantCt) { ++sampleVal; cumCt += hist[sampleVal]; } assert(sampleVal <= maxval); /* 'sampleVal' is the lowest sample value for which at least 'quantCt' pixels have that sample value or less. 'cumCt' is the number of pixels that have sample value 'sampleVal' or less. */ quantile[quantSeq-1] = sampleVal; } } static void countCumulative(unsigned int const hist[], gray const maxval, unsigned int ** const rcountP) { /*---------------------------------------------------------------------------- From the histogram hist[] (hist[N] is the number of pixels of sample value N), compute the cumulative distribution *rcountP ((*rcountP)[N] is the number of pixels of sample value N or higher). *rcountP is newly malloced memory. -----------------------------------------------------------------------------*/ unsigned int * rcount; unsigned int cumCount; int i; MALLOCARRAY(rcount, maxval + 1); if (rcount == NULL) pm_error("out of memory"); for (i = maxval, cumCount = 0; i >= 0; --i) { /* Because total pixels in image is limited: */ assert(UINT_MAX - hist[i] >= cumCount); cumCount += hist[i]; rcount[i] = cumCount; } *rcountP = rcount; } static void reportHistHumanFriendly(unsigned int const hist[], unsigned int const rcount[], gray const maxval) { unsigned int const totalPixels = rcount[0]; unsigned int cumCount; unsigned int i; printf("value count b%% w%% \n"); printf("----- ----- ------ ------\n"); for (i = 0, cumCount = 0; i <= maxval; ++i) { if (hist[i] > 0) { cumCount += hist[i]; printf( "%5d %5d %5.3g%% %5.3g%%\n", i, hist[i], (float) cumCount * 100.0 / totalPixels, (float) rcount[i] * 100.0 / totalPixels); } } } static void reportHistMachineFriendly(unsigned int const hist[], gray const maxval) { unsigned int i; for (i = 0; i <= maxval; ++i) { printf("%u %u\n", i, hist[i]); } } static void reportQuantilesMachineFriendly(gray const quantile[], unsigned int const n) { unsigned int i; for (i = 0; i < n; ++i) printf("%u\n", quantile[i]); } static void reportMedianHumanFriendly(gray const median) { printf("Median: %5u\n", median); } static void reportQuartilesHumanFriendly(gray const quartile[]) { unsigned int i; printf("Quartiles:\n"); printf("Q Value\n"); printf("---- -----\n"); for (i = 1; i <= 4; ++i) printf("%3u%% %5u\n", 25*i, quartile[i-1]); } static void reportDecilesHumanFriendly(gray const decile[]) { unsigned int i; printf("Deciles:\n"); printf("Q Value\n"); printf("--- -----\n"); for (i = 1; i <= 10; ++i) printf("%3u%% %5u\n", 10*i, decile[i-1]); } int main(int argc, const char ** argv) { struct cmdline_info cmdline; FILE * ifP; gray maxval; unsigned int * hist; /* malloc'ed array */ pm_proginit(&argc, argv); parseCommandLine(argc, argv, &cmdline); ifP = pm_openr(cmdline.inputFileName); buildHistogram(ifP, &hist, &maxval); if (cmdline.median) { gray median[2]; findQuantiles(2, hist, maxval, median); if (cmdline.machine) reportQuantilesMachineFriendly(median, 1); else reportMedianHumanFriendly(median[0]); } else if (cmdline.quartile) { gray quartile[4]; findQuantiles(4, hist, maxval, quartile); if (cmdline.machine) reportQuantilesMachineFriendly(quartile, 4); else reportQuartilesHumanFriendly(quartile); } else if (cmdline.decile) { gray decile[10]; findQuantiles(10, hist, maxval, decile); if (cmdline.machine) reportQuantilesMachineFriendly(decile, 10); else reportDecilesHumanFriendly(decile); } else { if (cmdline.machine) reportHistMachineFriendly(hist, maxval); else { unsigned int * rcount; /* malloc'ed array */ countCumulative(hist, maxval, &rcount); reportHistHumanFriendly(hist, rcount, maxval); free(rcount); } } free(hist); pm_close(ifP); return 0; }