Update to current Development release - 10.66.00

git-svn-id: http://svn.code.sf.net/p/netpbm/code/advanced@2172 9d0c8265-081b-0410-96cb-a4ca84ce46f8

2 files changed, 502 insertions, 140 deletions

diff --git a/analyzer/pgmhist.c b/analyzer/pgmhist.c
index f3a67383..1e779655 100644
--- a/analyzer/pgmhist.c
+++ b/analyzer/pgmhist.c
@@ -20,7 +20,7 @@
 
 
 
-struct cmdline_info {
+struct CmdlineInfo {
     /* All the information the user supplied in the command line,
        in a form easy for the program to use.
     */
@@ -29,13 +29,14 @@ struct cmdline_info {
     unsigned int median;
     unsigned int quartile;
     unsigned int decile;
+    unsigned int forensic;
 };
 
 
 
 static void
 parseCommandLine(int argc, const char ** argv,
-                 struct cmdline_info * const cmdlineP) {
+                 struct CmdlineInfo * const cmdlineP) {
 /*----------------------------------------------------------------------------
    Note that the file spec array we return is stored in the storage that
    was passed to us as the argv array.
@@ -43,7 +44,7 @@ parseCommandLine(int argc, const char ** argv,
     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 */
@@ -56,6 +57,8 @@ parseCommandLine(int argc, const char ** argv,
             &cmdlineP->quartile,            0);
     OPTENT3(0,   "decile",        OPT_FLAG,  NULL,
             &cmdlineP->decile,              0);
+    OPTENT3(0,   "forensic",      OPT_FLAG,  NULL,
+            &cmdlineP->forensic,            0);
 
     opt.opt_table = option_def;
     opt.short_allowed = FALSE;  /* We have no short (old-fashioned) options */
@@ -68,7 +71,7 @@ parseCommandLine(int argc, const char ** argv,
         pm_error("You may specify only one of -median, -quartile, "
                  "and -decile");
 
-    if (argc-1 == 0) 
+    if (argc-1 == 0)
         cmdlineP->inputFileName = "-";
     else if (argc-1 != 1)
         pm_error("Program takes zero or one argument (filename).  You "
@@ -81,39 +84,80 @@ parseCommandLine(int argc, const char ** argv,
 
 
 
+static gray
+universalMaxval(gray const maxval,
+                int  const format) {
+/*----------------------------------------------------------------------------
+  A maxval that makes it impossible for a pixel to be invalid in an image that
+  states it maxval as 'maxval' and has format 'format'.
+
+  E.g. in a one-byte-per-sample image, it's not possible to read a sample
+  value greater than 255, so a maxval of 255 makes it impossible for a sample
+  to be invalid.
+
+  But: we never go above 65535, which means our maxval isn't entirely
+  universal.  If the image is plain PGM, it could contain a pixel that
+  exceeds even that.
+-----------------------------------------------------------------------------*/
+    assert(0 < maxval && maxval < 65536);
+
+    if (format == RPGM_FORMAT) {
+        /*  Raw PGM stream has either one or two bytes per pixel, depending
+            upon its stated maxval.
+        */
+        if (maxval > 255)
+            return 65535;
+        else
+            return 255;
+    } else if (format == RPBM_FORMAT) {
+        /* A Raw PBM stream has one bit per pixel, which libnetpbm renders
+           as 0 or 255 when we read it.
+        */
+        assert(maxval == 255);
+        return 255;
+    } else {
+        /* A plain PGM or PBM stream has essentially unlimited range in the
+           tokens that are supposed to be sample values.  We arbitrarily draw
+           the line at 65535.
+        */
+        return 65535;
+    }
+}
+
+
+
 static void
-buildHistogram(FILE *          const ifP,
-               unsigned int ** const histP,
-               gray *          const maxvalP) {
+buildHistogram(FILE *               const ifP,
+               int                  const format,
+               unsigned int         const cols,
+               unsigned int         const rows,
+               gray                 const mmaxval,
+               unsigned long int ** const histP) {
+/*----------------------------------------------------------------------------
+   Compute the histogram of sample values in the input stream *ifP as *histP,
+   in newly malloced storage.
 
+   Assume the image maxval is 'mmaxval'.  Assume *ifP is positioned to the
+   start of the raster.
+-----------------------------------------------------------------------------*/
     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);
+    unsigned long int * hist;  /* malloc'ed array */
 
     grayrow = pgm_allocrow(cols);
 
-    MALLOCARRAY(hist, maxval + 1);
+    MALLOCARRAY(hist, mmaxval + 1);
     if (hist == NULL)
         pm_error("out of memory");
 
-    for (i = 0; i <= maxval; ++i)
+    for (i = 0; i <= mmaxval; ++i)
         hist[i] = 0;
 
     for (row = 0; row < rows; ++row) {
         unsigned int col;
 
-        pgm_readpgmrow(ifP, grayrow, cols, maxval, format);
+        pgm_readpgmrow(ifP, grayrow, cols, mmaxval, format);
 
         for (col = 0; col < cols; ++col) {
             /* Because total pixels in image is limited: */
@@ -124,49 +168,33 @@ buildHistogram(FILE *          const ifP,
     }
     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;
+    *histP = hist;
 }
 
 
 
 static void
-findQuantiles(unsigned int const n,
-              unsigned int const hist[],
-              gray         const maxval,
-              gray *       const quantile) {
+findQuantiles(unsigned int      const n,
+              unsigned long int const hist[],
+              unsigned long int const totalCt,
+              gray              const mmaxval,
+              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.
+   'mmaxval' is the highest index in hist[] (so its size is 'mmaxval' + 1,
+   and there are no pixels greater than 'mmaxval' in the image).
 
    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.  
+   are less than or equal to it.
 
    quantile[] must be allocated at least to size 'n'.
 
-   We return 
+   'n' must not be more than 100.
 -----------------------------------------------------------------------------*/
-    unsigned int const totalCt = sum(hist, maxval);
-
     unsigned int quantSeq;
         /* 0 is first quantile, 1 is second quantile, etc. */
 
@@ -183,9 +211,9 @@ findQuantiles(unsigned int const n,
     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;
+        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.
@@ -195,14 +223,11 @@ findQuantiles(unsigned int const n,
           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.
+        /* at sampleVal == mmaxval, cumCt == totalCt, so because
+           quantCt <= 'totalCt', 'sampleVal' cannot go above mmaxval.
         */
 
         while (cumCt < quantCt) {
@@ -210,7 +235,7 @@ findQuantiles(unsigned int const n,
             cumCt += hist[sampleVal];
         }
 
-        assert(sampleVal <= maxval);
+        assert(sampleVal <= mmaxval);
 
         /* 'sampleVal' is the lowest sample value for which at least 'quantCt'
            pixels have that sample value or less.  'cumCt' is the number
@@ -223,9 +248,10 @@ findQuantiles(unsigned int const n,
 
 
 static void
-countCumulative(unsigned int    const hist[],
-                gray            const maxval,
-                unsigned int ** const rcountP) {
+countCumulative(unsigned long int    const hist[],
+                gray                 const mmaxval,
+                unsigned long int    const totalPixelCt,
+                unsigned long 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]
@@ -233,17 +259,17 @@ countCumulative(unsigned int    const hist[],
 
    *rcountP is newly malloced memory.
 -----------------------------------------------------------------------------*/
-    unsigned int * rcount;
-    unsigned int cumCount;
+    unsigned long int * rcount;
+    unsigned long int cumCount;
     int i;
-    
-    MALLOCARRAY(rcount, maxval + 1);
+
+    MALLOCARRAY(rcount, mmaxval + 1);
     if (rcount == NULL)
         pm_error("out of memory");
 
-    for (i = maxval, cumCount = 0; i >= 0; --i) {
+    for (i = mmaxval, cumCount = 0; i >= 0; --i) {
         /* Because total pixels in image is limited: */
-        assert(UINT_MAX - hist[i] >= cumCount);
+        assert(ULONG_MAX - hist[i] >= cumCount);
 
         cumCount += hist[i];
         rcount[i] = cumCount;
@@ -255,11 +281,11 @@ countCumulative(unsigned int    const hist[],
 
 
 static void
-reportHistHumanFriendly(unsigned int const hist[],
-                        unsigned int const rcount[],
-                        gray         const maxval) {
+reportHistHumanFriendly(unsigned long int const hist[],
+                        unsigned long int const rcount[],
+                        gray              const maxval) {
 
-    unsigned int const totalPixels = rcount[0];
+    unsigned long int const totalPixelCt = rcount[0];
 
     unsigned int cumCount;
     unsigned int i;
@@ -271,9 +297,48 @@ reportHistHumanFriendly(unsigned int const hist[],
         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);
+                "%5d  %5ld  %5.3g%%  %5.3g%%\n", i, hist[i],
+                (float) cumCount * 100.0 / totalPixelCt,
+                (float) rcount[i] * 100.0 / totalPixelCt);
+        }
+    }
+}
+
+
+static void
+reportHistForensicHumanFriendly(unsigned long int const hist[],
+                                unsigned long int const rcount[],
+                                gray              const maxval,
+                                gray              const mmaxval) {
+
+    unsigned long int const totalPixelCt = rcount[0];
+
+    unsigned long 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  %5ld  %5.3g%%  %5.3g%%\n", i, hist[i],
+                (float) cumCount * 100.0 / totalPixelCt,
+                (float) rcount[i] * 100.0 / totalPixelCt);
+        }
+    }
+    if (totalPixelCt > cumCount) {
+        printf("-----  -----\n");
+
+        for (i = maxval; i <= mmaxval; ++i) {
+            if (hist[i] > 0) {
+                cumCount += hist[i];
+                printf(
+                    "%5d  %5ld  %5.3g%%  %5.3g%%\n", i, hist[i],
+                    (float) cumCount * 100.0 / totalPixelCt,
+                    (float) rcount[i] * 100.0 / totalPixelCt);
+            }
         }
     }
 }
@@ -281,13 +346,13 @@ reportHistHumanFriendly(unsigned int const hist[],
 
 
 static void
-reportHistMachineFriendly(unsigned int const hist[],
-                          gray         const maxval) {
+reportHistMachineFriendly(unsigned long int const hist[],
+                          gray              const maxval) {
 
     unsigned int i;
 
     for (i = 0; i <= maxval; ++i) {
-        printf("%u %u\n", i, hist[i]);
+        printf("%u %lu\n", i, hist[i]);
     }
 }
 
@@ -345,60 +410,132 @@ reportDecilesHumanFriendly(gray const decile[]) {
 
 
 
-int
-main(int argc, const char ** argv) {
+static void
+summarizeInvalidPixels(unsigned long int const hist[],
+                       unsigned long int const rcount[],
+                       gray              const mmaxval,
+                       gray              const maxval) {
+/*----------------------------------------------------------------------------
+  Print total count of valid and invalid pixels, if there are any
+  invalid ones.
+-----------------------------------------------------------------------------*/
+    unsigned long int const invalidPixelCt = 
+        mmaxval > maxval ? rcount[maxval+1] : 0;
+
+    if (invalidPixelCt > 0) {
+        unsigned long int const totalPixelCt = rcount[0];
+        unsigned long int const validPixelCt = totalPixelCt - invalidPixelCt;
+
+        printf("\n");
+        printf("** Image stream contains invalid sample values "
+               "(above maxval %u)\n", maxval);
+        printf("Valid sample values:   %lu (%5.4g%%)\n",
+               validPixelCt,   (float)validPixelCt / totalPixelCt * 100.0);
+        printf("Invalid sample values: %lu (%5.4g%%)\n",
+               invalidPixelCt, (float)invalidPixelCt / totalPixelCt * 100.0);
+    }
+}
 
-    struct cmdline_info cmdline;
-    FILE * ifP;
-    gray maxval;
-    unsigned int * hist;   /* malloc'ed array */
 
-    pm_proginit(&argc, argv);
 
-    parseCommandLine(argc, argv, &cmdline);
+static void
+reportFromHistogram(const unsigned long int * const hist,
+                    gray                      const mmaxval,
+                    gray                      const maxval,
+                    unsigned long int         const totalPixelCt,
+                    struct CmdlineInfo        const cmdline) {
+/*----------------------------------------------------------------------------
+   Analyze histogram 'hist', which has 'mmaxval' buckets, and report
+   what we find.
 
-    ifP = pm_openr(cmdline.inputFileName);
+   'maxval' is the maxval that the image states (but note that we tolerate
+   invalid sample values greater than maxval, which could be as high as
+   'mmaxval').
 
-    buildHistogram(ifP, &hist, &maxval);
+   'cmdline' tells what kind of reporting to do.
+-----------------------------------------------------------------------------*/
 
     if (cmdline.median) {
         gray median[2];
-        findQuantiles(2, hist, maxval, median); 
+        findQuantiles(2, hist, totalPixelCt, mmaxval, median);
         if (cmdline.machine)
             reportQuantilesMachineFriendly(median, 1);
         else
             reportMedianHumanFriendly(median[0]);
     } else if (cmdline.quartile) {
         gray quartile[4];
-        findQuantiles(4, hist, maxval, quartile);
+        findQuantiles(4, hist, totalPixelCt, mmaxval, quartile);
         if (cmdline.machine)
             reportQuantilesMachineFriendly(quartile, 4);
         else
             reportQuartilesHumanFriendly(quartile);
     } else if (cmdline.decile) {
         gray decile[10];
-        findQuantiles(10, hist, maxval, decile);
+        findQuantiles(10, hist, totalPixelCt, mmaxval, decile);
         if (cmdline.machine)
             reportQuantilesMachineFriendly(decile, 10);
         else
             reportDecilesHumanFriendly(decile);
     } else {
         if (cmdline.machine)
-            reportHistMachineFriendly(hist, maxval);
+            reportHistMachineFriendly(hist, mmaxval);
         else {
-            unsigned int * rcount; /* malloc'ed array */
-            countCumulative(hist, maxval, &rcount);
-            reportHistHumanFriendly(hist, rcount, maxval);
+            unsigned long int * rcount; /* malloc'ed array */
+            countCumulative(hist, mmaxval, totalPixelCt, &rcount);
+            if (cmdline.forensic)
+                reportHistForensicHumanFriendly(hist, rcount, maxval, mmaxval);
+            else
+                reportHistHumanFriendly(hist, rcount, maxval);
+
+            summarizeInvalidPixels(hist, rcount, mmaxval, maxval);
 
             free(rcount);
         }
     }
+}
+
+
+
+int
+main(int argc, const char ** argv) {
+
+    struct CmdlineInfo cmdline;
+    FILE * ifP;
+    int rows, cols;
+    int format;
+    gray maxval;
+        /* Stated maxval of the image */
+    gray mmaxval;
+        /* Maxval we assume, which may be greater than the stated maxval
+           so that we can process invalid pixels in the image that exceed
+           the maxval.
+        */
+    unsigned long int totalPixelCt;
+    unsigned long int * hist;   /* malloc'ed array */
+
+    pm_proginit(&argc, argv);
+
+    parseCommandLine(argc, argv, &cmdline);
+
+    ifP = pm_openr(cmdline.inputFileName);
+
+    pgm_readpgminit(ifP, &cols, &rows, &maxval, &format);
+
+    if (ULONG_MAX / cols < rows)
+        pm_error("Too many pixels (%u x %u) in image.  "
+                 "Maximum computable is %lu",
+                 cols, rows, ULONG_MAX);
+
+    totalPixelCt = cols * rows;
+
+    mmaxval = cmdline.forensic ? universalMaxval(maxval, format) : maxval;
+
+    buildHistogram(ifP, format, cols, rows, mmaxval, &hist);
+
+    reportFromHistogram(hist, mmaxval, maxval, totalPixelCt, cmdline);
 
     free(hist);
     pm_close(ifP);
 
     return 0;
 }
-
-
-
diff --git a/analyzer/ppmhist.c b/analyzer/ppmhist.c
index 2f6c9348..cc47bb82 100644
--- a/analyzer/ppmhist.c
+++ b/analyzer/ppmhist.c
@@ -20,17 +20,18 @@
 
 enum sort {SORT_BY_FREQUENCY, SORT_BY_RGB};
 
-enum colorFmt {FMT_DECIMAL, FMT_HEX, FMT_FLOAT, FMT_PPMPLAIN};
+enum ColorFmt {FMT_DECIMAL, FMT_HEX, FMT_FLOAT, FMT_PPMPLAIN};
 
 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;  /* Name of input file */
-    unsigned int noheader;    /* -noheader option */
-    enum colorFmt colorFmt;
-    unsigned int colorname;   /* -colorname option */
-    enum sort sort;           /* -sort option */
+    unsigned int noheader;
+    enum ColorFmt colorFmt;
+    unsigned int colorname;
+    enum sort sort;
+    unsigned int forensic;
 };
 
 
@@ -45,7 +46,7 @@ parseCommandLine(int argc, const char ** argv,
     optStruct3 opt;  /* set by OPTENT3 */
     optEntry * option_def;
     unsigned int option_def_index;
-    
+
     unsigned int hexcolorOpt, floatOpt, mapOpt, nomapOpt;
     const char * sort_type;
 
@@ -59,6 +60,7 @@ parseCommandLine(int argc, const char ** argv,
     OPTENT3(0,   "float",     OPT_FLAG, NULL,  &floatOpt,              0);
     OPTENT3(0,   "colorname", OPT_FLAG, NULL,  &cmdlineP->colorname,   0);
     OPTENT3(0,   "sort",      OPT_STRING, &sort_type, NULL,            0);
+    OPTENT3(0,   "forensic",  OPT_FLAG, NULL,  &cmdlineP->forensic,    0);
 
     opt.opt_table = option_def;
     opt.short_allowed = FALSE;  /* We have no short (old-fashioned) options */
@@ -69,8 +71,9 @@ parseCommandLine(int argc, const char ** argv,
 
     pm_optParseOptions3(&argc, (char **)argv, opt, sizeof(opt), 0);
         /* Uses and sets argc, argv, and some of *cmdlineP and others. */
+    free(option_def);
 
-    if (argc-1 == 0) 
+    if (argc-1 == 0)
         cmdlineP->inputFileName = "-";
     else if (argc-1 != 1)
         pm_error("Program takes zero or one argument (filename).  You "
@@ -84,9 +87,12 @@ parseCommandLine(int argc, const char ** argv,
         cmdlineP->colorFmt = FMT_HEX;
     else if (floatOpt)
         cmdlineP->colorFmt = FMT_FLOAT;
-    else if (mapOpt)
+    else if (mapOpt) {
+        if (cmdlineP->forensic)
+            pm_error("You cannot specify -map and -forensic together");
+
         cmdlineP->colorFmt = FMT_PPMPLAIN;
-    else 
+    } else
         cmdlineP->colorFmt = FMT_DECIMAL;
 
     if (strcmp(sort_type, "frequency") == 0)
@@ -161,8 +167,137 @@ rgbcompare(const void * const a,
 
 
 
+static pixval
+universalMaxval(pixval const maxval,
+                int    const format) {
+/*----------------------------------------------------------------------------
+  A maxval that makes it impossible for a pixel to be invalid in an image that
+  states it maxval as 'maxval' and has format 'format'.
+
+  E.g. in a one-byte-per-sample image, it's not possible to read a sample
+  value greater than 255, so a maxval of 255 makes it impossible for a sample
+  to be invalid.
+
+  But: we never go above 65535, which means our maxval isn't entirely
+  universal.  If the image is plain PPM, it could contain a pixel that
+  exceeds even that.
+-----------------------------------------------------------------------------*/
+    assert(0 < maxval && maxval < 65536);
+
+    if (format == RPPM_FORMAT || format == RPGM_FORMAT) {
+        /*  Raw PPM stream has either one or two bytes per pixel, depending
+            upon its stated maxval.
+        */
+        if (maxval > 255)
+            return 65535;
+        else
+            return 255;
+    } else if (format == RPBM_FORMAT) {
+        /* A Raw PBM stream has one bit per pixel, which libnetpbm renders
+           as 0 or 255 when we read it.
+        */
+        assert(maxval == 255);
+        return 255;
+    } else {
+        /* A plain PPM stream has essentially unlimited range in the
+           tokens that are supposed to be sample values.  We arbitrarily draw
+           the line at 65535.
+        */
+        return 65535;
+    }
+}
+
+
+
+static bool
+colorIsValid(pixel  const color,
+             pixval const maxval) {
+
+    pixval const r = PPM_GETR(color);
+    pixval const g = PPM_GETG(color);
+    pixval const b = PPM_GETB(color);
+
+    return r <= maxval && g <= maxval && b <= maxval;
+}
+
+
+
+static void
+separateInvalidItems(colorhist_vector const chv,
+                     colorhist_vector const chvInvalid,
+                     pixval           const maxval,
+                     unsigned int     const colorCt,
+                     unsigned int  *  const validColorCtP) {
+/*----------------------------------------------------------------------------
+  Move invalid color entries from chv to chvInvalid.
+  Count how many color entries are valid. 
+-----------------------------------------------------------------------------*/
+    unsigned int i;
+    unsigned int validCt;
+    unsigned int invalidCt;
+ 
+    for (i = 0, validCt = 0, invalidCt = 0; i < colorCt; ++i) {
+        if (!colorIsValid(chv[i].color, maxval))
+            chvInvalid[invalidCt++] = chv[i];
+        else
+            chv[validCt++] = chv[i];
+    } 
+    *validColorCtP = validCt;
+}
+  
+
+
+static void
+sortHistogramForensic(enum sort        const sortFn,
+                      colorhist_vector const chv,
+                      colorhist_vector const chvInvalid,
+                      pixval           const maxval,
+                      unsigned int     const colorCt,
+                      unsigned int *   const validColorCtP) {
+
+    unsigned int validColorCt;
+
+    separateInvalidItems(chv, chvInvalid, maxval, colorCt, &validColorCt);
+
+    {
+        int (*compare_function)(const void *, const void *);
+    
+        switch (sortFn) {
+        case SORT_BY_FREQUENCY: compare_function = countcompare; break;
+        case SORT_BY_RGB:       compare_function = rgbcompare;   break;
+        }
+
+        qsort((void*) chv, validColorCt,
+              sizeof(struct colorhist_item), compare_function);
+        
+        qsort((void*) chvInvalid, colorCt - validColorCt,
+              sizeof(struct colorhist_item), compare_function);
+    }
+    *validColorCtP = validColorCt;
+}
+
+
+
+static void
+sortHistogramNormal(enum sort        const sortFn,
+                    colorhist_vector const chv,
+                    unsigned int     const colorCt) {
+
+    int (*compare_function)(const void *, const void *);
+
+    switch (sortFn) {
+    case SORT_BY_FREQUENCY: compare_function = countcompare; break;
+    case SORT_BY_RGB:       compare_function = rgbcompare;   break;
+    }
+
+    qsort((void*) chv, colorCt, sizeof(struct colorhist_item),
+          compare_function);
+}
+
+
+
 static const char *
-colornameLabel(pixel        const color, 
+colornameLabel(pixel        const color,
                pixval       const maxval,
                unsigned int const nDictColor,
                pixel        const dictColors[],
@@ -172,15 +307,15 @@ colornameLabel(pixel        const color,
    dictionary to it.  If the name returned is not the exact color,
    prefix it with "*".  Otherwise, prefix it with " ".
 
-   'nDictColor', dictColors[], and dictColorNames[] are the color 
+   'nDictColor', dictColors[], and dictColorNames[] are the color
    dictionary.
 
    Return the name in static storage within this subroutine.
 -----------------------------------------------------------------------------*/
     static char retval[32];
     int colorIndex;
-    
-    pixel color255;  
+
+    pixel color255;
         /* The color, normalized to a maxval of 255: the maxval of a color
            dictionary.
         */
@@ -190,31 +325,31 @@ colornameLabel(pixel        const color,
     colorIndex = ppm_findclosestcolor(dictColors, nDictColor, &color);
 
     assert(colorIndex >= 0 && colorIndex < nDictColor);
-    
+
     if (PPM_EQUAL(dictColors[colorIndex], color))
         STRSCPY(retval, " ");
     else
         STRSCPY(retval, "*");
-    
+
     STRSCAT(retval, dictColornames[colorIndex]);
-    
+
     return retval;
 }
-                
+
 
 
 static void
-printColors(colorhist_vector const chv, 
-            int              const nColors,
+printColors(colorhist_vector const chv,
+            int              const colorCt,
             pixval           const maxval,
-            enum colorFmt    const colorFmt,
+            enum ColorFmt    const colorFmt,
             unsigned int     const nKnown,
             pixel            const knownColors[],
             const char *     const colornames[]) {
 
     int i;
 
-    for (i = 0; i < nColors; i++) {
+    for (i = 0; i < colorCt; i++) {
         pixval       const r          = PPM_GETR(chv[i].color);
         pixval       const g          = PPM_GETG(chv[i].color);
         pixval       const b          = PPM_GETB(chv[i].color);
@@ -226,7 +361,7 @@ printColors(colorhist_vector const chv,
         const char * colornameValue;
 
         if (colornames)
-            colornameValue = colornameLabel(chv[i].color, maxval, 
+            colornameValue = colornameLabel(chv[i].color, maxval,
                                             nKnown, knownColors, colornames);
         else
             colornameValue = "";
@@ -257,23 +392,97 @@ printColors(colorhist_vector const chv,
 
 
 
+static void
+summarizeInvalidPixels(unsigned long int const validPixelCt,
+                       unsigned long int const invalidPixelCt,
+                       pixval            const maxval) {
+/*----------------------------------------------------------------------------
+  Print total count of valid and invalid pixels, if there are any
+  invalid ones.
+-----------------------------------------------------------------------------*/
+    if (invalidPixelCt > 0) {
+        unsigned long int const totalPixelCt = validPixelCt + invalidPixelCt;
+
+        printf("\n");
+        printf("** Image stream contains invalid sample values "
+               "(above maxval %u)\n", maxval);
+        printf("** Valid sample values : %lu (%5.4g%%)\n",
+               validPixelCt,   (float) validPixelCt   / totalPixelCt * 100.0);
+        printf("** Invalid sample values : %lu (%5.4g%%)\n",
+               invalidPixelCt, (float) invalidPixelCt / totalPixelCt * 100.0);
+    }
+}
+
+
+
+static void
+printInvalidSamples(colorhist_vector const chv,
+                    colorhist_vector const chvInvalid,
+                    int              const totalColorCt,
+                    unsigned int     const validColorCt,
+                    pixval           const maxval,
+                    enum ColorFmt    const colorFmt) {
+
+    unsigned int const invalidColorCt = totalColorCt - validColorCt;
+
+    unsigned int i;
+    unsigned long int validPixelCt;
+    unsigned long int invalidPixelCt;
+
+    for (i = 0, validPixelCt = 0; i < validColorCt; ++i)
+        validPixelCt += chv[i].value; 
+
+    for (i = 0, invalidPixelCt = 0; i < invalidColorCt; ++i) {
+        pixval       const r     = PPM_GETR(chvInvalid[i].color);
+        pixval       const g     = PPM_GETG(chvInvalid[i].color);
+        pixval       const b     = PPM_GETB(chvInvalid[i].color);
+        unsigned int const count = chvInvalid[i].value;
+
+        invalidPixelCt += chvInvalid[i].value; 
+
+        switch(colorFmt) {
+        case FMT_FLOAT:
+            printf(" %1.3f %1.3f %1.3f\t\t%7d\n",
+                   (double)r / maxval,
+                   (double)g / maxval,
+                   (double)b / maxval,
+                   count);
+            break;
+        case FMT_HEX:
+            printf("  %04x  %04x  %04x\t\t%7d\n",
+                   r, g, b, count);
+            break;
+        case FMT_DECIMAL:
+            printf(" %5d %5d %5d\t\t%7d\n",
+                   r, g, b, count);
+            break;
+        case FMT_PPMPLAIN:
+            assert(false);
+            break;
+        }
+    }
+
+    summarizeInvalidPixels(validPixelCt, invalidPixelCt, maxval);
+}
+
+
+
 int
 main(int argc, const char *argv[]) {
 
     struct cmdline_info cmdline;
     FILE * ifP;
     colorhist_vector chv;
+    colorhist_vector chvInvalid;
     int rows, cols;
     pixval maxval;
+    pixval mmaxval;
     int format;
-    int nColors;
-    int (*compare_function)(const void *, const void *);
-        /* The compare function to be used with qsort() to sort the
-           histogram for output
-        */
-    unsigned int nDictColor;
+    int colorCt;
+    unsigned int dictColorCt;
     const char ** dictColornames;
     pixel * dictColors;
+    unsigned int validColorCt;
 
     pm_proginit(&argc, argv);
 
@@ -283,23 +492,29 @@ main(int argc, const char *argv[]) {
 
     ppm_readppminit(ifP, &cols, &rows, &maxval, &format);
 
-    chv = ppm_computecolorhist2(ifP, cols, rows, maxval, format, 0, &nColors);
+    mmaxval = cmdline.forensic ? universalMaxval(maxval, format) : maxval;
+
+    chv = ppm_computecolorhist2(ifP, cols, rows, mmaxval, format, 0, &colorCt);
 
     pm_close(ifP);
 
-    switch (cmdline.sort) {
-    case SORT_BY_FREQUENCY:
-        compare_function = countcompare; break;
-    case SORT_BY_RGB:
-        compare_function = rgbcompare; break;
-    }
+    /* Sort and produce histogram. */
+    if (cmdline.forensic) {
+        MALLOCARRAY(chvInvalid, colorCt);
+        if (chvInvalid == NULL)
+            pm_error("out of memory generating histogram");
 
-    qsort((char*) chv, nColors, sizeof(struct colorhist_item), 
-          compare_function);
+        sortHistogramForensic(cmdline.sort, chv, chvInvalid,
+                              maxval, colorCt, &validColorCt);
+    } else {
+        chvInvalid = NULL;
+        sortHistogramNormal(cmdline.sort, chv, colorCt);
+        validColorCt = colorCt;
+    }
 
     /* And print the histogram. */
-    if (cmdline.colorFmt == FMT_PPMPLAIN) 
-        printf("P3\n# color map\n%d 1\n%d\n", nColors, maxval);
+    if (cmdline.colorFmt == FMT_PPMPLAIN)
+        printf("P3\n# color map\n%d 1\n%d\n", colorCt, maxval);
 
     if (!cmdline.noheader) {
         const char commentDelim = cmdline.colorFmt == FMT_PPMPLAIN ? '#' : ' ';
@@ -309,16 +524,20 @@ main(int argc, const char *argv[]) {
                commentDelim, cmdline.colorname ? "----" : "");
     }
     if (cmdline.colorname) {
-        bool mustOpenTrue = TRUE;
-        ppm_readcolordict(NULL, mustOpenTrue, 
-                          &nDictColor, &dictColornames, &dictColors, NULL);
+        bool const mustOpenTrue = TRUE;
+        ppm_readcolordict(NULL, mustOpenTrue,
+                          &dictColorCt, &dictColornames, &dictColors, NULL);
     } else {
         dictColors = NULL;
         dictColornames = NULL;
     }
-        
-    printColors(chv, nColors, maxval,
-                cmdline.colorFmt, nDictColor, dictColors, dictColornames);
+
+    printColors(chv, validColorCt, maxval,
+                cmdline.colorFmt, dictColorCt, dictColors, dictColornames);
+
+    if (colorCt > validColorCt)
+        printInvalidSamples(chv, chvInvalid, colorCt, validColorCt,
+                            maxval, cmdline.colorFmt);
 
     if (dictColors)
         free(dictColors);
@@ -327,5 +546,11 @@ main(int argc, const char *argv[]) {
 
     ppm_freecolorhist(chv);
 
+    if (chvInvalid)
+        ppm_freecolorhist(chvInvalid);
+
     return 0;
 }
+
+
+