Promote Stable to Super_stable

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

13 files changed, 2377 insertions, 1060 deletions

diff --git a/analyzer/Makefile b/analyzer/Makefile
index ed180917..c96a6d0a 100644
--- a/analyzer/Makefile
+++ b/analyzer/Makefile
@@ -14,7 +14,7 @@ include $(BUILDDIR)/config.mk
 # This package is so big, it's useful even when some parts won't 
 # build.
 
-PORTBINARIES = pamfile pamslice pamsumm pamtilt \
+PORTBINARIES = pamfile pammosaicknit pamslice pamsumm pamtilt \
                pbmminkowski pgmhist pnmhistmap ppmhist pgmminkowski
 MATHBINARIES = pamsharpmap pamsharpness pgmtexture pnmpsnr 
 
diff --git a/analyzer/pamfile.c b/analyzer/pamfile.c
index 4fa81330..a7c3c5ac 100644
--- a/analyzer/pamfile.c
+++ b/analyzer/pamfile.c
@@ -16,7 +16,7 @@
 #include "shhopt.h"
 #include "pam.h"
 
-struct cmdlineInfo {
+struct CmdlineInfo {
     /* All the information the user supplied in the command line,
        in a form easy for the program to use.
     */
@@ -30,14 +30,14 @@ struct cmdlineInfo {
 
 
 static void
-parseCommandLine(int argc, char ** argv,
-                 struct cmdlineInfo * const cmdlineP) {
+parseCommandLine(int argc, const char ** argv,
+                 struct CmdlineInfo * const cmdlineP) {
 /*----------------------------------------------------------------------------
    Note that the file spec array we return is stored in the storage that
    was passed to as as the argv array.
 -----------------------------------------------------------------------------*/
     optEntry * option_def;
-        /* Instructions to optParseOptions3 on how to parse our options.
+        /* Instructions to pm_optParseOptions3 on how to parse our options.
          */
     optStruct3 opt;
 
@@ -54,11 +54,13 @@ parseCommandLine(int argc, char ** argv,
     opt.short_allowed = FALSE; /* We have no short (old-fashioned) options */
     opt.allowNegNum   = FALSE; /* We have no parms that are negative numbers */
     
-    optParseOptions3(&argc, argv, opt, sizeof(opt), 0);
+    pm_optParseOptions3(&argc, (char **)argv, opt, sizeof(opt), 0);
         /* Uses and sets argc, argv, and some of *cmdlineP and others */
 
     cmdlineP->inputFilespec = (const char **)&argv[1];
     cmdlineP->inputFileCount = argc - 1;
+
+    free(option_def);
 }
 
 
@@ -157,7 +159,7 @@ doOneImage(const char * const name,
         if (wantComments)
             dumpComments(comments);
     }
-    strfree(comments);
+    pm_strfree(comments);
 
     pnm_checkpam(&pam, PM_CHECK_BASIC, &checkRetval);
     if (allimages) {
@@ -171,7 +173,11 @@ doOneImage(const char * const name,
         
         pnm_freepamrow(tuplerow);
         
-        pnm_nextimage(fileP, eofP);
+        {
+            int eof;
+            pnm_nextimage(fileP, &eof);
+            *eofP = eof;
+        }
     }
 }
 
@@ -217,11 +223,11 @@ describeOneFile(const char * const name,
 
 
 int
-main(int argc, char *argv[]) {
+main(int argc, const char *argv[]) {
 
-    struct cmdlineInfo cmdline;
+    struct CmdlineInfo cmdline;
 
-    pnm_init(&argc, argv);
+    pm_proginit(&argc, argv);
 
     parseCommandLine(argc, argv, &cmdline);
 
diff --git a/analyzer/pammosaicknit.c b/analyzer/pammosaicknit.c
new file mode 100644
index 00000000..f0e4c7f9
--- /dev/null
+++ b/analyzer/pammosaicknit.c
@@ -0,0 +1,253 @@
+/* ----------------------------------------------------------------------
+ *
+ * Validate a mosaic knitting pattern
+ *
+ * By Scott Pakin <scott+pbm@pakin.org>
+ *
+ * ----------------------------------------------------------------------
+ *
+ * Copyright (C) 2010 Scott Pakin <scott+pbm@pakin.org>
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or (at
+ * your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful, but
+ * WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+ * General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program.  If not, see http://www.gnu.org/licenses/.
+ *
+ * ----------------------------------------------------------------------
+ */
+
+#include <stdlib.h>
+#include <stdio.h>
+#include <string.h>
+#include "mallocvar.h"
+#include "pam.h"
+
+int const max_skips = 3;      /* Maximum number of consecutive skips */
+
+/* Each pixel can be either black or white and can be valid or invalid. */
+typedef enum {
+    MK_VALID_BLACK,
+    MK_VALID_WHITE,
+    MK_INVALID_BLACK,
+    MK_INVALID_WHITE
+} mkPixelType;
+
+
+static void
+initializeInvalidColors(struct pam * const pamP,
+                        tuple *      const blackColorP,
+                        tuple *      const whiteColorP) {
+    tuple invalidBlack;
+    tuple invalidWhite;
+
+    invalidBlack = pnm_allocpamtuple(pamP);
+    invalidBlack[PAM_RED_PLANE] = pamP->maxval/3;
+    invalidBlack[PAM_GRN_PLANE] = 0;
+    invalidBlack[PAM_BLU_PLANE] = 0;
+    *blackColorP = invalidBlack;
+
+    invalidWhite = pnm_allocpamtuple(pamP);
+    invalidWhite[PAM_RED_PLANE] = pamP->maxval;
+    invalidWhite[PAM_GRN_PLANE] = pamP->maxval*3/4;
+    invalidWhite[PAM_BLU_PLANE] = pamP->maxval*3/4;
+    *whiteColorP = invalidWhite;
+}
+
+
+
+static void
+pamRowToKnitRow(struct pam *  const pamP,
+                const tuple * const pamRow,
+                mkPixelType * const knitRow) {
+
+    if (pamP->depth == 3) {
+        /* Convert from RGB. */
+        unsigned int col;
+
+        for (col = 0; col < pamP->width; ++col) {
+            double luminosity;
+
+            luminosity =
+                pamRow[col][0]*pnm_lumin_factor[0] +
+                pamRow[col][1]*pnm_lumin_factor[1] +
+                pamRow[col][2]*pnm_lumin_factor[2];
+            if (luminosity < pamP->maxval/2.0)
+                knitRow[col] = MK_VALID_BLACK;
+            else
+                knitRow[col] = MK_VALID_WHITE;
+        }
+    }
+    else {
+        /* Convert from either grayscale or black and white. */
+        unsigned int col;
+
+        for (col = 0; col < pamP->width; ++col) {
+            if (pamRow[col][0]*2 < pamP->maxval)
+                knitRow[col] = MK_VALID_BLACK;
+            else
+                knitRow[col] = MK_VALID_WHITE;
+        }
+    }
+}
+
+
+
+static void
+knitRowToPamRow(struct pam *        const pamP,
+                const mkPixelType * const knitRow,
+                tuple *             const pamRow,
+                tuple               const invalidBlack,
+                tuple               const invalidWhite) {
+
+    tuple        validBlack;
+    tuple        validWhite;
+    unsigned int col;
+
+    validBlack = pnm_allocpamtuple(pamP);
+    pnm_createBlackTuple(pamP, &validBlack);
+    validWhite = pnm_allocpamtuple(pamP);
+    validWhite[PAM_RED_PLANE] = pamP->maxval;
+    validWhite[PAM_GRN_PLANE] = pamP->maxval;
+    validWhite[PAM_BLU_PLANE] = pamP->maxval;
+
+    for (col = 0; col < pamP->width; ++col) {
+        switch (knitRow[col]) {
+          case MK_VALID_BLACK:
+              pnm_assigntuple(pamP, pamRow[col], validBlack);
+              break;
+
+          case MK_VALID_WHITE:
+              pnm_assigntuple(pamP, pamRow[col], validWhite);
+              break;
+
+          case MK_INVALID_BLACK:
+              pnm_assigntuple(pamP, pamRow[col], invalidBlack);
+              break;
+
+          case MK_INVALID_WHITE:
+              pnm_assigntuple(pamP, pamRow[col], invalidWhite);
+              break;
+
+          default:
+            abort();
+            break;
+        }
+    }
+
+    pnm_freepamtuple(validBlack);
+    pnm_freepamtuple(validWhite);
+}
+
+
+
+static void
+validateMosaicPattern(struct pam * const inPamP,
+                      struct pam * const outPamP) {
+/* --------------------------------------------------------------------------
+   Validate a mosaic knitting pattern.  A valid pattern starts with a
+   "black" row on the bottom and alternates "white" and "black" rows.
+   A "black" row can contain any arrangement of black pixels but no
+   more than max_skips consecutive white pixels.  A "white" row can
+   contain any arrangement of white pixels but no more than max_skips
+   consecutive black pixels.  Columns wrap horizontally, so that
+   characteristic is taken into consideration when tallying
+   consecutive pixels.  Invalid pixels are flagged with a red hue.  It
+   is assumed that the input image is fairly small (on the order of a
+   few tens of pixels in each dimension), so red pixels should stand
+   out fairly well when the image is zoomed in.
+--------------------------------------------------------------------------*/
+
+    tuple *       inPamRow;      /* One row of data read from the input file */
+    mkPixelType * knitRow;     /* Same as inPamRow but expressed as stitches */
+    tuple *       outPamRow;  /* One row of data to write to the output file */
+    tuple         invalidBlack; /* Color representing an invalid black input */
+    tuple         invalidWhite; /* Color representing an invalid white input */
+    mkPixelType   rowColor;     /* Base color of the current row */
+    unsigned int  row;
+
+    inPamRow = pnm_allocpamrow(inPamP);
+    outPamRow = pnm_allocpamrow(outPamP);
+    MALLOCARRAY(knitRow, inPamP->width);
+    initializeInvalidColors(outPamP, &invalidBlack, &invalidWhite);
+
+    rowColor = inPamP->height % 2 == 0 ? MK_VALID_WHITE : MK_VALID_BLACK;
+    for (row = 0; row < inPamP->height; ++row) {
+        unsigned int col;
+
+        pnm_readpamrow(inPamP, inPamRow);
+        pamRowToKnitRow(inPamP, inPamRow, knitRow);
+        for (col = 0; col < inPamP->width; ++col) {
+            if (knitRow[col] != rowColor) {
+                unsigned int runLength = 0;  /* Number of consecutive skips */
+
+                for (runLength = 0;
+                     runLength < inPamP->width &&
+                         knitRow[(col+runLength)%inPamP->width] != rowColor;
+                     ++runLength)
+                    ;
+                if (runLength > max_skips) {
+                    /* We have too many skips in a row -- mark them
+                       with the "invalid" color. */
+                    unsigned int badOffset;
+                    mkPixelType  badColor;
+
+                    badColor = rowColor == MK_VALID_WHITE ?
+                        MK_INVALID_BLACK : MK_INVALID_WHITE;
+                    for (badOffset = 0; badOffset < runLength; ++badOffset) {
+                        knitRow[(col+badOffset)%inPamP->width] = badColor;
+                    }
+                }
+                col += runLength - 1;
+            }
+        }
+        knitRowToPamRow(outPamP, knitRow, outPamRow,
+                        invalidBlack, invalidWhite);
+        pnm_writepamrow(outPamP, outPamRow);
+        rowColor = rowColor == MK_VALID_BLACK ?
+            MK_VALID_WHITE : MK_VALID_BLACK;
+    }
+
+    free(knitRow);
+    pnm_freepamrow(outPamRow);
+    pnm_freepamrow(inPamRow);
+}
+
+
+
+int
+main(int argc, const char *argv[]) {
+    struct pam   inPam;
+    struct pam   outPam;
+    const char * inputFilename;
+    FILE       * inFileP;
+
+    pm_proginit(&argc, argv);
+
+    inputFilename = (argc > 1) ? argv[1] : "-";
+    inFileP = pm_openr(inputFilename);
+
+    pnm_readpaminit(inFileP, &inPam, PAM_STRUCT_SIZE(tuple_type));
+
+    outPam = inPam;
+    outPam.file = stdout;
+    outPam.format = PAM_FORMAT;
+    outPam.depth = 3;
+    outPam.maxval = 255;
+    outPam.bytes_per_sample = 1;
+    strcpy(outPam.tuple_type, PAM_PPM_TUPLETYPE);
+    pnm_writepaminit(&outPam);
+
+    validateMosaicPattern(&inPam, &outPam);
+
+    pm_closer(inFileP);
+    return 0;
+}
+
diff --git a/analyzer/pamsharpmap.c b/analyzer/pamsharpmap.c
index ddeefe14..2f31835b 100644
--- a/analyzer/pamsharpmap.c
+++ b/analyzer/pamsharpmap.c
@@ -49,7 +49,7 @@ parseCommandLine ( int argc, char ** argv,
    was passed to us as the argv array.  We also trash *argv.
 -----------------------------------------------------------------------------*/
     optEntry *option_def = malloc(100*sizeof(optEntry));
-        /* Instructions to optParseOptions3 on how to parse our options.
+        /* Instructions to pm_optParseOptions3 on how to parse our options.
          */
     optStruct3 opt;
 
@@ -65,7 +65,7 @@ parseCommandLine ( int argc, char ** argv,
     opt.short_allowed = FALSE;  /* We have no short (old-fashioned) options */
     opt.allowNegNum = FALSE;  /* We have no parms that are negative numbers */
 
-    optParseOptions3( &argc, argv, opt, sizeof(opt), 0);
+    pm_optParseOptions3( &argc, argv, opt, sizeof(opt), 0);
         /* Uses and sets argc, argv, and some of *cmdlineP and others. */
 
     if (!contextSpec)
diff --git a/analyzer/pamsharpness.c b/analyzer/pamsharpness.c
index 8e2d9a3f..5943416f 100644
--- a/analyzer/pamsharpness.c
+++ b/analyzer/pamsharpness.c
@@ -21,8 +21,9 @@
 #include <math.h>
 
 #include "pm_c_util.h"
-#include "pam.h"
+#include "mallocvar.h"
 #include "shhopt.h"
+#include "pam.h"
 
 struct cmdlineInfo {
     /* All the information the user supplied in the command line,
@@ -35,8 +36,8 @@ struct cmdlineInfo {
 
 
 static void
-parseCommandLine ( int argc, char ** argv,
-                   struct cmdlineInfo *cmdlineP ) {
+parseCommandLine(int argc, char ** argv,
+                 struct cmdlineInfo *cmdlineP) {
 /*----------------------------------------------------------------------------
    parse program command line described in Unix standard form by argc
    and argv.  Return the information in the options as *cmdlineP.  
@@ -47,8 +48,8 @@ parseCommandLine ( int argc, char ** argv,
    Note that the strings we return are stored in the storage that
    was passed to us as the argv array.  We also trash *argv.
 -----------------------------------------------------------------------------*/
-    optEntry *option_def = malloc(100*sizeof(optEntry));
-        /* Instructions to optParseOptions3 on how to parse our options.
+    optEntry * option_def;
+        /* Instructions to pm_optParseOptions3 on how to parse our options.
          */
     optStruct3 opt;
 
@@ -56,6 +57,8 @@ parseCommandLine ( int argc, char ** argv,
 
     unsigned int contextSpec;
 
+    MALLOCARRAY_NOFAIL(option_def, 100);
+
     option_def_index = 0;   /* incremented by OPTENT3 */
     OPTENT3(0, "context",       OPT_UINT,   &cmdlineP->context,       
             &contextSpec,         0 );
@@ -64,7 +67,7 @@ parseCommandLine ( int argc, char ** argv,
     opt.short_allowed = FALSE;  /* We have no short (old-fashioned) options */
     opt.allowNegNum = FALSE;  /* We have no parms that are negative numbers */
 
-    optParseOptions3( &argc, argv, opt, sizeof(opt), 0);
+    pm_optParseOptions3( &argc, argv, opt, sizeof(opt), 0);
         /* Uses and sets argc, argv, and some of *cmdlineP and others. */
 
     if (!contextSpec)
@@ -80,6 +83,8 @@ parseCommandLine ( int argc, char ** argv,
         cmdlineP->inputFilespec = "-";
     else
         cmdlineP->inputFilespec = argv[1];
+
+    free(option_def);
 }
 
 
@@ -89,13 +94,13 @@ computeSharpness(struct pam * const inpamP,
                  tuplen **    const tuplenarray,
                  double *     const sharpnessP) {
 
-    int row;
+    unsigned int row;
     double totsharp;
     
     totsharp = 0.0;
 
     for (row = 1; row < inpamP->height-1; ++row) {
-        int col;
+        unsigned int col;
         for (col = 1; col < inpamP->width-1; ++col) {
             int dy;
             for (dy = -1; dy <= 1; ++dy) {
@@ -150,5 +155,6 @@ main(int argc, char **argv) {
 
 	pnm_freepamarrayn(tuplenarray, &inpam);
     pm_close(ifP);
+
 	return 0;
 }
diff --git a/analyzer/pamslice.c b/analyzer/pamslice.c
index ab372c8a..db0ab9c0 100644
--- a/analyzer/pamslice.c
+++ b/analyzer/pamslice.c
@@ -68,7 +68,7 @@ parseCommandLine(int argc, char ** const argv,
     opt.short_allowed = FALSE;  /* We have no short (old-fashioned) options */
     opt.allowNegNum = FALSE;  /* We have no parms that are negative numbers */
 
-    optParseOptions3(&argc, argv, opt, sizeof(opt), 0);
+    pm_optParseOptions3(&argc, argv, opt, sizeof(opt), 0);
         /* Uses and sets argc, argv, and some of *cmdlineP and others. */
 
     if (rowSpec && colSpec)
diff --git a/analyzer/pamsumm.c b/analyzer/pamsumm.c
index c9b8a7bf..c427fa7d 100644
--- a/analyzer/pamsumm.c
+++ b/analyzer/pamsumm.c
@@ -58,7 +58,7 @@ parseCommandLine(int argc, char ** const argv,
     opt.short_allowed = FALSE;  /* We have no short (old-fashioned) options */
     opt.allowNegNum = FALSE;  /* We have no parms that are negative numbers */
 
-    optParseOptions3(&argc, argv, opt, sizeof(opt), 0);
+    pm_optParseOptions3(&argc, argv, opt, sizeof(opt), 0);
         /* Uses and sets argc, argv, and some of *cmdlineP and others. */
 
     if (sumSpec + minSpec + maxSpec > 1)
@@ -73,7 +73,7 @@ parseCommandLine(int argc, char ** const argv,
     } else if (maxSpec) {
         cmdlineP->function = FN_MAX;
     } else 
-        pm_error("You must specify one of -sum, -min, or -max");
+        pm_error("You must specify one of -sum, -min, -max, or -mean");
         
     if (argc-1 > 1)
         pm_error("Too many arguments (%d).  File spec is the only argument.",
diff --git a/analyzer/pamtilt.c b/analyzer/pamtilt.c
index d70f491b..302c6247 100644
--- a/analyzer/pamtilt.c
+++ b/analyzer/pamtilt.c
@@ -47,7 +47,7 @@ abandon(void) {
 
 
 static void
-parseCommandLine(int argc, char *argv[],
+parseCommandLine(int argc, const char ** const argv,
                  struct cmdlineInfo * const cmdlineP) {
 
     static optEntry option_def[50];
@@ -76,7 +76,7 @@ parseCommandLine(int argc, char *argv[],
     opt.opt_table = option_def;
     opt.short_allowed = FALSE;          /* no short options used */
     opt.allowNegNum = FALSE;            /* don't allow negative values */
-    optParseOptions3(&argc, argv, opt, sizeof(opt), 0);
+    pm_optParseOptions3(&argc, (char **)argv, opt, sizeof(opt), 0);
 
     if (cmdlineP->hstep < 1)
         pm_error("-hstep must be at least 1 column.");
@@ -129,11 +129,13 @@ static void
 load(const struct pam * const pamP,
      unsigned int       const hstep,
      sample ***         const pixelsP,
-     unsigned int *     const hsamplesP) {
+     unsigned int *     const hsampleCtP) {
 /*----------------------------------------------------------------------------
-  read file into memory, returning array of rows
+  Read part of the raster - to wit, every hstep'th column of every row.
+  Return that sampling of the raster as *pixelsP, and the resulting
+  array width (pixels in each row) as *hsampleCtP.
 -----------------------------------------------------------------------------*/
-    unsigned int const hsamples = 1 + (pamP->width - 1) / hstep;
+    unsigned int const hsampleCt = 1 + (pamP->width - 1) / hstep;
         /* use only this many cols */
 
     unsigned int row;
@@ -143,9 +145,6 @@ load(const struct pam * const pamP,
     tuplerow = pnm_allocpamrow(pamP);
     
     MALLOCARRAY(pixels, pamP->height);
-    if (pixels == NULL)
-        pm_error("Unable to allocate array of %u pixel rows",
-                 pamP->height);
 
     if (pixels == NULL)
         pm_error("Unable to allocate array of %u rows", pamP->height);
@@ -156,19 +155,19 @@ load(const struct pam * const pamP,
 
         pnm_readpamrow(pamP, tuplerow);
 
-        MALLOCARRAY(pixels[row], hsamples);
+        MALLOCARRAY(pixels[row], hsampleCt);
         if (pixels[row] == NULL)
             pm_error("Unable to allocate %u-sample array for Row %u",
-                     hsamples, row);
+                     hsampleCt, row);
 
         /* save every hstep'th column */
-        for (i = col = 0; i < hsamples; ++i, col += hstep)
+        for (i = col = 0; i < hsampleCt; ++i, col += hstep)
             pixels[row][i] = tuplerow[col][0];
     }
 
     pnm_freepamrow(tuplerow);
 
-    *hsamplesP = hsamples;
+    *hsampleCtP = hsampleCt;
     *pixelsP   = pixels;
 }
 
@@ -192,7 +191,7 @@ static void
 replacePixelValuesWithScaledDiffs(
     const struct pam * const pamP,
     sample **          const pixels,
-    unsigned int       const hsamples,
+    unsigned int       const hsampleCt,
     unsigned int       const dstep) {
 /*----------------------------------------------------------------------------
   Replace pixel values with scaled differences used for all
@@ -204,7 +203,7 @@ replacePixelValuesWithScaledDiffs(
 
     for (row = dstep; row < pamP->height; ++row) {
         unsigned int col;
-        for (col = 0; col < hsamples; ++col) {
+        for (col = 0; col < hsampleCt; ++col) {
             int const d = pixels[row - dstep][col] - pixels[row][col];
             unsigned int const absd = d < 0 ? -d : d;
             pixels[row - dstep][col] = m * absd;  /* scale the difference */
@@ -214,9 +213,46 @@ replacePixelValuesWithScaledDiffs(
 
 
 
+static unsigned long
+totalBrightness(sample **    const pixels,
+                unsigned int const hsampleCt,
+                unsigned int const startRow,
+                float        const dy) {
+/*----------------------------------------------------------------------------
+   Total brightness of samples in the line that goes from the left edge
+   of Row 'startRow' of 'pixels' down to the right at 'dy' rows per column.
+   
+   Note that 'dy' can be negative.
+
+   Assume that whatever 'dy' is, the sloping line thus described remains
+   within 'pixels'.
+-----------------------------------------------------------------------------*/
+    unsigned long total;
+    unsigned int x;
+    float y;
+        /* The vertical location in the 'pixels' region of the currently
+           analyzed pixel.  0 is the top of the image.  Note that while
+           'pixels' is discrete pixels, this is a location in continuous
+           space.  We consider the brightness of the pixel at 1.5 to be
+           the mean of the brightness of the pixel in row 1 and the pixel
+           in row 2.
+        */
+    for (x = 0, y = startRow + 0.5, total = 0;
+         x < hsampleCt;
+         ++x, y += dy) {
+
+        assert(y >= 0.0);  /* Because of entry conditions */
+
+        total += pixels[(unsigned)y][x];
+    }
+    return total;
+}
+
+
+
 static void
 scoreAngleRegion(sample **    const pixels,
-                 unsigned int const hsamples,
+                 unsigned int const hsampleCt,
                  unsigned int const startRow,
                  unsigned int const endRow,
                  unsigned int const vstep,
@@ -234,7 +270,7 @@ scoreAngleRegion(sample **    const pixels,
    Instead of a tilt angle, we have 'dy', the slope (downward) of the lines
    in our assumed tilt.
 -----------------------------------------------------------------------------*/
-    double const tscale  = 1.0 / hsamples;
+    double const tscale  = 1.0 / hsampleCt;
 
     unsigned int row;
     double sum;
@@ -245,17 +281,10 @@ scoreAngleRegion(sample **    const pixels,
         /* Number of lines that went into 'total' */
 
     for (row = startRow, sum = 0.0, n = 0; row < endRow; row += vstep) {
-        float o;
-        long t;     /* total brightness of the samples in the line */
-        double dt;  /* mean brightness of the samples in the line */
+        double const dt =
+            tscale * totalBrightness(pixels, hsampleCt, row, dy);
+            /* mean brightness of the samples in the line */
 
-        unsigned int i;
-
-        for (i = 0, t = 0, o = 0.5;
-             i < hsamples;
-             t += pixels[(int)(row + o)][i], ++i, o += dy) {
-        }
-        dt = tscale * t;
         sum += dt * dt;
         n += 1;
     }
@@ -270,20 +299,21 @@ scoreAngle(const struct pam * const pamP,
            sample **          const pixels,
            unsigned int       const hstep,
            unsigned int       const vstep,
-           unsigned int       const hsamples,
+           unsigned int       const hsampleCt,
            float              const angle,
            float *            const scoreP) {
 /*----------------------------------------------------------------------------
-  Calculate the score for assuming the image described by *pamP and
-  'pixels' is tilted down by angle 'angle' (in degrees) from what it
-  should be.  I.e. the angle from the top edge of the paper to the
-  lines of text in the image is 'angle'.
+  Calculate the score for assuming the image described by *pamP and 'pixels'
+  is tilted down by angle 'angle' (in degrees) from what it should be.
+  I.e. the angle from the top edge of the paper to the lines of text in the
+  image is 'angle'.  Positive means the text slopes down; negative means it
+  slopes up.
 
   The score is a measure of how bright the lines of the image are if
   we assume this angle.  If the angle is right, there should be lots
   of lines that are all white, because they are between lines of text.
   If the angle is wrong, many lines will cross over lines of text and
-  thus be less that all white.  A higher score indicates 'angle' is more
+  thus be less than all white.  A higher score indicates 'angle' is more
   likely to be the angle at which the image is tilted.
 
   If 'angle' is so great that not a single line goes all the way across the
@@ -291,7 +321,7 @@ scoreAngle(const struct pam * const pamP,
   every other case, it is nonnegative.
   
   'pixels' is NOT all the pixels in the image; it is just a sampling.
-  In each row, it contains only 'hsamples' pixels, sampled from the
+  In each row, it contains only 'hsampleCt' pixels, sampled from the
   image at intervals of 'hstep' pixels.  E.g if the image is 1000
   pixels wide, pixels might be only 10 pixels wide, containing columns
   0, 100, 200, etc. of the image.
@@ -300,10 +330,10 @@ scoreAngle(const struct pam * const pamP,
 -----------------------------------------------------------------------------*/
     float  const radians = (float)angle/360 * 2 * M_PI;
     float  const dy      = hstep * tan(radians);
-        /* How much a line sinks due to the tilt when we move one sample
+        /* How much a line sinks because of the tilt when we move one sample
            ('hstep' columns of the image) to the right.
         */
-    if (fabs(dy * hsamples) > pamP->height) {
+    if (fabs(dy * hsampleCt) > pamP->height) {
         /* This is so tilted that not a single line of the image fits
            entirely on the page, so we can't do the measurement.
         */
@@ -316,17 +346,17 @@ scoreAngle(const struct pam * const pamP,
                off the page and we can't follow them.
             */
             startRow = 0;
-            endRow = pamP->height - dy * hsamples;
+            endRow = (unsigned int)(pamP->height - dy * hsampleCt + 0.5);
         } else {
             /* Lines of image rise as you go right, so the topmost lines go
                off the page and we can't follow them.
             */
-            startRow = 0 - dy * hsamples;
+            startRow = (unsigned int)(0 - dy * hsampleCt + 0.5);
             endRow = pamP->height;
         }
         assert(endRow > startRow);  /* because of 'if (fabs(dy ...' */
 
-        scoreAngleRegion(pixels, hsamples, startRow, endRow, vstep, dy,
+        scoreAngleRegion(pixels, hsampleCt, startRow, endRow, vstep, dy,
                          scoreP);
     }
 }
@@ -339,7 +369,7 @@ getBestAngleLocal(
     sample **          const pixels,
     unsigned int       const hstep,
     unsigned int       const vstep,
-    unsigned int       const hsamples,
+    unsigned int       const hsampleCt,
     float              const minangle,
     float              const maxangle,
     float              const incr,
@@ -347,7 +377,10 @@ getBestAngleLocal(
     float *            const bestAngleP,
     float *            const qualityP) {
 /*----------------------------------------------------------------------------
-  find angle of highest score within a range
+  Find the angle that gives the highest score within the range
+  [minangle, maxangle].  Those are in degrees, with positive meaning
+  the subject is rotated clockwise on the page (so a horizontal line in
+  the subject would slope down across the page).
 -----------------------------------------------------------------------------*/
     int const nsamples = ((maxangle - minangle) / incr + 1.5);
 
@@ -369,7 +402,7 @@ getBestAngleLocal(
     total = 0;
     for (i = 0; i < nsamples; i++) {
         angle = minangle + i * incr;
-        scoreAngle(pamP, pixels, hstep, vstep, hsamples, angle, &score);
+        scoreAngle(pamP, pixels, hstep, vstep, hsampleCt, angle, &score);
         results[i] = score;
         if (score > bestscore ||
             (score == bestscore && fabs(angle) < fabs(bestangle))) {
@@ -402,26 +435,36 @@ getBestAngleLocal(
 
 
 static void
-readRelevantPixels(const char *   const inputFilename,
-                   unsigned int   const hstepReq,
-                   unsigned int   const vstepReq,
-                   unsigned int * const hstepP,
-                   unsigned int * const vstepP,
-                   sample ***     const pixelsP,
-                   struct pam *   const pamP,
-                   unsigned int * const hsamplesP) {
+readSampledPixels(const char *   const inputFilename,
+                  unsigned int   const hstepReq,
+                  unsigned int   const vstepReq,
+                  unsigned int * const hstepP,
+                  unsigned int * const vstepP,
+                  sample ***     const pixelsP,
+                  struct pam *   const pamP,
+                  unsigned int * const hsampleCtP) {
 /*----------------------------------------------------------------------------
-  load the image, saving only the pixels we might actually inspect
+  Read the image.
+
+  Sample the image and return the selected pixels as *pixelsP, which is
+  an array the same height as the image, but with only certain columns
+  sampled.  Return as *hsampleCtP the number of columns sampled (i.e. the
+  width of the *pixelsP array).
+
+  Return the horizontal step size used in the sampling as *hstepP.
+  Return the appropriate vertical step size as *vstepP.
 -----------------------------------------------------------------------------*/
     FILE * ifP;
     unsigned int hstep;
     unsigned int vstep;
 
     ifP = pm_openr(inputFilename);
+
     pnm_readpaminit(ifP, pamP, PAM_STRUCT_SIZE(tuple_type));
+
     computeSteps(pamP, hstepReq, vstepReq, &hstep, &vstep);
 
-    load(pamP, hstep, pixelsP, hsamplesP);
+    load(pamP, hstep, pixelsP, hsampleCtP);
 
     *hstepP = hstep;
     *vstepP = vstep;
@@ -436,7 +479,7 @@ getAngle(const struct pam * const pamP,
          sample **          const pixels,
          unsigned int       const hstep,
          unsigned int       const vstep,
-         unsigned int       const hsamples,
+         unsigned int       const hsampleCt,
          float              const maxangle,
          float              const astep,
          float              const qmin,
@@ -448,7 +491,7 @@ getAngle(const struct pam * const pamP,
     float da;
     float lastq;        /* quality (s/n ratio) of last measurement */
     
-    getBestAngleLocal(pamP, pixels, hstep, vstep, hsamples,
+    getBestAngleLocal(pamP, pixels, hstep, vstep, hsampleCt,
                       -maxangle, maxangle, astep, verbose,
                       &a, &lastq);
 
@@ -461,14 +504,14 @@ getAngle(const struct pam * const pamP,
 
     /* make a finer search in the neighborhood */
     da = astep / 10;
-    getBestAngleLocal(pamP, pixels, hstep, vstep, hsamples,
+    getBestAngleLocal(pamP, pixels, hstep, vstep, hsampleCt,
                       a - 9 * da, a + 9 * da, da, verbose,
                       &a, &lastq);
 
     /* iterate once more unless we don't need that much accuracy */
     if (!fast) {
         da /= 10;
-        getBestAngleLocal(pamP, pixels, hstep, vstep, hsamples,
+        getBestAngleLocal(pamP, pixels, hstep, vstep, hsampleCt,
                           a - 9 * da, a + 9 * da, da, verbose,
                           &a, &lastq);
     }
@@ -478,26 +521,26 @@ getAngle(const struct pam * const pamP,
 
 
 int
-main(int argc, char *argv[]) {
+main(int argc, const char ** argv) {
 
     struct cmdlineInfo cmdline;
     struct pam pam;
     sample ** pixels;       /* pixel data */
-    unsigned int hsamples; /* horizontal samples used */
+    unsigned int hsampleCt; /* number of horizontal samples used */
     unsigned int hstep;    /* horizontal step size */
     unsigned int vstep;    /* vertical step size */
     float angle;
 
-    pgm_init(&argc, argv);              /* initialize netpbm system */
+    pm_proginit(&argc, argv);
 
     parseCommandLine(argc, argv, &cmdline);
 
-    readRelevantPixels(cmdline.inputFilename, cmdline.hstep, cmdline.vstep,
-                       &hstep, &vstep, &pixels, &pam, &hsamples);
+    readSampledPixels(cmdline.inputFilename, cmdline.hstep, cmdline.vstep,
+                      &hstep, &vstep, &pixels, &pam, &hsampleCt);
 
-    replacePixelValuesWithScaledDiffs(&pam, pixels, hsamples, cmdline.dstep);
+    replacePixelValuesWithScaledDiffs(&pam, pixels, hsampleCt, cmdline.dstep);
 
-    getAngle(&pam, pixels, hstep, vstep, hsamples,
+    getAngle(&pam, pixels, hstep, vstep, hsampleCt,
              cmdline.maxangle, cmdline.astep, cmdline.qmin,
              cmdline.fast, cmdline.verbose, &angle);
 
diff --git a/analyzer/pgmhist.c b/analyzer/pgmhist.c
index 4790ecba..1e779655 100644
--- a/analyzer/pgmhist.c
+++ b/analyzer/pgmhist.c
@@ -20,18 +20,23 @@
 
 
 
-struct cmdline_info {
+struct CmdlineInfo {
     /* 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;
+    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.
@@ -39,64 +44,120 @@ 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 */
 
-    OPTENTINIT;
+    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);
+    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 */
     opt.allowNegNum = FALSE;  /* We have no parms that are negative numbers */
 
-    optParseOptions3(&argc, (char **)argv, opt, sizeof(opt), 0);
+    pm_optParseOptions3(&argc, (char **)argv, opt, sizeof(opt), 0);
         /* Uses and sets argc, argv, and some of *cmdlineP and others. */
 
-    if (argc-1 == 0) 
+    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 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: */
@@ -107,28 +168,108 @@ buildHistogram(FILE *          const ifP,
     }
     pgm_freerow(grayrow);
 
-    *histP   = hist;
-    *maxvalP = maxval;
+    *histP = hist;
 }
 
 
 
 static void
-countCumulative(unsigned int    const hist[],
-                gray            const maxval,
-                unsigned int ** const rcountP) {
+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).
 
-    unsigned int * rcount;
-    unsigned int cumCount;
+   '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.
+
+   quantile[] must be allocated at least to size 'n'.
+
+   'n' must not be more than 100.
+-----------------------------------------------------------------------------*/
+    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.
+       */
+        assert(quantCt <= totalCt);
+
+        /* at sampleVal == mmaxval, cumCt == totalCt, so because
+           quantCt <= 'totalCt', 'sampleVal' cannot go above mmaxval.
+        */
+
+        while (cumCt < quantCt) {
+            ++sampleVal;
+            cumCt += hist[sampleVal];
+        }
+
+        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
+           of pixels that have sample value 'sampleVal' or less.
+        */
+        quantile[quantSeq-1] = sampleVal;
+    }
+}
+
+
+
+static void
+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]
+   is the number of pixels of sample value N or higher).
+
+   *rcountP is newly malloced memory.
+-----------------------------------------------------------------------------*/
+    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;
@@ -140,26 +281,215 @@ countCumulative(unsigned int    const hist[],
 
 
 static void
-report(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 int count;
+    unsigned long int const totalPixelCt = rcount[0];
+
+    unsigned int cumCount;
     unsigned int i;
 
-    printf("value  count  b%%      w%%   \n");
+    printf("value  count  b%%     w%%   \n");
     printf("-----  -----  ------  ------\n");
 
-    count = 0;
+    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);
+        }
+    }
+}
 
-    for (i = 0; i <= maxval; ++i) {
+
+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) {
-            count += hist[i];
+            cumCount += hist[i];
             printf(
-                "%5d  %5d  %5.3g%%  %5.3g%%\n", i, hist[i],
-                (float) count * 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);
+        }
+    }
+    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);
+            }
+        }
+    }
+}
+
+
+
+static void
+reportHistMachineFriendly(unsigned long int const hist[],
+                          gray              const maxval) {
+
+    unsigned int i;
+
+    for (i = 0; i <= maxval; ++i) {
+        printf("%u %lu\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]);
+}
+
+
+
+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);
+    }
+}
+
+
+
+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.
+
+   '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').
+
+   'cmdline' tells what kind of reporting to do.
+-----------------------------------------------------------------------------*/
+
+    if (cmdline.median) {
+        gray median[2];
+        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, totalPixelCt, mmaxval, quartile);
+        if (cmdline.machine)
+            reportQuantilesMachineFriendly(quartile, 4);
+        else
+            reportQuartilesHumanFriendly(quartile);
+    } else if (cmdline.decile) {
+        gray decile[10];
+        findQuantiles(10, hist, totalPixelCt, mmaxval, decile);
+        if (cmdline.machine)
+            reportQuantilesMachineFriendly(decile, 10);
+        else
+            reportDecilesHumanFriendly(decile);
+    } else {
+        if (cmdline.machine)
+            reportHistMachineFriendly(hist, mmaxval);
+        else {
+            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);
         }
     }
 }
@@ -169,11 +499,19 @@ report(unsigned int const hist[],
 int
 main(int argc, const char ** argv) {
 
-    struct cmdline_info cmdline;
+    struct CmdlineInfo cmdline;
     FILE * ifP;
+    int rows, cols;
+    int format;
     gray maxval;
-    unsigned int * rcount; /* malloc'ed array */
-    unsigned int * hist;   /* malloc'ed array */
+        /* 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);
 
@@ -181,18 +519,23 @@ main(int argc, const char ** argv) {
 
     ifP = pm_openr(cmdline.inputFileName);
 
-    buildHistogram(ifP, &hist, &maxval);
+    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;
 
-    countCumulative(hist, maxval, &rcount);
+    mmaxval = cmdline.forensic ? universalMaxval(maxval, format) : maxval;
 
-    report(hist, rcount, maxval);
+    buildHistogram(ifP, format, cols, rows, mmaxval, &hist);
+
+    reportFromHistogram(hist, mmaxval, maxval, totalPixelCt, cmdline);
 
-    free(rcount);
     free(hist);
     pm_close(ifP);
 
     return 0;
 }
-
-
-
diff --git a/analyzer/pgmtexture.c b/analyzer/pgmtexture.c
index e9a03af5..07317336 100644
--- a/analyzer/pgmtexture.c
+++ b/analyzer/pgmtexture.c
@@ -1,4 +1,4 @@
-/* pgmtexture.c - calculate textural features on a portable graymap
+/* pgmtexture.c - calculate textural features of a PGM image
 **
 ** Author: James Darrell McCauley
 **         Texas Agricultural Experiment Station
@@ -41,13 +41,19 @@
 **
 ** 05 Oct 05 - Marc Breithecker <Marc.Breithecker@informatik.uni-erlangen.de>
 **             Fix calculation or normalizing constants for d > 1.
+** 9 Jul 11  - Francois P. S. Luus <fpsluus@gmail.com> supplied fix for sum
+**             variance calculation (use F6:savg instead of F8:sentropy in
+**             F7:svar equation).
+
+
 */
 
+#include <assert.h>
 #include <math.h>
 
 #include "pm_c_util.h"
-#include "pgm.h"
 #include "mallocvar.h"
+#include "pgm.h"
 
 #define RADIX 2.0
 #define EPSILON 0.000000001
@@ -67,32 +73,60 @@
 #define F13 "Meas of Correlation-2 "
 #define F14 "Max Correlation Coeff "
 
-#define SIGN(x,y) ((y)<0 ? -fabs(x) : fabs(x))
-#define DOT fprintf(stderr,".")
-#define SWAP(a,b) {y=(a);(a)=(b);(b)=y;}
+#define SWAP(a,b) do {float const y=(a);(a)=(b);(b)=y;} while (0)
+
+
+
+static float
+sign(float const x,
+     float const y) {
+
+    return y < 0 ? -fabs(x) : fabs(x);
+}
+
+
+
+static bool const sortit = FALSE;
 
 
-static bool sortit = FALSE;
 
 static float *
-vector (int nl, int nh)
-{
-    float *v;
+vector(unsigned int const nl,
+       unsigned int const nh) {
+
+    float * v;
+
+    assert(nh >= nl);
 
     MALLOCARRAY(v, (unsigned) (nh - nl + 1));
+
     if (v == NULL)
         pm_error("Unable to allocate memory for a vector.");
+
     return v - nl;
 }
 
 
+
 static float **
-matrix (int nrl, int nrh, int ncl, int nch)
+matrix (unsigned int const nrl,
+        unsigned int const nrh,
+        unsigned int const ncl,
+        unsigned int const nch) {
+/*----------------------------------------------------------------------------
+  Allocate a float matrix with range [nrl..nrh][ncl..nch]
+
+  We do some seedy C here, subtracting an arbitrary integer from a pointer and
+  calling the result a pointer.  It normally works because the only way we'll
+  use that pointer is by adding that same integer or something greater to it.
 
-/* Allocates a float matrix with range [nrl..nrh][ncl..nch] */
-{
-    int i;
-    float **m;
+  The point of this is not to allocate memory for matrix elements that will
+  never be referenced (row < nrl or column < ncl).
+-----------------------------------------------------------------------------*/
+    unsigned int i;
+    float ** m;
+
+    assert(nrh >= nrl);
 
     /* allocate pointers to rows */
     MALLOCARRAY(m, (unsigned) (nrh - nrl + 1));
@@ -101,92 +135,102 @@ matrix (int nrl, int nrh, int ncl, int nch)
 
     m -= ncl;
 
+    assert (nch >= ncl);
+
     /* allocate rows and set pointers to them */
-    for (i = nrl; i <= nrh; i++)
-    {
+    for (i = nrl; i <= nrh; ++i) {
         MALLOCARRAY(m[i], (unsigned) (nch - ncl + 1));
         if (m[i] == NULL)
             pm_error("Unable to allocate memory for a matrix row.");
         m[i] -= ncl;
     }
+
     /* return pointer to array of pointers to rows */
     return m;
 }
 
+
+
 static void 
-results (const char * const c, const float * const a)
-{
-    int i;
+results (const char *  const name,
+         const float * const a) {
+
+    unsigned int i;
+
+    fprintf(stdout, "%s", name);
 
-    DOT;
-    fprintf (stdout, "%s", c);
     for (i = 0; i < 4; ++i)
-        fprintf (stdout, "% 1.3e ", a[i]);
-    fprintf (stdout, "% 1.3e\n", (a[0] + a[1] + a[2] + a[3]) / 4);
+        fprintf(stdout, "% 1.3e ", a[i]);
+
+    fprintf(stdout, "% 1.3e\n", (a[0] + a[1] + a[2] + a[3]) / 4);
 }
 
+
+
 static void 
-simplesrt (int n, float arr[])
-{
-    int i, j;
+simplesrt (unsigned int  const n,
+           float *       const arr) {
+
+    unsigned int j;
     float a;
 
-    for (j = 2; j <= n; j++)
-    {
+    for (j = 2; j <= n; ++j) {
+        unsigned int i;
+
         a = arr[j];
-        i = j - 1;
-        while (i > 0 && arr[i] > a)
-        {
-            arr[i + 1] = arr[i];
-            i--;
+        i = j;
+
+        while (i > 1 && arr[i-1] > a) {
+            arr[i] = arr[i-1];
+            --i;
         }
-        arr[i + 1] = a;
+        arr[i] = a;
     }
 }
 
+
+
 static void 
-mkbalanced (float **a, int n)
-{
-    int last, j, i;
-    float s, r, g, f, c, sqrdx;
+mkbalanced (float **     const a,
+            unsigned int const n) {
+
+    float const sqrdx = SQR(RADIX);
+
+    unsigned int last, i;
+    float s, r, g, f, c;
 
-    sqrdx = RADIX * RADIX;
     last = 0;
-    while (last == 0)
-    {
+    while (last == 0) {
         last = 1;
-        for (i = 1; i <= n; i++)
-        {
+        for (i = 1; i <= n; ++i) {
+            unsigned int j;
             r = c = 0.0;
-            for (j = 1; j <= n; j++)
-                if (j != i)
-                {
+            for (j = 1; j <= n; ++j) {
+                if (j != i) {
                     c += fabs (a[j][i]);
                     r += fabs (a[i][j]);
                 }
-            if (c && r)
-            {
+            }
+            if (c && r) {
                 g = r / RADIX;
                 f = 1.0;
                 s = c + r;
-                while (c < g)
-                {
+                while (c < g) {
                     f *= RADIX;
                     c *= sqrdx;
                 }
                 g = r * RADIX;
-                while (c > g)
-                {
+                while (c > g) {
                     f /= RADIX;
                     c /= sqrdx;
                 }
-                if ((c + r) / f < 0.95 * s)
-                {
+                if ((c + r) / f < 0.95 * s) {
+                    unsigned int j;
                     last = 0;
                     g = 1.0 / f;
-                    for (j = 1; j <= n; j++)
+                    for (j = 1; j <= n; ++j)
                         a[i][j] *= g;
-                    for (j = 1; j <= n; j++)
+                    for (j = 1; j <= n; ++j)
                         a[j][i] *= f;
                 }
             }
@@ -195,43 +239,43 @@ mkbalanced (float **a, int n)
 }
 
 
+
 static void 
-reduction (float **a, int n)
-{
-    int m, j, i;
-    float y, x;
+reduction (float **     const a,
+           unsigned int const n) {
 
-    for (m = 2; m < n; m++)
-    {
+    unsigned int m;
+
+    for (m = 2; m < n; ++m) {
+        unsigned int j;
+        unsigned int i;
+        float x;
         x = 0.0;
         i = m;
-        for (j = m; j <= n; j++)
-        {
-            if (fabs (a[j][m - 1]) > fabs (x))
-            {
+        for (j = m; j <= n; ++j) {
+            if (fabs(a[j][m - 1]) > fabs(x)) {
                 x = a[j][m - 1];
                 i = j;
             }
         }
-        if (i != m)
-        {
-            for (j = m - 1; j <= n; j++)
-                SWAP (a[i][j], a[m][j])  
-                    for (j = 1; j <= n; j++)
-                        SWAP (a[j][i], a[j][m]) 
-                            a[j][i] = a[j][i];
+        if (i != m) {
+            for (j = m - 1; j <= n; ++j)
+                SWAP(a[i][j], a[m][j]);
+            for (j = 1; j <= n; j++)
+                SWAP(a[j][i], a[j][m]); 
+            a[j][i] = a[j][i];
         }
-        if (x)
-        {
-            for (i = m + 1; i <= n; i++)
-            {
-                if ((y = a[i][m - 1]))
-                {
+        if (x != 0.0) {
+            unsigned int i;
+            for (i = m + 1; i <= n; ++i) {
+                float y;
+                y = a[i][m - 1];
+                if (y) {
                     y /= x;
                     a[i][m - 1] = y;
-                    for (j = m; j <= n; j++)
+                    for (j = m; j <= n; ++j)
                         a[i][j] -= y * a[m][j];
-                    for (j = 1; j <= n; j++)
+                    for (j = 1; j <= n; ++j)
                         a[j][m] += y * a[j][i];
                 }
             }
@@ -241,128 +285,140 @@ reduction (float **a, int n)
 
 
 
+static float
+norm(float **     const a,
+     unsigned int const n) {
+
+    float anorm;
+    unsigned int i;
+
+    for (i = 2, anorm = fabs(a[1][1]); i <= n; ++i) {
+        unsigned int j;
+        for (j = (i - 1); j <= n; ++j)
+            anorm += fabs(a[i][j]);
+    }
+    return anorm;
+}
+
+
+
 static void 
-hessenberg (float **a, int n, float wr[], float wi[])
-
-{
-    int nn, m, l, k, j, its, i, mmin;
-    float z, y, x, w, v, u, t, s, r, q, p, anorm;
-
-    anorm = fabs (a[1][1]);
-    for (i = 2; i <= n; i++)
-        for (j = (i - 1); j <= n; j++)
-            anorm += fabs (a[i][j]);
-    nn = n;
-    t = 0.0;
-    while (nn >= 1)
-    {
+hessenberg(float **     const a,
+           unsigned int const n,
+           float *      const wr,
+           float *      const wi) {
+
+    float const anorm = norm(a, n);
+
+    int nn;
+    float t;
+
+    assert(n >= 1);
+
+    for (nn = n, t = 0.0; nn >= 1; ) {
+        unsigned int its;
+        int l;
         its = 0;
-        do
-        {
-            for (l = nn; l >= 2; l--)
-            {
+        do {
+            float x;
+            for (l = nn; l >= 2; --l) {
+                float s;
                 s = fabs (a[l - 1][l - 1]) + fabs (a[l][l]);
                 if (s == 0.0)
                     s = anorm;
                 if ((float) (fabs (a[l][l - 1]) + s) == s)
                     break;
             }
+            assert(nn >= 1);
             x = a[nn][nn];
-            if (l == nn)
-            {
+            if (l == nn) {
                 wr[nn] = x + t;
                 wi[nn--] = 0.0;
-            }
-            else
-            {
-                y = a[nn - 1][nn - 1];
-                w = a[nn][nn - 1] * a[nn - 1][nn];
-                if (l == (nn - 1))
-                {
-                    p = 0.5 * (y - x);
-                    q = p * p + w;
-                    z = sqrt (fabs (q));
+            } else {
+                float w, y;
+                y = a[nn - 1][nn - 1];  /* initial value */
+                w = a[nn][nn - 1] * a[nn - 1][nn];  /* initial value */
+                if (l == (nn - 1)) {
+                    float const p = 0.5 * (y - x);
+                    float const q = p * p + w;
+                    float const z = sqrt(fabs(q));
                     x += t;
-                    if (q >= 0.0)
-                    {
-                        z = p + SIGN (z, p); 
-                        wr[nn - 1] = wr[nn] = x + z;
-                        if (z)
-                            wr[nn] = x - w / z;
+                    if (q >= 0.0) {
+                        float const z2 = p + sign(z, p); 
+                        wr[nn - 1] = wr[nn] = x + z2;
+                        if (z2)
+                            wr[nn] = x - w / z2;
                         wi[nn - 1] = wi[nn] = 0.0;
-                    }
-                    else
-                    {
+                    } else {
                         wr[nn - 1] = wr[nn] = x + p;
                         wi[nn - 1] = -(wi[nn] = z);
                     }
                     nn -= 2;
-                }
-                else
-                {
+                } else {
+                    int i, k, m;
+                    float p, q, r;
                     if (its == 30)
                         pm_error("Too many iterations to required "
                                  "to find %s.  Giving up", F14);
-                    if (its == 10 || its == 20)
-                    {
+                    if (its == 10 || its == 20) {
+                        int i;
+                        float s;
                         t += x;
-                        for (i = 1; i <= nn; i++)
+                        for (i = 1; i <= nn; ++i)
                             a[i][i] -= x;
-                        s = fabs (a[nn][nn - 1]) + fabs (a[nn - 1][nn - 2]);
+                        s = fabs(a[nn][nn - 1]) + fabs(a[nn - 1][nn - 2]);
                         y = x = 0.75 * s;
                         w = -0.4375 * s * s;
                     }
                     ++its;
-                    for (m = (nn - 2); m >= l; m--)
-                    {
-                        z = a[m][m];
+                    for (m = (nn - 2); m >= l; --m) {
+                        float const z = a[m][m];
+                        float s, u, v;
                         r = x - z;
                         s = y - z;
                         p = (r * s - w) / a[m + 1][m] + a[m][m + 1];
                         q = a[m + 1][m + 1] - z - r - s;
                         r = a[m + 2][m + 1];
-                        s = fabs (p) + fabs (q) + fabs (r);
+                        s = fabs(p) + fabs(q) + fabs(r);
                         p /= s;
                         q /= s;
                         r /= s;
                         if (m == l)
                             break;
-                        u = fabs (a[m][m - 1]) * (fabs (q) + fabs (r));
-                        v = fabs (p) * (fabs (a[m - 1][m - 1]) + fabs (z) + 
-                                        fabs (a[m + 1][m + 1]));
-                        if ((float) (u + v) == v)
+                        u = fabs(a[m][m - 1]) * (fabs(q) + fabs(r));
+                        v = fabs(p) * (fabs(a[m - 1][m - 1]) + fabs(z) + 
+                                       fabs(a[m + 1][m + 1]));
+                        if (u + v == v)
                             break;
                     }
-                    for (i = m + 2; i <= nn; i++)
-                    {
+                    for (i = m + 2; i <= nn; ++i) {
                         a[i][i - 2] = 0.0;
                         if (i != (m + 2))
                             a[i][i - 3] = 0.0;
                     }
-                    for (k = m; k <= nn - 1; k++)
-                    {
-                        if (k != m)
-                        {
+                    for (k = m; k <= nn - 1; ++k) {
+                        float s;
+                        if (k != m) {
                             p = a[k][k - 1];
                             q = a[k + 1][k - 1];
                             r = 0.0;
                             if (k != (nn - 1))
                                 r = a[k + 2][k - 1];
-                            if ((x = fabs (p) + fabs (q) + fabs (r)))
-                            {
+                            if ((x = fabs(p) + fabs(q) + fabs(r))) {
                                 p /= x;
                                 q /= x;
                                 r /= x;
                             }
                         }
-                        if ((s = SIGN (sqrt (p * p + q * q + r * r), p))) 
-                        {
-                            if (k == m)
-                            {
+                        s = sign(sqrt(SQR(p) + SQR(q) + SQR(r)), p);
+                        if (s) {
+                            int const mmin = nn < k + 3 ? nn : k + 3;
+                            float z;
+                            int j;
+                            if (k == m) {
                                 if (l != m)
                                     a[k][k - 1] = -a[k][k - 1];
-                            }
-                            else
+                            } else
                                 a[k][k - 1] = -s * x;
                             p += s;
                             x = p / s;
@@ -370,23 +426,18 @@ hessenberg (float **a, int n, float wr[], float wi[])
                             z = r / s;
                             q /= p;
                             r /= p;
-                            for (j = k; j <= nn; j++)
-                            {
+                            for (j = k; j <= nn; ++j) {
                                 p = a[k][j] + q * a[k + 1][j];
-                                if (k != (nn - 1))
-                                {
+                                if (k != (nn - 1)) {
                                     p += r * a[k + 2][j];
                                     a[k + 2][j] -= p * z;
                                 }
                                 a[k + 1][j] -= p * y;
                                 a[k][j] -= p * x;
                             }
-                            mmin = nn < k + 3 ? nn : k + 3;
-                            for (i = l; i <= mmin; i++)
-                            {
+                            for (i = l; i <= mmin; ++i) {
                                 p = x * a[i][k] + y * a[i][k + 1];
-                                if (k != (nn - 1))
-                                {
+                                if (k != (nn - 1)) {
                                     p += z * a[i][k + 2];
                                     a[i][k + 2] -= p * r;
                                 }
@@ -404,428 +455,497 @@ hessenberg (float **a, int n, float wr[], float wi[])
 
 
 static float 
-f1_asm (float **P, int Ng)
-
-/* Angular Second Moment */
-{
-    int i, j;
-    float sum = 0;
-
-    for (i = 0; i < Ng; ++i)
-        for (j = 0; j < Ng; ++j)
-            sum += P[i][j] * P[i][j];
-
+f1_a2m(float **     const p,
+       unsigned int const ng) {
+/*----------------------------------------------------------------------------
+  Angular Second Moment
+
+  The angular second-moment feature (ASM) f1 is a measure of homogeneity of
+  the image. In a homogeneous image, there are very few dominant gray-tone
+  transitions. Hence the P matrix for such an image will have fewer entries of
+  large magnitude.
+-----------------------------------------------------------------------------*/
+    unsigned int i;
+    float sum;
+
+    for (i = 0, sum = 0.0; i < ng; ++i) {
+        unsigned int j;
+        for (j = 0; j < ng; ++j)
+            sum += p[i][j] * p[i][j];
+    }
     return sum;
-
-    /*
-     * The angular second-moment feature (ASM) f1 is a measure of homogeneity
-     * of the image. In a homogeneous image, there are very few dominant
-     * gray-tone transitions. Hence the P matrix for such an image will have
-     * fewer entries of large magnitude.
-     */
 }
 
 
-static float 
-f2_contrast (float **P, int Ng)
-
-/* Contrast */
-{
-    int i, j, n;
-    float sum = 0, bigsum = 0;
 
-    for (n = 0; n < Ng; ++n)
-    {
-        for (i = 0; i < Ng; ++i)
-            for (j = 0; j < Ng; ++j)
+static float 
+f2_contrast(float **     const p,
+            unsigned int const ng) {
+/*----------------------------------------------------------------------------
+   Contrast
+   
+   The contrast feature is a difference moment of the P matrix and is a
+   measure of the contrast or the amount of local variations present in an
+   image.
+-----------------------------------------------------------------------------*/
+    unsigned int n;
+    float bigsum;
+
+    for (n = 0, bigsum = 0.0; n < ng; ++n) {
+        unsigned int i;
+        float sum;
+        for (i = 0, sum = 0.0; i < ng; ++i) {
+            unsigned int j;
+            for (j = 0; j < ng; ++j) {
                 if ((i - j) == n || (j - i) == n)
-                    sum += P[i][j];
-        bigsum += n * n * sum;
-
-        sum = 0;
+                    sum += p[i][j];
+            }
+        }
+        bigsum += SQR(n) * sum;
     }
     return bigsum;
-
-    /*
-     * The contrast feature is a difference moment of the P matrix and is a
-     * measure of the contrast or the amount of local variations present in an
-     * image.
-     */
 }
 
-static float 
-f3_corr (float **P, int Ng)
 
-/* Correlation */
-{
-    int i, j;
-    float sum_sqrx = 0, sum_sqry = 0, tmp, *px;
-    float meanx =0 , meany = 0 , stddevx, stddevy;
 
-    px = vector (0, Ng);
-    for (i = 0; i < Ng; ++i)
+static float 
+f3_corr(float **     const p,
+        unsigned int const ng) {
+/*----------------------------------------------------------------------------
+   Correlation
+
+   This correlation feature is a measure of gray-tone linear-dependencies in
+   the image.
+-----------------------------------------------------------------------------*/
+    unsigned int i;
+    float sumSqrx, sumSqry, tmp;
+    float * px;
+    float meanx, meany, stddevx, stddevy;
+
+    sumSqrx = 0.0; sumSqry = 0.0;
+    meanx = 0.0; meany = 0.0;
+
+    px = vector(0, ng);
+    for (i = 0; i < ng; ++i)
         px[i] = 0;
 
-    /*
-     * px[i] is the (i-1)th entry in the marginal probability matrix obtained
-     * by summing the rows of p[i][j]
-     */
-    for (i = 0; i < Ng; ++i)
-        for (j = 0; j < Ng; ++j)
-            px[i] += P[i][j];
-
+    /* px[i] is the (i-1)th entry in the marginal probability matrix obtained
+       by summing the rows of p[i][j]
+    */
+    for (i = 0; i < ng; ++i) {
+        unsigned int j;
+        for (j = 0; j < ng; ++j)
+            px[i] += p[i][j];
+    }
 
     /* Now calculate the means and standard deviations of px and py */
-    /*- fix supplied by J. Michael Christensen, 21 Jun 1991 */
-    /*- further modified by James Darrell McCauley, 16 Aug 1991 
-     *     after realizing that meanx=meany and stddevx=stddevy
-     */
-    for (i = 0; i < Ng; ++i)
-    {
-        meanx += px[i]*i;
-        sum_sqrx += px[i]*i*i;
+    for (i = 0; i < ng; ++i) {
+        meanx += px[i] * i;
+        sumSqrx += px[i] * SQR(i);
     }
+
     meany = meanx;
-    sum_sqry = sum_sqrx;
-    stddevx = sqrt (sum_sqrx - (meanx * meanx));
+    sumSqry = sumSqrx;
+    stddevx = sqrt(sumSqrx - (SQR(meanx)));
     stddevy = stddevx;
 
     /* Finally, the correlation ... */
-    for (tmp = 0, i = 0; i < Ng; ++i)
-        for (j = 0; j < Ng; ++j)
-            tmp += i*j*P[i][j];
-
+    for (i = 0, tmp = 0; i < ng; ++i) {
+        unsigned int j;
+        for (j = 0; j < ng; ++j)
+            tmp += i * j * p[i][j];
+    }
     return (tmp - meanx * meany) / (stddevx * stddevy);
-    /*
-     * This correlation feature is a measure of gray-tone linear-dependencies
-     * in the image.
-     */
 }
 
 
-static float 
-f4_var (float **P, int Ng)
-
-/* Sum of Squares: Variance */
-{
-    int i, j;
-    float mean = 0, var = 0;
-
-    /*- Corrected by James Darrell McCauley, 16 Aug 1991
-     *  calculates the mean intensity level instead of the mean of
-     *  cooccurrence matrix elements 
-     */
-    for (i = 0; i < Ng; ++i)
-        for (j = 0; j < Ng; ++j)
-            mean += i * P[i][j];
-
-    for (i = 0; i < Ng; ++i)
-        for (j = 0; j < Ng; ++j)
-            var += (i + 1 - mean) * (i + 1 - mean) * P[i][j];
 
+static float 
+f4_var (float **     const p,
+        unsigned int const ng) {
+/*----------------------------------------------------------------------------
+  Sum of Squares: Variance
+-----------------------------------------------------------------------------*/
+    unsigned int i;
+    float mean, var;
+
+    for (i = 0, mean = 0.0; i < ng; ++i) {
+        unsigned int j;
+        for (j = 0; j < ng; ++j)
+            mean += i * p[i][j];
+    }
+    for (i = 0, var = 0.0; i < ng; ++i) {
+        unsigned int j;
+        for (j = 0; j < ng; ++j)
+            var += (i + 1 - mean) * (i + 1 - mean) * p[i][j];
+    }
     return var;
 }
 
-static float 
-f5_idm (float **P, int Ng)
-
-/* Inverse Difference Moment */
-{
-    int i, j;
-    float idm = 0;
 
-    for (i = 0; i < Ng; ++i)
-        for (j = 0; j < Ng; ++j)
-            idm += P[i][j] / (1 + (i - j) * (i - j));
 
+static float 
+f5_idm (float **     const p,
+        unsigned int const ng) {
+/*----------------------------------------------------------------------------
+  Inverse Difference Moment
+-----------------------------------------------------------------------------*/
+    unsigned int i;
+    float idm;
+
+    for (i = 0, idm = 0.0; i < ng; ++i) {
+        unsigned int j;
+        for (j = 0; j < ng; ++j)
+            idm += p[i][j] / (1 + (i - j) * (i - j));
+    }
     return idm;
 }
 
-static float 
-Pxpy[2 * (PGM_MAXMAXVAL+1) + 1];
-
-static float 
-f6_savg (float **P, int Ng)
-
-/* Sum Average */
-{
-    int i, j;
-    float savg = 0;
 
-    for (i = 0; i <= 2 * Ng; ++i)
-        Pxpy[i] = 0;
 
-    for (i = 0; i < Ng; ++i)
-        for (j = 0; j < Ng; ++j)
-            Pxpy[i + j + 2] += P[i][j];
-    for (i = 2; i <= 2 * Ng; ++i)
-        savg += i * Pxpy[i];
+static float 
+f6_savg (float **     const p,
+         unsigned int const ng) {
+/*----------------------------------------------------------------------------
+   Sum Average
+-----------------------------------------------------------------------------*/
+    float pxpy[2 * (PGM_MAXMAXVAL+1) + 1];
+    unsigned int i;
+    float savg;
+
+    assert(2*ng < ARRAY_SIZE(pxpy));
+
+    for (i = 0; i <= 2 * ng; ++i)
+        pxpy[i] = 0.0;
+
+    for (i = 0; i < ng; ++i) {
+        unsigned int j;
+        for (j = 0; j < ng; ++j)
+            pxpy[i + j + 2] += p[i][j];
+    }
+    for (i = 2, savg = 0.0; i <= 2 * ng; ++i)
+        savg += i * pxpy[i];
 
     return savg;
 }
 
 
-static float 
-f7_svar (float **P, int Ng, float S) {
-/* Sum Variance */
-    int i, j;
-    float var = 0;
-
-    for (i = 0; i <= 2 * Ng; ++i)
-        Pxpy[i] = 0;
-
-    for (i = 0; i < Ng; ++i)
-        for (j = 0; j < Ng; ++j)
-            Pxpy[i + j + 2] += P[i][j];
 
-    for (i = 2; i <= 2 * Ng; ++i)
-        var += (i - S) * (i - S) * Pxpy[i];
+static float 
+f7_svar (float **     const p,
+         unsigned int const ng,
+         float        const s) {
+/*----------------------------------------------------------------------------
+   Sum Variance
+-----------------------------------------------------------------------------*/
+    float pxpy[2 * (PGM_MAXMAXVAL+1) + 1];
+    unsigned int i;
+    float var;
+
+    assert(2*ng < ARRAY_SIZE(pxpy));
+
+    for (i = 0; i <= 2 * ng; ++i)
+        pxpy[i] = 0;
+
+    for (i = 0; i < ng; ++i) {
+        unsigned int j;
+        for (j = 0; j < ng; ++j)
+            pxpy[i + j + 2] += p[i][j];
+    }
+    for (i = 2, var = 0.0; i <= 2 * ng; ++i)
+        var += (i - s) * (i - s) * pxpy[i];
 
     return var;
 }
 
-static float 
-f8_sentropy (float **P, int Ng)
-
-/* Sum Entropy */
-{
-    int i, j;
-    float sentropy = 0;
 
-    for (i = 0; i <= 2 * Ng; ++i)
-        Pxpy[i] = 0;
 
-    for (i = 0; i < Ng; ++i)
-        for (j = 0; j < Ng; ++j)
-            Pxpy[i + j + 2] += P[i][j];
-
-    for (i = 2; i <= 2 * Ng; ++i)
-        sentropy -= Pxpy[i] * log10 (Pxpy[i] + EPSILON);
+static float 
+f8_sentropy (float **     const p,
+             unsigned int const ng) {
+/*----------------------------------------------------------------------------
+   Sum Entropy
+-----------------------------------------------------------------------------*/
+    float pxpy[2 * (PGM_MAXMAXVAL+1) + 1];
+    unsigned int i;
+    float sentropy;
+
+    assert(2*ng < ARRAY_SIZE(pxpy));
+
+    for (i = 0; i <= 2 * ng; ++i)
+        pxpy[i] = 0;
+
+    for (i = 0; i < ng; ++i) {
+        unsigned int j;
+        for (j = 0; j < ng; ++j)
+            pxpy[i + j + 2] += p[i][j];
+    }
+    for (i = 2, sentropy = 0.0; i <= 2 * ng; ++i)
+        sentropy -= pxpy[i] * log10(pxpy[i] + EPSILON);
 
     return sentropy;
 }
 
 
-static float 
-f9_entropy (float **P, int Ng)
-
-/* Entropy */
-{
-    int i, j;
-    float entropy = 0;
-
-    for (i = 0; i < Ng; ++i)
-        for (j = 0; j < Ng; ++j)
-            entropy += P[i][j] * log10 (P[i][j] + EPSILON);
 
+static float 
+f9_entropy (float **     const p,
+            unsigned int const ng) {
+/*----------------------------------------------------------------------------
+   Entropy
+-----------------------------------------------------------------------------*/
+    unsigned int i;
+    float entropy;
+
+    for (i = 0, entropy = 0.0; i < ng; ++i) {
+        unsigned int j;
+        for (j = 0; j < ng; ++j)
+            entropy += p[i][j] * log10(p[i][j] + EPSILON);
+    }
     return -entropy;
 }
 
 
-static float 
-f10_dvar (float **P, int Ng)
-
-/* Difference Variance */
-{
-    int i, j;
-    double tmp;
-    double sum = 0, sum_sqr = 0, var = 0;
-
-    for (i = 0; i <= 2 * Ng; ++i)
-        Pxpy[i] = 0;
-
-    for (i = 0; i < Ng; ++i)
-        for (j = 0; j < Ng; ++j)
-            Pxpy[abs (i - j)] += P[i][j];
 
+static float 
+f10_dvar (float **     const p,
+          unsigned int const ng) {
+/*----------------------------------------------------------------------------
+   Difference Variance
+-----------------------------------------------------------------------------*/
+    double pxpy[PGM_MAXMAXVAL + 1];
+    unsigned int i;
+    double sqrNg;  /* Square of 'ng' */
+    double sum;
+    double sumSqr;
+    double var;
+
+    assert(ng <= ARRAY_SIZE(pxpy));
+
+    for (i = 0; i < ng; ++i)
+        pxpy[i] = 0;
+
+    for (i = 0; i < ng; ++i) {
+        unsigned int j;
+        for (j = 0; j < ng; ++j)
+            pxpy[abs(i - j)] += p[i][j];
+    }
     /* Now calculate the variance of Pxpy (Px-y) */
-    for (i = 0; i < Ng; ++i)
-    {
-        sum += Pxpy[i];
-        sum_sqr += Pxpy[i] * Pxpy[i];
+    for (i = 0, sum = 0.0, sumSqr = 0.0; i < ng; ++i) {
+        sum += pxpy[i];
+        sumSqr += SQR(pxpy[i]);
     }
-    tmp = Ng * Ng;
-    var = ((tmp * sum_sqr) - (sum * sum)) / (tmp * tmp);
+    sqrNg = SQR(ng);
+    var = (sqrNg * sumSqr - SQR(sum)) / SQR(sqrNg);
 
     return var;
 }
 
-static float 
-f11_dentropy (float **P, int Ng)
-    
-/* Difference Entropy */
-{
-    int i, j;
-    float sum = 0;
-
-    for (i = 0; i <= 2 * Ng; ++i)
-        Pxpy[i] = 0;
 
-    for (i = 0; i < Ng; ++i)
-        for (j = 0; j < Ng; ++j)
-            Pxpy[abs (i - j)] += P[i][j];
 
-    for (i = 0; i < Ng; ++i)
-        sum += Pxpy[i] * log10 (Pxpy[i] + EPSILON);
+static float 
+f11_dentropy (float **     const p,
+              unsigned int const ng) {
+/*----------------------------------------------------------------------------
+   Difference Entropy
+-----------------------------------------------------------------------------*/
+    float pxpy[2 * (PGM_MAXMAXVAL+1) + 1];
+    unsigned int i;
+    float sum;
+
+    assert(2*ng < ARRAY_SIZE(pxpy));
+
+    for (i = 0; i <= 2 * ng; ++i)
+        pxpy[i] = 0;
+
+    for (i = 0; i < ng; ++i) {
+        unsigned int j;
+        for (j = 0; j < ng; ++j)
+            pxpy[abs(i - j)] += p[i][j];
+    }
+    for (i = 0, sum = 0.0; i < ng; ++i)
+        sum += pxpy[i] * log10(pxpy[i] + EPSILON);
 
     return -sum;
 }
 
-static float 
-f12_icorr (float **P, int Ng)
 
-/* Information Measures of Correlation */
-{
-    int i, j;
-    float *px, *py;
-    float hx = 0, hy = 0, hxy = 0, hxy1 = 0, hxy2 = 0;
 
-    px = vector (0, Ng);
-    py = vector (0, Ng);
+static float 
+f12_icorr (float **     const p,
+           unsigned int const ng) {
+/*----------------------------------------------------------------------------
+  Information Measures of Correlation
+-----------------------------------------------------------------------------*/
+    unsigned int i;
+    float * px;
+    float * py;
+    float hx, hy, hxy, hxy1, hxy2;
+
+    px = vector(0, ng);
+    py = vector(0, ng);
 
     /*
      * px[i] is the (i-1)th entry in the marginal probability matrix obtained
      * by summing the rows of p[i][j]
      */
-    for (i = 0; i < Ng; ++i)
-    {
-        for (j = 0; j < Ng; ++j)
-        {
-            px[i] += P[i][j];
-            py[j] += P[i][j];
+    for (i = 0; i < ng; ++i) {
+        unsigned int j;
+        for (j = 0; j < ng; ++j) {
+            px[i] += p[i][j];
+            py[j] += p[i][j];
         }
     }
 
-    for (i = 0; i < Ng; ++i)
-        for (j = 0; j < Ng; ++j)
-        {
-            hxy1 -= P[i][j] * log10 (px[i] * py[j] + EPSILON);
-            hxy2 -= px[i] * py[j] * log10 (px[i] * py[j] + EPSILON);
-            hxy -= P[i][j] * log10 (P[i][j] + EPSILON);
-        }
+    hx = 0.0; hy = 0.0; hxy = 0.0; hxy1 = 0.0; hxy2 = 0.0;
 
+    for (i = 0; i < ng; ++i) {
+        unsigned int j;
+        for (j = 0; j < ng; ++j) {
+            hxy1 -= p[i][j] * log10(px[i] * py[j] + EPSILON);
+            hxy2 -= px[i] * py[j] * log10(px[i] * py[j] + EPSILON);
+            hxy  -= p[i][j] * log10 (p[i][j] + EPSILON);
+        }
+    }
     /* Calculate entropies of px and py - is this right? */
-    for (i = 0; i < Ng; ++i)
-    {
-        hx -= px[i] * log10 (px[i] + EPSILON);
-        hy -= py[i] * log10 (py[i] + EPSILON);
+    for (i = 0; i < ng; ++i) {
+        hx -= px[i] * log10(px[i] + EPSILON);
+        hy -= py[i] * log10(py[i] + EPSILON);
     }
-/*  fprintf(stderr,"hxy1=%f\thxy=%f\thx=%f\thy=%f\n",hxy1,hxy,hx,hy); */
-    return ((hxy - hxy1) / (hx > hy ? hx : hy));
+    return (hxy - hxy1) / (hx > hy ? hx : hy);
 }
 
-static float 
-f13_icorr (float **P, int Ng)
 
-/* Information Measures of Correlation */
-{
-    int i, j;
-    float *px, *py;
-    float hx = 0, hy = 0, hxy = 0, hxy1 = 0, hxy2 = 0;
 
-    px = vector (0, Ng);
-    py = vector (0, Ng);
+static float 
+f13_icorr (float **     const p, 
+           unsigned int const ng) {
+/*----------------------------------------------------------------------------
+  Information Measures of Correlation
+-----------------------------------------------------------------------------*/
+    unsigned int i;
+    float * px;
+    float * py;
+    float hx, hy, hxy, hxy1, hxy2;
+
+    px = vector(0, ng);
+    py = vector(0, ng);
 
     /*
      * px[i] is the (i-1)th entry in the marginal probability matrix obtained
      * by summing the rows of p[i][j]
      */
-    for (i = 0; i < Ng; ++i)
-    {
-        for (j = 0; j < Ng; ++j)
-        {
-            px[i] += P[i][j];
-            py[j] += P[i][j];
+    for (i = 0; i < ng; ++i) {
+        unsigned int j;
+        for (j = 0; j < ng; ++j) {
+            px[i] += p[i][j];
+            py[j] += p[i][j];
         }
     }
 
-    for (i = 0; i < Ng; ++i)
-        for (j = 0; j < Ng; ++j)
-        {
-            hxy1 -= P[i][j] * log10 (px[i] * py[j] + EPSILON);
-            hxy2 -= px[i] * py[j] * log10 (px[i] * py[j] + EPSILON);
-            hxy -= P[i][j] * log10 (P[i][j] + EPSILON);
-        }
+    hx = 0.0; hy = 0.0; hxy = 0.0; hxy1 = 0.0; hxy2 = 0.0;
 
+    for (i = 0; i < ng; ++i) {
+        unsigned int j;
+        for (j = 0; j < ng; ++j) {
+            hxy1 -= p[i][j] * log10(px[i] * py[j] + EPSILON);
+            hxy2 -= px[i] * py[j] * log10(px[i] * py[j] + EPSILON);
+            hxy  -= p[i][j] * log10(p[i][j] + EPSILON);
+        }
+    }
     /* Calculate entropies of px and py */
-    for (i = 0; i < Ng; ++i)
-    {
+    for (i = 0; i < ng; ++i) {
         hx -= px[i] * log10 (px[i] + EPSILON);
         hy -= py[i] * log10 (py[i] + EPSILON);
     }
-    /* fprintf(stderr,"hx=%f\thxy2=%f\n",hx,hxy2); */
-    return (sqrt (fabs (1 - exp (-2.0 * (hxy2 - hxy)))));
+    return sqrt(fabs(1 - exp (-2.0 * (hxy2 - hxy))));
 }
 
-static float 
-f14_maxcorr (float **P, int Ng)
 
-/* Returns the Maximal Correlation Coefficient */
-{
-    int i, j, k;
-    float *px, *py, **Q;
-    float *x, *iy, tmp;
 
-    px = vector (0, Ng);
-    py = vector (0, Ng);
-    Q = matrix (1, Ng + 1, 1, Ng + 1);
-    x = vector (1, Ng);
-    iy = vector (1, Ng);
+static float 
+f14_maxcorr (float **     const p,
+             unsigned int const ng) {
+/*----------------------------------------------------------------------------
+  The Maximal Correlation Coefficient
+-----------------------------------------------------------------------------*/
+    unsigned int i;
+    float *px, *py;
+    float ** q;
+    float * x;
+    float * iy;
+    float tmp;
+
+    px = vector(0, ng);
+    py = vector(0, ng);
+    q = matrix(1, ng + 1, 1, ng + 1);
+    x = vector(1, ng);
+    iy = vector(1, ng);
 
     /*
      * px[i] is the (i-1)th entry in the marginal probability matrix obtained
      * by summing the rows of p[i][j]
      */
-    for (i = 0; i < Ng; ++i)
-    {
-        for (j = 0; j < Ng; ++j)
-        {
-            px[i] += P[i][j];
-            py[j] += P[i][j];
+    for (i = 0; i < ng; ++i) {
+        unsigned int j;
+        for (j = 0; j < ng; ++j) {
+            px[i] += p[i][j];
+            py[j] += p[i][j];
         }
     }
 
-    /* Find the Q matrix */
-    for (i = 0; i < Ng; ++i)
-    {
-        for (j = 0; j < Ng; ++j)
-        {
-            Q[i + 1][j + 1] = 0;
-            for (k = 0; k < Ng; ++k)
-                Q[i + 1][j + 1] += P[i][k] * P[j][k] / px[i] / py[k];
+    /* Compute the Q matrix */
+    for (i = 0; i < ng; ++i) {
+        unsigned int j;
+        for (j = 0; j < ng; ++j) {
+            unsigned int k;
+            q[i + 1][j + 1] = 0;
+            for (k = 0; k < ng; ++k)
+                q[i + 1][j + 1] += p[i][k] * p[j][k] / px[i] / py[k];
         }
     }
 
     /* Balance the matrix */
-    mkbalanced (Q, Ng);
+    mkbalanced(q, ng);
     /* Reduction to Hessenberg Form */
-    reduction (Q, Ng);
+    reduction(q, ng);
     /* Finding eigenvalue for nonsymetric matrix using QR algorithm */
-    hessenberg (Q, Ng, x, iy);
+    hessenberg(q, ng, x, iy);
     if (sortit)
-        simplesrt(Ng,x);
-    /* Returns the sqrt of the second largest eigenvalue of Q */
-    for (i = 2, tmp = x[1]; i <= Ng; ++i)
+        simplesrt(ng, x);
+
+    /* Return the sqrt of the second largest eigenvalue of q */
+    for (i = 2, tmp = x[1]; i <= ng; ++i)
         tmp = (tmp > x[i]) ? tmp : x[i];
-    return sqrt (x[Ng - 1]);
+
+    return sqrt(x[ng - 1]);
 }
 
+
+
 int
-main (int argc, char *argv[]) {
-    FILE *ifp;
-    register gray **grays;
-    int tone[PGM_MAXMAXVAL+1], R0, R45, R90, angle, d = 1, x, y;
-    int argn, rows, cols, row, col;
-    int itone, jtone, tones;
-    float **P_matrix0, **P_matrix45, **P_matrix90, **P_matrix135;
-    float ASM[4], contrast[4], corr[4], var[4], idm[4], savg[4];
+main (int argc, const char ** argv) {
+
+    FILE * ifP;
+    gray ** grays;
+    unsigned int tone[PGM_MAXMAXVAL+1];
+    unsigned int r0, r45, r90;
+    unsigned int d;
+    unsigned int x, y;
+    unsigned int row;
+    int rows, cols;
+    int argn;
+    unsigned int itone;
+    unsigned int toneCt;
+    float ** p_matrix0, ** p_matrix45, ** p_matrix90, ** p_matrix135;
+    float a2m[4], contrast[4], corr[4], var[4], idm[4], savg[4];
     float sentropy[4], svar[4], entropy[4], dvar[4], dentropy[4];
     float icorr[4], maxcorr[4];
     gray maxval;
+    unsigned int i;
     const char * const usage = "[-d <d>] [pgmfile]";
 
-
-    pgm_init( &argc, argv );
+    pm_proginit(&argc, argv);
 
     argn = 1;
 
@@ -835,7 +955,7 @@ main (int argc, char *argv[]) {
         if ( argv[argn][1] == 'd' )
         {
             ++argn;
-            if ( argn == argc || sscanf( argv[argn], "%d", &d ) != 1 )
+            if ( argn == argc || sscanf( argv[argn], "%u", &d ) != 1 )
                 pm_usage( usage );
         }
         else
@@ -845,204 +965,210 @@ main (int argc, char *argv[]) {
 
     if ( argn < argc )
     {
-        ifp = pm_openr( argv[argn] );
+        ifP = pm_openr( argv[argn] );
         ++argn;
     }
     else
-        ifp = stdin;
+        ifP = stdin;
 
     if ( argn != argc )
         pm_usage( usage );
 
-    grays = pgm_readpgm (ifp, &cols, &rows, &maxval);
-    pm_close (ifp);
+    d = 1;
 
-    if (maxval > PGM_MAXMAXVAL) 
-        pm_error("The maxval of the image (%d) is too high.  \n"
-                 "This program's maximum is %d.", maxval, PGM_MAXMAXVAL);
+    grays = pgm_readpgm(ifP, &cols, &rows, &maxval);
+    pm_close (ifP);
 
     /* Determine the number of different gray scales (not maxval) */
-    for (row = PGM_MAXMAXVAL; row >= 0; --row)
-        tone[row] = -1;
-    for (row = rows - 1; row >= 0; --row)
+    for (i = 0; i <= PGM_MAXMAXVAL; ++i)
+        tone[i] = -1;
+    for (row = 0; row < rows; ++row) {
+        unsigned int col;
         for (col = 0; col < cols; ++col)
             tone[grays[row][col]] = grays[row][col];
-    for (row = PGM_MAXMAXVAL, tones = 0; row >= 0; --row)
-        if (tone[row] != -1)
-            tones++;
-    pm_message("(Image has %d graylevels.)", tones);
+    }
+    for (i = 0, toneCt = 0; i <= PGM_MAXMAXVAL; ++i) {
+        if (tone[i] != -1)
+            ++toneCt;
+    }
+    pm_message("(Image has %u gray levels.)", toneCt);
 
     /* Collapse array, taking out all zero values */
-    for (row = 0, itone = 0; row <= PGM_MAXMAXVAL; row++)
+    for (row = 0, itone = 0; row <= PGM_MAXMAXVAL; ++row)
         if (tone[row] != -1)
             tone[itone++] = tone[row];
     /* Now array contains only the gray levels present (in ascending order) */
 
     /* Allocate memory for gray-tone spatial dependence matrix */
-    P_matrix0 = matrix (0, tones, 0, tones);
-    P_matrix45 = matrix (0, tones, 0, tones);
-    P_matrix90 = matrix (0, tones, 0, tones);
-    P_matrix135 = matrix (0, tones, 0, tones);
-    for (row = 0; row < tones; ++row)
-        for (col = 0; col < tones; ++col)
-        {
-            P_matrix0[row][col] = P_matrix45[row][col] = 0;
-            P_matrix90[row][col] = P_matrix135[row][col] = 0;
+    p_matrix0   = matrix (0, toneCt, 0, toneCt);
+    p_matrix45  = matrix (0, toneCt, 0, toneCt);
+    p_matrix90  = matrix (0, toneCt, 0, toneCt);
+    p_matrix135 = matrix (0, toneCt, 0, toneCt);
+
+    for (row = 0; row < toneCt; ++row) {
+        unsigned int col;
+        for (col = 0; col < toneCt; ++col) {
+            p_matrix0 [row][col] = p_matrix45 [row][col] = 0;
+            p_matrix90[row][col] = p_matrix135[row][col] = 0;
         }
+    }
+    if (d > cols)
+        pm_error("Image is narrower (%u columns) "
+                 "than specified distance (%u)", cols, d);
 
     /* Find gray-tone spatial dependence matrix */
-    fprintf (stderr, "(Computing spatial dependence matrix...");
-    for (row = 0; row < rows; ++row)
-        for (col = 0; col < cols; ++col)
-            for (x = 0, angle = 0; angle <= 135; angle += 45)
-            {
+    pm_message("Computing spatial dependence matrix...");
+    for (row = 0; row < rows; ++row) {
+        unsigned int col;
+        for (col = 0; col < cols; ++col) {
+            unsigned int angle;
+            for (angle = 0, x = 0; angle <= 135; angle += 45) {
                 while (tone[x] != grays[row][col])
-                    x++;
-                if (angle == 0 && col + d < cols)
-                {
+                    ++x;
+                if (angle == 0 && col + d < cols) {
                     y = 0;
                     while (tone[y] != grays[row][col + d])
-                        y++;
-                    P_matrix0[x][y]++;
-                    P_matrix0[y][x]++;
+                        ++y;
+                    ++p_matrix0[x][y];
+                    ++p_matrix0[y][x];
                 }
-                if (angle == 90 && row + d < rows)
-                {
+                if (angle == 90 && row + d < rows) {
                     y = 0;
                     while (tone[y] != grays[row + d][col])
-                        y++;
-                    P_matrix90[x][y]++;
-                    P_matrix90[y][x]++;
+                        ++y;
+                    ++p_matrix90[x][y];
+                    ++p_matrix90[y][x];
                 }
-                if (angle == 45 && row + d < rows && col - d >= 0)
-                {
+                if (angle == 45 && row + d < rows && col >= d) {
                     y = 0;
                     while (tone[y] != grays[row + d][col - d])
-                        y++;
-                    P_matrix45[x][y]++;
-                    P_matrix45[y][x]++;
+                        ++y;
+                    ++p_matrix45[x][y];
+                    ++p_matrix45[y][x];
                 }
-                if (angle == 135 && row + d < rows && col + d < cols)
-                {
+                if (angle == 135 && row + d < rows && col + d < cols) {
                     y = 0;
                     while (tone[y] != grays[row + d][col + d])
-                        y++;
-                    P_matrix135[x][y]++;
-                    P_matrix135[y][x]++;
+                        ++y;
+                    ++p_matrix135[x][y];
+                    ++p_matrix135[y][x];
                 }
             }
+        }
+    }
     /* Gray-tone spatial dependence matrices are complete */
 
     /* Find normalizing constants */
-    R0 = 2 * rows * (cols - d);
-    R45 = 2 * (rows - d) * (cols - d);
-    R90 = 2 * (rows - d) * cols;
+    r0  = 2 * rows * (cols - d);
+    r45 = 2 * (rows - d) * (cols - d);
+    r90 = 2 * (rows - d) * cols;
 
     /* Normalize gray-tone spatial dependence matrix */
-    for (itone = 0; itone < tones; ++itone)
-        for (jtone = 0; jtone < tones; ++jtone)
-        {
-            P_matrix0[itone][jtone] /= R0;
-            P_matrix45[itone][jtone] /= R45;
-            P_matrix90[itone][jtone] /= R90;
-            P_matrix135[itone][jtone] /= R45;
+    for (itone = 0; itone < toneCt; ++itone) {
+        unsigned int jtone;
+        for (jtone = 0; jtone < toneCt; ++jtone) {
+            p_matrix0[itone][jtone]   /= r0;
+            p_matrix45[itone][jtone]  /= r45;
+            p_matrix90[itone][jtone]  /= r90;
+            p_matrix135[itone][jtone] /= r45;
         }
-
-    fprintf (stderr, " done.)\n");
-    fprintf (stderr, "(Computing textural features");
-    fprintf (stdout, "\n");
-    DOT;
-    fprintf (stdout,
-             "%s         0         45         90        135        Avg\n",
-             BL);
-
-    ASM[0] = f1_asm (P_matrix0, tones);
-    ASM[1] = f1_asm (P_matrix45, tones);
-    ASM[2] = f1_asm (P_matrix90, tones);
-    ASM[3] = f1_asm (P_matrix135, tones);
-    results (F1, ASM);
-
-    contrast[0] = f2_contrast (P_matrix0, tones);
-    contrast[1] = f2_contrast (P_matrix45, tones);
-    contrast[2] = f2_contrast (P_matrix90, tones);
-    contrast[3] = f2_contrast (P_matrix135, tones);
-    results (F2, contrast);
-
-
-    corr[0] = f3_corr (P_matrix0, tones);
-    corr[1] = f3_corr (P_matrix45, tones);
-    corr[2] = f3_corr (P_matrix90, tones);
-    corr[3] = f3_corr (P_matrix135, tones);
-    results (F3, corr);
-
-    var[0] = f4_var (P_matrix0, tones);
-    var[1] = f4_var (P_matrix45, tones);
-    var[2] = f4_var (P_matrix90, tones);
-    var[3] = f4_var (P_matrix135, tones);
-    results (F4, var);
-
-
-    idm[0] = f5_idm (P_matrix0, tones);
-    idm[1] = f5_idm (P_matrix45, tones);
-    idm[2] = f5_idm (P_matrix90, tones);
-    idm[3] = f5_idm (P_matrix135, tones);
-    results (F5, idm);
-
-    savg[0] = f6_savg (P_matrix0, tones);
-    savg[1] = f6_savg (P_matrix45, tones);
-    savg[2] = f6_savg (P_matrix90, tones);
-    savg[3] = f6_savg (P_matrix135, tones);
-    results (F6, savg);
-
-    sentropy[0] = f8_sentropy (P_matrix0, tones);
-    sentropy[1] = f8_sentropy (P_matrix45, tones);
-    sentropy[2] = f8_sentropy (P_matrix90, tones);
-    sentropy[3] = f8_sentropy (P_matrix135, tones);
-    svar[0] = f7_svar (P_matrix0, tones, savg[0]);
-    svar[1] = f7_svar (P_matrix45, tones, savg[1]);
-    svar[2] = f7_svar (P_matrix90, tones, savg[2]);
-    svar[3] = f7_svar (P_matrix135, tones, savg[3]);
-    results (F7, svar);
-    results (F8, sentropy);
-
-    entropy[0] = f9_entropy (P_matrix0, tones);
-    entropy[1] = f9_entropy (P_matrix45, tones);
-    entropy[2] = f9_entropy (P_matrix90, tones);
-    entropy[3] = f9_entropy (P_matrix135, tones);
-    results (F9, entropy);
-
-    dvar[0] = f10_dvar (P_matrix0, tones);
-    dvar[1] = f10_dvar (P_matrix45, tones);
-    dvar[2] = f10_dvar (P_matrix90, tones);
-    dvar[3] = f10_dvar (P_matrix135, tones);
-    results (F10, dvar);
-
-    dentropy[0] = f11_dentropy (P_matrix0, tones);
-    dentropy[1] = f11_dentropy (P_matrix45, tones);
-    dentropy[2] = f11_dentropy (P_matrix90, tones);
-    dentropy[3] = f11_dentropy (P_matrix135, tones);
+    }
+    pm_message(" ...done.");
+
+    pm_message("Computing textural features ...");
+
+    fprintf(stdout, "\n");
+    fprintf(stdout,
+            "%s         0         45         90        135        Avg\n",
+            BL);
+
+    a2m[0] = f1_a2m(p_matrix0,   toneCt);
+    a2m[1] = f1_a2m(p_matrix45,  toneCt);
+    a2m[2] = f1_a2m(p_matrix90,  toneCt);
+    a2m[3] = f1_a2m(p_matrix135, toneCt);
+    results(F1, a2m);
+
+    contrast[0] = f2_contrast(p_matrix0,   toneCt);
+    contrast[1] = f2_contrast(p_matrix45,  toneCt);
+    contrast[2] = f2_contrast(p_matrix90,  toneCt);
+    contrast[3] = f2_contrast(p_matrix135, toneCt);
+    results(F2, contrast);
+
+
+    corr[0] = f3_corr(p_matrix0,   toneCt);
+    corr[1] = f3_corr(p_matrix45,  toneCt);
+    corr[2] = f3_corr(p_matrix90,  toneCt);
+    corr[3] = f3_corr(p_matrix135, toneCt);
+    results(F3, corr);
+
+    var[0] = f4_var(p_matrix0,   toneCt);
+    var[1] = f4_var(p_matrix45,  toneCt);
+    var[2] = f4_var(p_matrix90,  toneCt);
+    var[3] = f4_var(p_matrix135, toneCt);
+    results(F4, var);
+
+
+    idm[0] = f5_idm(p_matrix0,   toneCt);
+    idm[1] = f5_idm(p_matrix45,  toneCt);
+    idm[2] = f5_idm(p_matrix90,  toneCt);
+    idm[3] = f5_idm(p_matrix135, toneCt);
+    results(F5, idm);
+
+    savg[0] = f6_savg(p_matrix0,  toneCt);
+    savg[1] = f6_savg(p_matrix45,  toneCt);
+    savg[2] = f6_savg(p_matrix90,  toneCt);
+    savg[3] = f6_savg(p_matrix135, toneCt);
+    results(F6, savg);
+
+    svar[0] = f7_svar(p_matrix0,   toneCt, savg[0]);
+    svar[1] = f7_svar(p_matrix45,  toneCt, savg[1]);
+    svar[2] = f7_svar(p_matrix90,  toneCt, savg[2]);
+    svar[3] = f7_svar(p_matrix135, toneCt, savg[3]);
+    results(F7, svar);
+
+    sentropy[0] = f8_sentropy(p_matrix0,   toneCt);
+    sentropy[1] = f8_sentropy(p_matrix45,  toneCt);
+    sentropy[2] = f8_sentropy(p_matrix90,  toneCt);
+    sentropy[3] = f8_sentropy(p_matrix135, toneCt);
+    results(F8, sentropy);
+
+    entropy[0] = f9_entropy(p_matrix0,   toneCt);
+    entropy[1] = f9_entropy(p_matrix45,  toneCt);
+    entropy[2] = f9_entropy(p_matrix90,  toneCt);
+    entropy[3] = f9_entropy(p_matrix135, toneCt);
+    results(F9, entropy);
+
+    dvar[0] = f10_dvar(p_matrix0,   toneCt);
+    dvar[1] = f10_dvar(p_matrix45,  toneCt);
+    dvar[2] = f10_dvar(p_matrix90,  toneCt);
+    dvar[3] = f10_dvar(p_matrix135, toneCt);
+    results(F10, dvar);
+
+    dentropy[0] = f11_dentropy(p_matrix0,   toneCt);
+    dentropy[1] = f11_dentropy(p_matrix45,  toneCt);
+    dentropy[2] = f11_dentropy(p_matrix90,  toneCt);
+    dentropy[3] = f11_dentropy(p_matrix135, toneCt);
     results (F11, dentropy);
 
-    icorr[0] = f12_icorr (P_matrix0, tones);
-    icorr[1] = f12_icorr (P_matrix45, tones);
-    icorr[2] = f12_icorr (P_matrix90, tones);
-    icorr[3] = f12_icorr (P_matrix135, tones);
-    results (F12, icorr);
-
-    icorr[0] = f13_icorr (P_matrix0, tones);
-    icorr[1] = f13_icorr (P_matrix45, tones);
-    icorr[2] = f13_icorr (P_matrix90, tones);
-    icorr[3] = f13_icorr (P_matrix135, tones);
-    results (F13, icorr);
-
-    maxcorr[0] = f14_maxcorr (P_matrix0, tones);
-    maxcorr[1] = f14_maxcorr (P_matrix45, tones);
-    maxcorr[2] = f14_maxcorr (P_matrix90, tones);
-    maxcorr[3] = f14_maxcorr (P_matrix135, tones);
-    results (F14, maxcorr);
-
-
-    fprintf (stderr, " done.)\n");
+    icorr[0] = f12_icorr(p_matrix0,   toneCt);
+    icorr[1] = f12_icorr(p_matrix45,  toneCt);
+    icorr[2] = f12_icorr(p_matrix90,  toneCt);
+    icorr[3] = f12_icorr(p_matrix135, toneCt);
+    results(F12, icorr);
+
+    icorr[0] = f13_icorr(p_matrix0,   toneCt);
+    icorr[1] = f13_icorr(p_matrix45,  toneCt);
+    icorr[2] = f13_icorr(p_matrix90,  toneCt);
+    icorr[3] = f13_icorr(p_matrix135, toneCt);
+    results(F13, icorr);
+
+    maxcorr[0] = f14_maxcorr(p_matrix0,   toneCt);
+    maxcorr[1] = f14_maxcorr(p_matrix45,  toneCt);
+    maxcorr[2] = f14_maxcorr(p_matrix90,  toneCt);
+    maxcorr[3] = f14_maxcorr(p_matrix135, toneCt);
+    results(F14, maxcorr);
+
+    pm_message(" ...done.");
 
     return 0;
 }
diff --git a/analyzer/pnmhistmap.c b/analyzer/pnmhistmap.c
index 7e504734..fc1bbbde 100644
--- a/analyzer/pnmhistmap.c
+++ b/analyzer/pnmhistmap.c
@@ -17,6 +17,7 @@
  * - Deal properly with maxvals other than 256
  */
 
+#include <assert.h>
 #include <string.h>
 
 #include "pm_c_util.h"
@@ -26,8 +27,6 @@
 
 static double const epsilon = .00001;
 
-#define SCALE_H(value) (hscale_unity ? (value) : (int)((value) * hscale))
-
 enum wantedColor {WANT_RED=0, WANT_GRN=1, WANT_BLU=2};
 
 struct cmdlineInfo {
@@ -53,14 +52,14 @@ struct cmdlineInfo {
 
 
 static void
-parseCommandLine(int argc, char ** argv,
+parseCommandLine(int argc, const char ** argv,
                  struct cmdlineInfo *cmdlineP) {
 /*----------------------------------------------------------------------------
    Note that the file spec array we return is stored in the storage that
    was passed to us as the argv array.
 -----------------------------------------------------------------------------*/
     optEntry *option_def;
-        /* Instructions to optParseOptions3 on how to parse our options.
+        /* Instructions to pm_optParseOptions3 on how to parse our options.
          */
     optStruct3 opt;
 
@@ -93,7 +92,7 @@ parseCommandLine(int argc, char ** argv,
     opt.short_allowed = FALSE;  /* We have no short (old-fashioned) options */
     opt.allowNegNum = FALSE;  /* We may have parms that are negative numbers */
 
-    optParseOptions3(&argc, argv, opt, sizeof(opt), 0);
+    pm_optParseOptions3(&argc, (char **)argv, opt, sizeof(opt), 0);
         /* Uses and sets argc, argv, and some of *cmdlineP and others. */
 
     if (!lvalSpec)
@@ -118,7 +117,7 @@ parseCommandLine(int argc, char ** argv,
         cmdlineP->inputFilespec = "-";
     else if (argc-1 != 1)
         pm_error("Program takes zero or one argument (filename).  You "
-                 "specified %d", argc-1);
+                 "specified %u", argc-1);
     else
         cmdlineP->inputFilespec = argv[1];
 }
@@ -127,7 +126,7 @@ parseCommandLine(int argc, char ** argv,
 
 static unsigned int
 maxSlotCount(const unsigned int * const hist,
-             unsigned int         const hist_width,
+             unsigned int         const histWidth,
              bool                 const no_white,
              bool                 const no_black) {
 /*----------------------------------------------------------------------------
@@ -138,7 +137,7 @@ maxSlotCount(const unsigned int * const hist,
     unsigned int i;
 
     unsigned int const start = (no_black ? 1 : 0);
-    unsigned int const finish = (no_white ? hist_width - 1 : hist_width);
+    unsigned int const finish = (no_white ? histWidth - 1 : histWidth);
     for (hmax = 0, i = start; i < finish; ++i)
         if (hmax < hist[i])
             hmax = hist[i];
@@ -150,36 +149,52 @@ maxSlotCount(const unsigned int * const hist,
 
 static void
 clipHistogram(unsigned int * const hist,
-              unsigned int   const hist_width,
+              unsigned int   const histWidth,
               unsigned int   const hmax) {
 
     unsigned int i;
 
-    for (i = 0; i < hist_width; ++i)
+    for (i = 0; i < histWidth; ++i)
         hist[i] = MIN(hmax, hist[i]);
 }
 
 
 
 static void
-pgm_hist(FILE *       const ifP,
-         int          const cols,
-         int          const rows,
-         xelval       const maxval,
-         int          const format,
-         bool         const dots,
-         bool         const no_white,
-         bool         const no_black,
-         bool         const verbose,
-         xelval       const startval,
-         xelval       const endval,
-         unsigned int const hist_width,
-         unsigned int const hist_height,
-         bool         const clipSpec,
-         unsigned int const clipCount,
-         double       const hscale) {
-
-    bool const hscale_unity = hscale - 1 < epsilon;
+countComp(xelval         const value,
+          xelval         const startval,
+          xelval         const endval,
+          unsigned int   const histWidth,
+          unsigned int * const hist) {
+
+    double const hscale = (float)(histWidth-1) / (endval - startval - 1);
+
+    if (value >= startval && value < endval) {
+        unsigned int const bin = ROUNDU((value-startval) * hscale);
+
+        assert(bin < histWidth);
+        ++hist[bin];
+    }
+}
+
+
+
+static void
+pgmHist(FILE *       const ifP,
+        int          const cols,
+        int          const rows,
+        xelval       const maxval,
+        int          const format,
+        bool         const dots,
+        bool         const no_white,
+        bool         const no_black,
+        bool         const verbose,
+        xelval       const startval,
+        xelval       const endval,
+        unsigned int const histWidth,
+        unsigned int const histHeight,
+        bool         const clipSpec,
+        unsigned int const clipCount) {
 
     gray * grayrow;
     bit ** bits;
@@ -188,15 +203,15 @@ pgm_hist(FILE *       const ifP,
     double vscale;
     unsigned int hmax;
     
-    MALLOCARRAY(ghist, hist_width);
+    MALLOCARRAY(ghist, histWidth);
     if (ghist == NULL)
-        pm_error("Not enough memory for histogram array (%d bytes)",
-                  hist_width * sizeof(int));
-    bits = pbm_allocarray(hist_width, hist_height);
+        pm_error("Not enough memory for histogram array (%u bytes)",
+                 histWidth * (unsigned)sizeof(int));
+    bits = pbm_allocarray(histWidth, histHeight);
     if (bits == NULL)
-        pm_error("no space for output array (%d bits)",
-                 hist_width * hist_height);
-    memset(ghist, 0, hist_width * sizeof(ghist[0]));
+        pm_error("no space for output array (%u bits)",
+                 histWidth * histHeight);
+    memset(ghist, 0, histWidth * sizeof(ghist[0]));
 
     /* read the pixel values into the histogram arrays */
     grayrow = pgm_allocrow(cols);
@@ -206,12 +221,8 @@ pgm_hist(FILE *       const ifP,
 
     for (i = rows; i > 0; --i) {
         pgm_readpgmrow (ifP, grayrow, cols, maxval, format);
-        for (j = cols-1; j >= 0; --j) {
-            int const value = grayrow[j];
-
-            if (value >= startval && value <= endval)
-                ++ghist[SCALE_H(value-startval)];
-        }
+        for (j = cols-1; j >= 0; --j)
+            countComp(grayrow[j], startval, endval, histWidth, ghist);
     }
     pgm_freerow(grayrow);
 
@@ -221,33 +232,35 @@ pgm_hist(FILE *       const ifP,
     if (clipSpec)
         hmax = clipCount;
     else 
-        hmax = maxSlotCount(ghist, hist_width, no_white, no_black);
+        hmax = maxSlotCount(ghist, histWidth, no_white, no_black);
+
+    assert(hmax > 0);
 
     if (verbose)
         pm_message("Done: height = %u", hmax);
 
-    clipHistogram(ghist, hist_width, hmax);
+    clipHistogram(ghist, histWidth, hmax);
 
-    vscale = (double) hist_height / hmax;
+    vscale = (double) histHeight / hmax;
 
-    for (i = 0; i < hist_width; ++i) {
-        int mark = hist_height - (int)(vscale * ghist[i]);
+    for (i = 0; i < histWidth; ++i) {
+        int mark = histHeight - (int)(vscale * ghist[i]);
         for (j = 0; j < mark; ++j)
             bits[j][i] = PBM_BLACK;
-        if (j < hist_height)
+        if (j < histHeight)
             bits[j++][i] = PBM_WHITE;
-        for ( ; j < hist_height; ++j)
+        for ( ; j < histHeight; ++j)
             bits[j][i] = dots ? PBM_BLACK : PBM_WHITE;
     }
 
-    pbm_writepbm(stdout, bits, hist_width, hist_height, 0);
+    pbm_writepbm(stdout, bits, histWidth, histHeight, 0);
 }
 
 
 
 static unsigned int
 maxSlotCountAll(unsigned int *       const hist[3],
-                unsigned int         const hist_width,
+                unsigned int         const histWidth,
                 bool                 const no_white,
                 bool                 const no_black) {
 /*----------------------------------------------------------------------------
@@ -264,7 +277,7 @@ maxSlotCountAll(unsigned int *       const hist[3],
         if (hist[color])
             hmax = MAX(hmax, 
                        maxSlotCount(hist[color], 
-                                    hist_width, no_white, no_black));
+                                    histWidth, no_white, no_black));
     
     return hmax;
 }
@@ -273,107 +286,132 @@ maxSlotCountAll(unsigned int *       const hist[3],
 
 static void
 createHist(bool             const colorWanted[3],
-           unsigned int     const hist_width,
+           unsigned int     const histWidth,
            unsigned int * (* const histP)[3]) {
 /*----------------------------------------------------------------------------
    Allocate the histogram arrays and set each slot count to zero.
 -----------------------------------------------------------------------------*/
     unsigned int color;
 
-    for (color = 0; color < 3; ++color)
+    for (color = 0; color < 3; ++color) {
         if (colorWanted[color]) {
             unsigned int * hist;
             unsigned int i;
-            MALLOCARRAY(hist, hist_width);
+            MALLOCARRAY(hist, histWidth);
             if (hist == NULL)
-                pm_error ("Not enough memory for histogram arrays (%u bytes)",
-                          hist_width * sizeof(int) * 3);
+                pm_error("Not enough memory for histogram arrays (%u bytes)",
+                         histWidth * (unsigned)sizeof(hist[0]) * 3);
 
-            for (i = 0; i < hist_width; ++i)
+            for (i = 0; i < histWidth; ++i)
                 hist[i] = 0;
             (*histP)[color] = hist;
         } else
             (*histP)[color] = NULL;
+    }
 }
 
 
 
 static void
 clipHistogramAll(unsigned int * const hist[3],
-                 unsigned int   const hist_width,
+                 unsigned int   const histWidth,
                  unsigned int   const hmax) {
 
     unsigned int color;
 
     for (color = 0; color < 3; ++color)
         if (hist[color])
-            clipHistogram(hist[color], hist_width, hmax);
+            clipHistogram(hist[color], histWidth, hmax);
 }
 
 
 
 static void
-ppm_hist(FILE *       const ifP,
-         int          const cols,
-         int          const rows,
-         xelval       const maxval,
-         int          const format,
-         bool         const dots,
-         bool         const no_white,
-         bool         const no_black,
-         bool         const colorWanted[3],
-         bool         const verbose,
-         xelval       const startval,
-         xelval       const endval,
-         unsigned int const hist_width,
-         unsigned int const hist_height,
-         bool         const clipSpec,
-         unsigned int const clipCount,
-         double       const hscale) {
-
-    bool const hscale_unity = hscale - 1 < epsilon;
-
+fillPpmBins(FILE *          const ifP,
+            unsigned int    const cols,
+            unsigned int    const rows,
+            xelval          const maxval,
+            int             const format,
+            bool            const colorWanted[3],
+            bool            const verbose,
+            xelval          const startval,
+            xelval          const endval,
+            unsigned int    const histWidth,
+            unsigned int ** const hist) {
+/*----------------------------------------------------------------------------
+   For each wanted color component, given by colorWanted[], hist[color] is the
+   histogram.  Each histogram as 'histWidth' bins; we ignore color component
+   values less than 'startval' and greater than or equal to 'endval' and
+   spread the rest evenly across the 'histWidth' bins.
+
+   We get the color component values from the PNM image on *ifP,
+   which is positioned to the raster, whose format is described
+   by 'cols', 'rows', 'maxval', and 'format'.
+-----------------------------------------------------------------------------*/
     pixel * pixrow;
-    pixel ** pixels;
-    int i, j;
-    unsigned int * hist[3];  /* Subscript is enum wantedColor */
-    double vscale;
-    unsigned int hmax;
-
-    createHist(colorWanted, hist_width, &hist);
+    unsigned int row;
 
-    if ((pixels = ppm_allocarray (hist_width, hist_height)) == NULL)
-        pm_error("no space for output array (%d pixels)",
-                 hist_width * hist_height);
-    for (i = 0; i < hist_height; ++i)
-        memset(pixels[i], 0, hist_width * sizeof(pixels[i][0]));
-
-    /* read the pixel values into the histogram arrays */
     pixrow = ppm_allocrow(cols);
 
     if (verbose)
         pm_message("making histogram...");
 
-    for (i = rows; i > 0; --i) {
+    for (row = 0; row < rows; ++row) {
+        unsigned int col;
         ppm_readppmrow(ifP, pixrow, cols, maxval, format);
-        for (j = cols-1; j >= 0; --j) {
-            int value;
-
-            if (colorWanted[WANT_RED] && 
-                (value = PPM_GETR(pixrow[j])) >= startval && 
-                value <= endval)
-                hist[WANT_RED][SCALE_H(value-startval)]++;
-            if (colorWanted[WANT_GRN] && 
-                (value = PPM_GETG(pixrow[j])) >= startval && 
-                value <= endval)
-                hist[WANT_GRN][SCALE_H(value-startval)]++;
-            if (colorWanted[WANT_BLU] && 
-                (value = PPM_GETB(pixrow[j])) >= startval && 
-                value <= endval)
-                hist[WANT_BLU][SCALE_H(value-startval)]++;
+        for (col = 0; col < cols; ++col) {
+            if (colorWanted[WANT_RED])
+                countComp(PPM_GETR(pixrow[col]),
+                          startval, endval, histWidth, hist[WANT_RED]);
+
+            if (colorWanted[WANT_GRN])
+                countComp(PPM_GETG(pixrow[col]),
+                          startval, endval, histWidth, hist[WANT_GRN]);
+
+            if (colorWanted[WANT_BLU])
+                countComp(PPM_GETB(pixrow[col]),
+                          startval, endval, histWidth, hist[WANT_BLU]);
         }
     }
     ppm_freerow(pixrow);
+}
+
+
+
+static void
+ppmHist(FILE *       const ifP,
+        unsigned int const cols,
+        unsigned int const rows,
+        xelval       const maxval,
+        int          const format,
+        bool         const dots,
+        bool         const no_white,
+        bool         const no_black,
+        bool         const colorWanted[3],
+        bool         const verbose,
+        xelval       const startval,
+        xelval       const endval,
+        unsigned int const histWidth,
+        unsigned int const histHeight,
+        bool         const clipSpec,
+        unsigned int const clipCount) {
+
+    pixel ** pixels;
+    unsigned int i;
+    unsigned int * hist[3];  /* Subscript is enum wantedColor */
+    double vscale;
+    unsigned int hmax;
+
+    createHist(colorWanted, histWidth, &hist);
+
+    if ((pixels = ppm_allocarray (histWidth, histHeight)) == NULL)
+        pm_error("no space for output array (%u pixels)",
+                 histWidth * histHeight);
+    for (i = 0; i < histHeight; ++i)
+        memset(pixels[i], 0, histWidth * sizeof(pixels[i][0]));
+
+    fillPpmBins(ifP, cols, rows, maxval, format, colorWanted, verbose,
+                startval, endval, histWidth, hist);
 
     /* find the highest-valued slot and set the vertical scale value */
     if (verbose)
@@ -381,22 +419,24 @@ ppm_hist(FILE *       const ifP,
     if (clipSpec)
         hmax = clipCount;
     else 
-        hmax = maxSlotCountAll(hist, hist_width, no_white, no_black);
+        hmax = maxSlotCountAll(hist, histWidth, no_white, no_black);
 
-    clipHistogramAll(hist, hist_width, hmax);
+    assert(hmax > 0);
 
-    vscale = (double) hist_height / hmax;
-    if (verbose)
-        pm_message("Done: height = %d, vertical scale factor = %g", 
+    clipHistogramAll(hist, histWidth, hmax);
+
+    vscale = (double) histHeight / hmax;
+    if (verbose && pm_have_float_format())
+        pm_message("Done: height = %u, vertical scale factor = %g", 
                    hmax, vscale);
 
-    for (i = 0; i < hist_width; ++i) {
+    for (i = 0; i < histWidth; ++i) {
         if (hist[WANT_RED]) {
             unsigned int j;
             bool plotted;
             plotted = FALSE;
-            for (j = hist_height - (int)(vscale * hist[WANT_RED][i]); 
-                 j < hist_height && !plotted; 
+            for (j = histHeight - (int)(vscale * hist[WANT_RED][i]); 
+                 j < histHeight && !plotted; 
                  ++j) {
                 PPM_PUTR(pixels[j][i], maxval);
                 plotted = dots;
@@ -406,8 +446,8 @@ ppm_hist(FILE *       const ifP,
             unsigned int j;
             bool plotted;
             plotted = FALSE;
-            for (j = hist_height - (int)(vscale * hist[WANT_GRN][i]); 
-                 j < hist_height && !plotted; 
+            for (j = histHeight - (int)(vscale * hist[WANT_GRN][i]); 
+                 j < histHeight && !plotted; 
                  ++j) {
                 PPM_PUTG(pixels[j][i], maxval);
                 plotted = dots;
@@ -417,33 +457,46 @@ ppm_hist(FILE *       const ifP,
             unsigned int j;
             bool plotted;
             plotted = FALSE;
-            for (j = hist_height - (int)(vscale * hist[WANT_BLU][i]); 
-                 j < hist_height && !plotted; 
+            for (j = histHeight - (int)(vscale * hist[WANT_BLU][i]); 
+                 j < histHeight && !plotted; 
                  ++j) {
                 PPM_PUTB(pixels[j][i], maxval);
                 plotted = dots;
             }
         }
     }
-    ppm_writeppm(stdout, pixels, hist_width, hist_height, maxval, 0);
+    ppm_writeppm(stdout, pixels, histWidth, histHeight, maxval, 0);
+}
+
+
+
+static void
+reportScale(unsigned int const histWidth,
+            unsigned int const range,
+            bool         const verbose) {
+
+    double const hscale = (float)(histWidth-1) / (range-1);
+
+    if (hscale - 1.0 < epsilon && verbose && pm_have_float_format())
+        pm_message("Horizontal scale factor: %g", hscale);
 }
 
 
 
 int
-main(int argc, char ** argv) {
+main(int argc, const char ** argv) {
 
     struct cmdlineInfo cmdline;
     FILE * ifP;
     int cols, rows;
     xelval maxval;
     int format;
-    unsigned int hist_width;
+    unsigned int histWidth;
     unsigned int range;
-    double hscale;
-    int hmax;
+    unsigned int hmax;
+    xelval startval, endval;
 
-    pnm_init (&argc, argv);
+    pm_proginit(&argc, argv);
 
     parseCommandLine(argc, argv, &cmdline);
 
@@ -451,34 +504,33 @@ main(int argc, char ** argv) {
 
     pnm_readpnminit(ifP, &cols, &rows, &maxval, &format);
 
-    range = MIN(maxval, cmdline.rval) - cmdline.lval + 1;
+    startval = cmdline.lval;
+    endval   = MIN(maxval, cmdline.rval) + 1;
+
+    range = endval - startval;
 
     if (cmdline.widthSpec)
-        hist_width = cmdline.width;
+        histWidth = cmdline.width;
     else
-        hist_width = range;
-
-    hscale = (float)(hist_width-1) / (range-1);
-    if (hscale - 1.0 < epsilon && cmdline.verbose)
-        pm_message("Horizontal scale factor: %g (maxval = %u)", 
-                   hscale, maxval);
+        histWidth = range;
 
+    reportScale(histWidth, range, cmdline.verbose);
     if (cmdline.nmaxSpec)
-        hmax = cols * rows / hist_width * cmdline.nmax;
+        hmax = cols * rows / histWidth * cmdline.nmax;
 
     switch (PNM_FORMAT_TYPE(format)) {
     case PPM_TYPE:
-        ppm_hist(ifP, cols, rows, maxval, format,
-                 cmdline.dots, cmdline.white, cmdline.black,
-                 cmdline.colorWanted,
-                 cmdline.verbose, cmdline.lval, cmdline.rval, 
-                 hist_width, cmdline.height, cmdline.nmaxSpec, hmax, hscale);
+        ppmHist(ifP, cols, rows, maxval, format,
+                cmdline.dots, cmdline.white, cmdline.black,
+                cmdline.colorWanted,
+                cmdline.verbose, startval, endval,
+                histWidth, cmdline.height, cmdline.nmaxSpec, hmax);
         break;
     case PGM_TYPE:
-        pgm_hist(ifP, cols, rows, maxval, format,
-                 cmdline.dots, cmdline.white, cmdline.black,
-                 cmdline.verbose, cmdline.lval, cmdline.rval,
-                 hist_width, cmdline.height, cmdline.nmaxSpec, hmax, hscale);
+        pgmHist(ifP, cols, rows, maxval, format,
+                cmdline.dots, cmdline.white, cmdline.black,
+                cmdline.verbose, startval, endval,
+                histWidth, cmdline.height, cmdline.nmaxSpec, hmax);
         break;
     case PBM_TYPE:
         pm_error("Cannot do a histogram of a a PBM file");
diff --git a/analyzer/pnmpsnr.c b/analyzer/pnmpsnr.c
index 1160fff6..af74e8c8 100644
--- a/analyzer/pnmpsnr.c
+++ b/analyzer/pnmpsnr.c
@@ -8,30 +8,91 @@
  *  Copyright (C) 1994-2000 Ullrich Hafner <hafner@bigfoot.de>
  */
 
+#include <assert.h>
 #include <string.h>
 #include <stdio.h>
 #include <math.h>
 
 #include "pm_c_util.h"
+#include "mallocvar.h"
 #include "nstring.h"
 #include "pam.h"
+#include "shhopt.h"
+
+
+
+struct CmdlineInfo {
+    /* All the information the user supplied in the command line,
+       in a form easy for the program to use.
+    */
+    const char * inputFile1Name;  /* Name of first input file */
+    const char * inputFile2Name;  /* Name of second input file */
+    unsigned int rgb;
+};
+
+
+
+static void
+parseCommandLine(int argc, const char ** argv,
+                 struct CmdlineInfo * const cmdlineP) {
+/*----------------------------------------------------------------------------
+   Note that the file spec array we return is stored in the storage that
+   was passed to as as the argv array.
+-----------------------------------------------------------------------------*/
+    optEntry * option_def;
+        /* Instructions to pm_optParseOptions3 on how to parse our options.
+         */
+    optStruct3 opt;
+
+    unsigned int option_def_index;
+
+    MALLOCARRAY_NOFAIL(option_def, 100);
+    
+    option_def_index = 0;   /* incremented by OPTENT3 */
+    OPTENT3(0,   "rgb",      OPT_FLAG,  NULL, &cmdlineP->rgb,       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 (argc-1 < 2) 
+        pm_error("Takes two arguments:  names of the two files to compare");
+    else {
+        cmdlineP->inputFile1Name = argv[1];
+        cmdlineP->inputFile2Name = argv[2];
+
+        if (argc-1 > 2)
+            pm_error("Too many arguments (%u).  The only arguments are "
+                     "the names of the two files to compare", argc-1);
+    }
+
+    free(option_def);
+}
+
 
-#define MAXFILES 16
 
 static int
-udiff(unsigned int const subtrahend, unsigned int const subtractor) {
-    return subtrahend-subtractor;
+udiff(unsigned int const subtrahend,
+      unsigned int const subtractor) {
+
+    return subtrahend - subtractor;
 }
 
 
+
 static double
 square(double const arg) {
-    return(arg*arg);
+    return(arg * arg);
 }
 
 
+
 static void
-validate_input(const struct pam pam1, const struct pam pam2) {
+validateInput(struct pam const pam1,
+              struct pam const pam2) {
 
     if (pam1.width != pam2.width)
         pm_error("images are not the same width, so can't be compared.  "
@@ -52,159 +113,325 @@ validate_input(const struct pam pam1, const struct pam pam2) {
                  "maxval of one of them.",
                  (unsigned int) pam1.maxval, (unsigned int) pam2.maxval);
 
-    if (strcmp(pam1.tuple_type, pam2.tuple_type) != 0)
+    if (!streq(pam1.tuple_type, pam2.tuple_type))
         pm_error("images are not of the same type.  The tuple types are "
                  "'%s' and '%s', respectively.",
                  pam1.tuple_type, pam2.tuple_type);
 
-    if (strcmp(pam1.tuple_type, PAM_PBM_TUPLETYPE) != 0 &&
-        strcmp(pam1.tuple_type, PAM_PGM_TUPLETYPE) != 0 &&
-        strcmp(pam1.tuple_type, PAM_PPM_TUPLETYPE) != 0)
+    if (!streq(pam1.tuple_type, PAM_PBM_TUPLETYPE) &&
+        !streq(pam1.tuple_type, PAM_PGM_TUPLETYPE) &&
+        !streq(pam1.tuple_type, PAM_PPM_TUPLETYPE))
         pm_error("Images are not of a PNM type.  Tuple type is '%s'",
                  pam1.tuple_type);
 }
 
 
+enum ColorSpaceId {
+    COLORSPACE_GRAYSCALE,
+    COLORSPACE_YCBCR,
+    COLORSPACE_RGB
+};
 
-static void
-psnr_color(tuple    const tuple1,
-           tuple    const tuple2,
-           double * const ySqDiffP, 
-           double * const cbSqDiffP,
-           double * const crSqDiffP) {
+typedef struct {
+
+    enum ColorSpaceId id;
+
+    unsigned int componentCt;
+
+    const char * componentName[3];
+        /* Only first 'componentCt' elements are valid */
+
+} ColorSpace;
+
+
+struct SqDiff {
+/*----------------------------------------------------------------------------
+   The square-differences of the components of two pixels, for some
+   component set.
+-----------------------------------------------------------------------------*/
+    double sqDiff[3];
+};
+
+
+
+static struct SqDiff
+zeroSqDiff() {
+
+    struct SqDiff retval;
+    unsigned int i;
+
+    for (i = 0; i < 3; ++i)
+        retval.sqDiff[i] = 0.0;
+
+    return retval;
+}
+
+
+
+static struct SqDiff
+sqDiffSum(ColorSpace    const colorSpace,
+          struct SqDiff const addend,
+          struct SqDiff const adder) {
+
+    struct SqDiff retval;
+    unsigned int i;
+
+    for (i = 0; i < colorSpace.componentCt; ++i)
+        retval.sqDiff[i] = addend.sqDiff[i] + adder.sqDiff[i];
+
+    return retval;
+}
+
+
+
+#define Y_INDEX  0
+#define CB_INDEX 1
+#define CR_INDEX 2
+
+static ColorSpace
+yCbCrColorSpace() {
+
+    ColorSpace retval;
+
+    retval.id = COLORSPACE_YCBCR;
+
+    retval.componentCt = 3;
+
+    retval.componentName[Y_INDEX]  = "Y";
+    retval.componentName[CR_INDEX] = "CR";
+    retval.componentName[CB_INDEX] = "CB";
+
+    return retval;
+}
+
+
+
+static struct SqDiff
+sqDiffYCbCr(tuple    const tuple1,
+            tuple    const tuple2) {
+
+    struct SqDiff retval;
 
     double y1, y2, cb1, cb2, cr1, cr2;
     
     pnm_YCbCrtuple(tuple1, &y1, &cb1, &cr1);
     pnm_YCbCrtuple(tuple2, &y2, &cb2, &cr2);
     
-    *ySqDiffP  = square(y1  - y2);
-    *cbSqDiffP = square(cb1 - cb2);
-    *crSqDiffP = square(cr1 - cr2);
+    retval.sqDiff[Y_INDEX]  = square(y1  - y2);
+    retval.sqDiff[CB_INDEX] = square(cb1 - cb2);
+    retval.sqDiff[CR_INDEX] = square(cr1 - cr2);
+
+    return retval;
+}
+
+
+
+#define R_INDEX 0
+#define G_INDEX 1
+#define B_INDEX 2
+
+
+
+static ColorSpace
+rgbColorSpace() {
+
+    ColorSpace retval;
+
+    retval.id = COLORSPACE_RGB;
+
+    retval.componentCt = 3;
+
+    retval.componentName[R_INDEX] = "Red";
+    retval.componentName[G_INDEX] = "Green";
+    retval.componentName[B_INDEX] = "Blue";
+
+    return retval;
+}
+
+
+
+static struct SqDiff
+sqDiffRgb(tuple    const tuple1,
+          tuple    const tuple2) {
+
+    struct SqDiff retval;
+
+    retval.sqDiff[R_INDEX] =
+        square((int)tuple1[PAM_RED_PLANE]  - (int)tuple2[PAM_RED_PLANE]);
+    retval.sqDiff[G_INDEX] =
+        square((int)tuple1[PAM_GRN_PLANE]  - (int)tuple2[PAM_GRN_PLANE]);
+    retval.sqDiff[B_INDEX] =
+        square((int)tuple1[PAM_BLU_PLANE]  - (int)tuple2[PAM_BLU_PLANE]);
+
+    return retval;
+}
+
+
+
+static ColorSpace
+grayscaleColorSpace() {
+
+    ColorSpace retval;
+
+    retval.id = COLORSPACE_GRAYSCALE;
+
+    retval.componentCt = 1;
+
+    retval.componentName[Y_INDEX]  = "luminance";
+
+    return retval;
+}
+
+
+
+static struct SqDiff
+sqDiffGrayscale(tuple    const tuple1,
+                tuple    const tuple2) {
+
+    struct SqDiff sqDiff;
+
+    sqDiff.sqDiff[Y_INDEX] = square(udiff(tuple1[0], tuple2[0]));
+
+    return sqDiff;
+}
+
+
+
+static struct SqDiff
+sumSqDiffFromRaster(struct pam * const pam1P,
+                    struct pam * const pam2P,
+                    ColorSpace   const colorSpace) {
+
+    struct SqDiff sumSqDiff;
+    tuple *tuplerow1, *tuplerow2;  /* malloc'ed */
+    unsigned int row;
+
+    tuplerow1 = pnm_allocpamrow(pam1P);
+    tuplerow2 = pnm_allocpamrow(pam2P);
+    
+    sumSqDiff = zeroSqDiff();
+
+    for (row = 0; row < pam1P->height; ++row) {
+        unsigned int col;
+        
+        pnm_readpamrow(pam1P, tuplerow1);
+        pnm_readpamrow(pam2P, tuplerow2);
+
+        assert(pam1P->width == pam2P->width);
+
+        for (col = 0; col < pam1P->width; ++col) {
+            struct SqDiff sqDiff;
+
+            switch (colorSpace.id) {
+            case COLORSPACE_GRAYSCALE:
+                sqDiff = sqDiffGrayscale(tuplerow1[col], tuplerow2[col]);
+                break;
+            case COLORSPACE_YCBCR:
+                sqDiff = sqDiffYCbCr(tuplerow1[col], tuplerow2[col]);
+                break;
+            case COLORSPACE_RGB:
+                sqDiff = sqDiffRgb(tuplerow1[col], tuplerow2[col]);
+                break;
+            }
+            sumSqDiff = sqDiffSum(colorSpace, sumSqDiff, sqDiff);
+        }
+    }
+
+    pnm_freepamrow(tuplerow1);
+    pnm_freepamrow(tuplerow2);
+
+    return sumSqDiff;
 }
 
 
 
 static void
-reportPsnr(struct pam const pam,
-           double     const ySumSqDiff, 
-           double     const crSumSqDiff,
-           double     const cbSumSqDiff,
-           char       const filespec1[],
-           char       const filespec2[]) {
+reportPsnr(struct pam    const pam,
+           struct SqDiff const sumSqDiff,
+           ColorSpace    const colorSpace,
+           const char *  const fileName1,
+           const char *  const fileName2) {
+
+    double const maxSumSqDiff = square(pam.maxval) * pam.width * pam.height;
+        /* Maximum possible sum square difference, i.e. the sum of the squares
+           of the sample differences between an entirely white image and
+           entirely black image of the given dimensions.
+        */
+
+    unsigned int i;
 
-    bool const color = streq(pam.tuple_type, PAM_PPM_TUPLETYPE);
 
     /* The PSNR is the ratio of the maximum possible mean square difference
-       to the actual mean square difference.
+       to the actual mean square difference, which is also the ratio of
+       the maximum possible sum square difference to the actual sum square
+       difference.
+   
+       Note that in the important special case that the images are
+       identical, the sum square differences are identically 0.0.
+       No precision error; no rounding error.
     */
-    double const yPsnr =
-        square(pam.maxval) / (ySumSqDiff / (pam.width * pam.height));
 
-    /* Note that in the important special case that the images are
-       identical, the sum square differences are identically 0.0.  No
-       precision error; no rounding error.
-    */
+    pm_message("PSNR between '%s' and '%s':", fileName1, fileName2);
 
-    if (color) {
-        double const cbPsnr =
-            square(pam.maxval) / (cbSumSqDiff / (pam.width * pam.height));
-        double const crPsnr =
-            square(pam.maxval) / (crSumSqDiff / (pam.width * pam.height));
+    for (i = 0; i < colorSpace.componentCt; ++i) {
+        const char * label;
 
-        pm_message("PSNR between %s and %s:", filespec1, filespec2);
-        if (ySumSqDiff > 0)
-            pm_message("Y  color component: %.2f dB", 10 * log10(yPsnr));
-        else
-            pm_message("Y color component does not differ.");
-        if (cbSumSqDiff > 0)
-            pm_message("Cb color component: %.2f dB", 10 * log10(cbPsnr));
-        else
-            pm_message("Cb color component does not differ.");
-        if (crSumSqDiff > 0)
-            pm_message("Cr color component: %.2f dB", 10 * log10(crPsnr));
-        else
-            pm_message("Cr color component does not differ.");
-    } else {
-        if (ySumSqDiff > 0)
-            pm_message("PSNR between %s and %s: %.2f dB",
-                       filespec1, filespec2, 10 * log10(yPsnr));
+        pm_asprintf(&label, "%s:", colorSpace.componentName[i]);
+
+        if (sumSqDiff.sqDiff[i] > 0)
+            pm_message("  %-6.6s %.2f dB",
+                       label,
+                       10 * log10(maxSumSqDiff/sumSqDiff.sqDiff[i]));
         else
-            pm_message("Images %s and %s don't differ.",
-                       filespec1, filespec2);
+            pm_message("  %-6.6s no difference", label);
+
+        pm_strfree(label);
     }
 }
 
 
 
 int
-main (int argc, char **argv) {
-    char *filespec1, *filespec2;  /* specs of two files to compare */
-    FILE *file1, *file2;
+main (int argc, const char **argv) {
+    FILE * if1P;
+    FILE * if2P;
     struct pam pam1, pam2;
-    bool color;
-        /* It's a color image */
-    double ySumSqDiff, crSumSqDiff, cbSumSqDiff;
-    tuple *tuplerow1, *tuplerow2;  /* malloc'ed */
-    int row;
+    ColorSpace colorSpace;
     
-    pnm_init(&argc, argv);
+    struct CmdlineInfo cmdline;
 
-    if (argc-1 < 2) 
-        pm_error("Takes two arguments:  specifications of the two files.");
-    else {
-        filespec1 = argv[1];
-        filespec2 = argv[2];
-    }
-    
-    file1 = pm_openr(filespec1);
-    file2 = pm_openr(filespec2);
+    pm_proginit(&argc, argv);
 
-    pnm_readpaminit(file1, &pam1, PAM_STRUCT_SIZE(tuple_type));
-    pnm_readpaminit(file2, &pam2, PAM_STRUCT_SIZE(tuple_type));
+    parseCommandLine(argc, argv, &cmdline);
 
-    validate_input(pam1, pam2);
+    if1P = pm_openr(cmdline.inputFile1Name);
+    if2P = pm_openr(cmdline.inputFile2Name);
 
-    if (strcmp(pam1.tuple_type, PAM_PPM_TUPLETYPE) == 0) 
-        color = TRUE;
-    else
-        color = FALSE;
+    pnm_readpaminit(if1P, &pam1, PAM_STRUCT_SIZE(tuple_type));
+    pnm_readpaminit(if2P, &pam2, PAM_STRUCT_SIZE(tuple_type));
 
-    tuplerow1 = pnm_allocpamrow(&pam1);
-    tuplerow2 = pnm_allocpamrow(&pam2);
-    
-    ySumSqDiff  = 0.0;
-    cbSumSqDiff = 0.0;
-    crSumSqDiff = 0.0;
+    validateInput(pam1, pam2);
 
-    for (row = 0; row < pam1.height; ++row) {
-        int col;
-        
-        pnm_readpamrow(&pam1, tuplerow1);
-        pnm_readpamrow(&pam2, tuplerow2);
-
-        for (col = 0; col < pam1.width; ++col) {
-            if (color) {
-                double ySqDiff, cbSqDiff, crSqDiff;
-                psnr_color(tuplerow1[col], tuplerow2[col], 
-                           &ySqDiff, &cbSqDiff, &crSqDiff);
-                ySumSqDiff  += ySqDiff;
-                cbSumSqDiff += cbSqDiff;
-                crSumSqDiff += crSqDiff;
-                
-            } else {
-                unsigned int const yDiffSq =
-                    square(udiff(tuplerow1[col][0], tuplerow2[col][0]));
-                ySumSqDiff += yDiffSq;
-            }
-        }
-    }
+    if (streq(pam1.tuple_type, PAM_PPM_TUPLETYPE)) {
+        if (cmdline.rgb)
+            colorSpace = rgbColorSpace();
+        else
+            colorSpace = yCbCrColorSpace();
+    } else
+        colorSpace = grayscaleColorSpace();
 
-    reportPsnr(pam1, ySumSqDiff, crSumSqDiff, cbSumSqDiff,
-               filespec1, filespec2);
+    {
+        struct SqDiff const sumSqDiff =
+            sumSqDiffFromRaster(&pam1, &pam2, colorSpace);
 
-    pnm_freepamrow(tuplerow1);
-    pnm_freepamrow(tuplerow2);
+        reportPsnr(pam1, sumSqDiff, colorSpace,
+                   cmdline.inputFile1Name, cmdline.inputFile2Name);
+    }
+    pm_close(if2P);
+    pm_close(if1P);
 
     return 0;
 }
+
+
+
diff --git a/analyzer/ppmhist.c b/analyzer/ppmhist.c
index f42adab5..9b526606 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 */
@@ -67,10 +69,11 @@ parseCommandLine(int argc, const char ** argv,
     /* Set defaults */
     sort_type = "frequency";
 
-    optParseOptions3(&argc, (char **)argv, opt, sizeof(opt), 0);
+    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)
@@ -101,32 +107,197 @@ parseCommandLine(int argc, const char ** argv,
 
 
 static int
-countcompare(const void *ch1, const void *ch2) {
-    return ((colorhist_vector)ch2)->value - ((colorhist_vector)ch1)->value;
+cmpUint(unsigned int const a,
+        unsigned int const b) {
+/*----------------------------------------------------------------------------
+   Return 1 if 'a' > 'b'; -1 if 'a' is < 'b'; 0 if 'a' == 'b'.
+
+   I.e. what a libc qsort() comparison function returns.
+-----------------------------------------------------------------------------*/
+    return a > b ? 1 : a < b ? -1 : 0;
 }
 
 
+
+#ifndef LITERAL_FN_DEF_MATCH
+static qsort_comparison_fn countcompare;
+#endif
+
+
 static int
-rgbcompare(const void * arg1, const void * arg2) {
+countcompare(const void * const a,
+             const void * const b) {
+/*----------------------------------------------------------------------------
+   This is a 'qsort' collation function.
+-----------------------------------------------------------------------------*/
+    const struct colorhist_item * const histItem1P = a;
+    const struct colorhist_item * const histItem2P = b;
+
+    return cmpUint(histItem2P->value, histItem1P->value);
+}
+
+
+
+#ifndef LITERAL_FN_DEF_MATCH
+static qsort_comparison_fn rgbcompare;
+#endif
 
-    colorhist_vector const ch1 = (colorhist_vector) arg1;
-    colorhist_vector const ch2 = (colorhist_vector) arg2;
+static int
+rgbcompare(const void * const a,
+           const void * const b) {
+/*----------------------------------------------------------------------------
+   This is a 'qsort' collation function.
+-----------------------------------------------------------------------------*/
+    const struct colorhist_item * const histItem1P = a;
+    const struct colorhist_item * const histItem2P = b;
 
     int retval;
 
-    retval = (PPM_GETR(ch1->color) - PPM_GETR(ch2->color));
+    retval = cmpUint(PPM_GETR(histItem1P->color),
+                     PPM_GETR(histItem2P->color));
     if (retval == 0) {
-        retval = (PPM_GETG(ch1->color) - PPM_GETG(ch2->color));
+        retval = cmpUint(PPM_GETG(histItem1P->color),
+                         PPM_GETG(histItem2P->color));
         if (retval == 0)
-            retval = (PPM_GETB(ch1->color) - PPM_GETB(ch2->color));
+            retval = cmpUint(PPM_GETB(histItem1P->color),
+                             PPM_GETB(histItem2P->color));
     }
     return retval;
 }
 
 
 
+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[],
@@ -136,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.
         */
@@ -154,31 +325,31 @@ colornameLabel(pixel        const color,
     colorIndex = ppm_findclosestcolor(dictColors, nDictColor, &color255);
 
     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);
@@ -190,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 = "";
@@ -221,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);
 
@@ -247,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 ? '#' : ' ';
@@ -273,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);
@@ -291,5 +546,11 @@ main(int argc, const char *argv[]) {
 
     ppm_freecolorhist(chv);
 
+    if (chvInvalid)
+        ppm_freecolorhist(chvInvalid);
+
     return 0;
 }
+
+
+