promote Advanced to Stable

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

11 files changed, 1145 insertions, 555 deletions

diff --git a/other/Makefile b/other/Makefile
index aa0d75bb..b3884666 100644
--- a/other/Makefile
+++ b/other/Makefile
@@ -24,10 +24,14 @@ endif
 # build.
 
 PORTBINARIES = pamarith pambayer pamchannel pamdepth \
-	pamendian pamexec pamfix pamlookup pampick pamsplit \
+	pamendian pamfix pamlookup pampick pamsplit \
 	pamstack pamsummcol pamunlookup pamvalidate pnmcolormap \
 	ppmdcfont ppmddumpfont ppmdmkfont 
 
+ifneq ($(DONT_HAVE_PROCESS_MGMT),Y)
+  PORTBINARIES += pamexec
+endif
+
 ifneq ($(LINUXSVGALIB),NONE)
   PORTBINARIES += ppmsvgalib
 endif
diff --git a/other/pamarith.c b/other/pamarith.c
index 3f44dc73..31f52a59 100644
--- a/other/pamarith.c
+++ b/other/pamarith.c
@@ -1,14 +1,18 @@
 #include <assert.h>
 #include <string.h>
 #include <limits.h>
+#include <math.h>
 
 #include "pm_c_util.h"
 #include "mallocvar.h"
+#include "nstring.h"
 #include "shhopt.h"
 #include "pam.h"
 
-enum function {FN_ADD, FN_SUBTRACT, FN_MULTIPLY, FN_DIVIDE, FN_DIFFERENCE,
-               FN_MINIMUM, FN_MAXIMUM, FN_MEAN, FN_COMPARE,
+static double const EPSILON = 1.0e-5;
+
+enum Function {FN_ADD, FN_SUBTRACT, FN_MULTIPLY, FN_DIVIDE, FN_DIFFERENCE,
+               FN_MINIMUM, FN_MAXIMUM, FN_MEAN, FN_EQUAL, FN_COMPARE,
                FN_AND, FN_OR, FN_NAND, FN_NOR, FN_XOR,
                FN_SHIFTLEFT, FN_SHIFTRIGHT
               };
@@ -16,16 +20,15 @@ enum function {FN_ADD, FN_SUBTRACT, FN_MULTIPLY, FN_DIVIDE, FN_DIFFERENCE,
 
 
 static bool
-isDyadic(enum function const function) {
+isDyadic(enum Function const function) {
 
     bool retval;
-    
+
     switch (function) {
     case FN_ADD:
     case FN_MULTIPLY:
     case FN_MINIMUM:
     case FN_MAXIMUM:
-    case FN_MEAN:
     case FN_AND:
     case FN_NAND:
     case FN_OR:
@@ -37,8 +40,10 @@ isDyadic(enum function const function) {
     case FN_DIFFERENCE:
     case FN_COMPARE:
     case FN_DIVIDE:
+    case FN_MEAN:
     case FN_SHIFTLEFT:
     case FN_SHIFTRIGHT:
+    case FN_EQUAL:
         retval = TRUE;
         break;
     }
@@ -51,9 +56,10 @@ struct CmdlineInfo {
     /* All the information the user supplied in the command line,
        in a form easy for the program to use.
     */
-    enum function function;
+    enum Function function;
     unsigned int operandCt;
     const char ** operandFileNames;
+    double closeness;
 };
 
 
@@ -71,12 +77,14 @@ parseCommandLine(int argc, const char ** const argv,
     optStruct3 opt;
 
     unsigned int option_def_index;
-    
+
     unsigned int addSpec, subtractSpec, multiplySpec, divideSpec,
         differenceSpec,
-        minimumSpec, maximumSpec, meanSpec, compareSpec,
+        minimumSpec, maximumSpec, meanSpec, equalSpec, compareSpec,
         andSpec, orSpec, nandSpec, norSpec, xorSpec,
-        shiftleftSpec, shiftrightSpec;
+        shiftleftSpec, shiftrightSpec, closenessSpec;
+
+    float closeness;
 
     MALLOCARRAY_NOFAIL(option_def, 100);
 
@@ -89,6 +97,7 @@ parseCommandLine(int argc, const char ** const argv,
     OPTENT3(0, "minimum",     OPT_FLAG,   NULL, &minimumSpec,    0);
     OPTENT3(0, "maximum",     OPT_FLAG,   NULL, &maximumSpec,    0);
     OPTENT3(0, "mean",        OPT_FLAG,   NULL, &meanSpec,       0);
+    OPTENT3(0, "equal",       OPT_FLAG,   NULL, &equalSpec,      0);
     OPTENT3(0, "compare",     OPT_FLAG,   NULL, &compareSpec,    0);
     OPTENT3(0, "and",         OPT_FLAG,   NULL, &andSpec,        0);
     OPTENT3(0, "or",          OPT_FLAG,   NULL, &orSpec,         0);
@@ -97,6 +106,7 @@ parseCommandLine(int argc, const char ** const argv,
     OPTENT3(0, "xor",         OPT_FLAG,   NULL, &xorSpec,        0);
     OPTENT3(0, "shiftleft",   OPT_FLAG,   NULL, &shiftleftSpec,  0);
     OPTENT3(0, "shiftright",  OPT_FLAG,   NULL, &shiftrightSpec, 0);
+    OPTENT3(0, "closeness",   OPT_FLOAT,  &closeness, &closenessSpec, 0);
 
     opt.opt_table = option_def;
     opt.short_allowed = FALSE;  /* We have no short (old-fashioned) options */
@@ -106,7 +116,7 @@ parseCommandLine(int argc, const char ** const argv,
         /* Uses and sets argc, argv, and some of *cmdlineP and others. */
 
     if (addSpec + subtractSpec + multiplySpec + divideSpec + differenceSpec +
-        minimumSpec + maximumSpec + meanSpec + compareSpec +
+        minimumSpec + maximumSpec + meanSpec + equalSpec + compareSpec +
         andSpec + orSpec + nandSpec + norSpec + xorSpec +
         shiftleftSpec + shiftrightSpec > 1)
         pm_error("You may specify only one function");
@@ -127,6 +137,8 @@ parseCommandLine(int argc, const char ** const argv,
         cmdlineP->function = FN_MAXIMUM;
     else if (meanSpec)
         cmdlineP->function = FN_MEAN;
+    else if (equalSpec)
+        cmdlineP->function = FN_EQUAL;
     else if (compareSpec)
         cmdlineP->function = FN_COMPARE;
     else if (andSpec)
@@ -146,21 +158,35 @@ parseCommandLine(int argc, const char ** const argv,
     else
         pm_error("You must specify a function (e.g. '-add')");
 
+    if (closenessSpec) {
+        if (cmdlineP->function != FN_EQUAL)
+            pm_error("-closeness is valid only with -equal");
+        else {
+            if (closeness < 0 || closeness > 100)
+                pm_error("-closeness value %f is not a valid percentage",
+                         closeness);
+            cmdlineP->closeness = (double)closeness/100;
+        }
+    } else
+        cmdlineP->closeness = EPSILON;
+
     if (argc-1 < 2)
         pm_error("You must specify at least two arguments: the files which "
                  "are the operands of the function.  You specified %u",
                  argc-1);
     else {
         if (isDyadic(cmdlineP->function) && argc-1 > 2)
-            pm_error("You specified %u arguments, but a dyadic function.  "
-                     "For a dyadic function, you must specify 2 arguments:  "
-                     "the operands of the function", argc-1);
+            pm_error("You specified %u arguments, but a function "
+                     "which is not associative and commutative.  "
+                     "Unless a function is associative and commutative "
+                     "(like -add), you must specify exactly two arguments. ",
+                     argc-1);
         else {
             cmdlineP->operandCt = argc-1;
             cmdlineP->operandFileNames = &argv[1];
         }
     }
-}        
+}
 
 
 
@@ -179,11 +205,69 @@ enum category {
 
 
 
+static void
+validateSameWidthHeight(const struct pam * const inpam  /* array */,
+                        unsigned int       const operandCt) {
+    unsigned int i;
+
+    for (i = 1; i < operandCt; ++i) {
+
+        if (i > 0) {
+            if (inpam[i].width  != inpam[0].width ||
+                inpam[i].height != inpam[0].height) {
+                pm_error("The images must be the same width and height.  "
+                         "The first is %ux%ux%u, but another is %ux%ux%u",
+                         inpam[0].width, inpam[0].height, inpam[0].depth,
+                         inpam[i].width, inpam[i].height, inpam[i].depth);
+            }
+        }
+    }
+}
+
+
+
+static void
+validateCompatibleDepth(const struct pam * const inpam  /* array */,
+                        unsigned int       const operandCt) {
+
+    unsigned int i;
+    bool         haveNonUnityDepth;
+    unsigned int nonUnityDepth;
+
+    for (i = 0, haveNonUnityDepth = false; i < operandCt; ++i) {
+        if (inpam[i].depth != 1) {
+            if (haveNonUnityDepth) {
+                if (inpam[i].depth != nonUnityDepth)
+                    pm_error("The images must have the same depth "
+                             "or depth 1.  "
+                             "But one has depth %u and another %u",
+                             nonUnityDepth, inpam[i].depth);
+            } else {
+                haveNonUnityDepth = true;
+                nonUnityDepth = inpam[i].depth;
+            }
+        }
+    }
+}
+
+
+
+static void
+validateConsistentDimensions(const struct pam * const inpam  /* array */,
+                             unsigned int       const operandCt) {
+
+    validateSameWidthHeight(inpam, operandCt);
+
+    validateCompatibleDepth(inpam, operandCt);
+}
+
+
+
 static enum category
-functionCategory(enum function const function) {
+functionCategory(enum Function const function) {
 
     enum category retval;
-    
+
     switch (function) {
     case FN_ADD:
     case FN_SUBTRACT:
@@ -191,6 +275,7 @@ functionCategory(enum function const function) {
     case FN_MINIMUM:
     case FN_MAXIMUM:
     case FN_MEAN:
+    case FN_EQUAL:
     case FN_COMPARE:
     case FN_MULTIPLY:
     case FN_DIVIDE:
@@ -237,58 +322,147 @@ outFmtForCompare(int const format1,
 
 
 
+static unsigned int
+maxDepth(const struct pam * const pam,
+         unsigned int       const pamCt) {
+
+    unsigned int maxDepthSoFar;
+    unsigned int i;
+
+    assert(pamCt >= 1);
+
+    maxDepthSoFar = pam[0].depth;
+    for (i = 1; i < pamCt; ++i) {
+        if (pam[i].depth > maxDepthSoFar)
+            maxDepthSoFar = pam[i].depth;
+    }
+    return maxDepthSoFar;
+}
+
+
+
+static int
+maxFormat(const struct pam * const pam,
+          unsigned int       const pamCt) {
+
+    int maxFormatSoFar;
+    unsigned int i;
+
+    assert(pamCt >= 1);
+
+    maxFormatSoFar = pam[0].format;
+    for (i = 1; i < pamCt; ++i) {
+        if (pam[i].format > maxFormatSoFar)
+            maxFormatSoFar = pam[i].format;
+    }
+    return maxFormatSoFar;
+}
+
+
+
+static sample
+maxMaxval(const struct pam * const pam,
+          unsigned int       const pamCt) {
+
+    sample maxMaxvalSoFar;
+    unsigned int i;
+
+    assert(pamCt >= 1);
+
+    maxMaxvalSoFar = pam[0].maxval;
+    for (i = 1; i < pamCt; ++i) {
+        if (pam[i].maxval > maxMaxvalSoFar)
+            maxMaxvalSoFar = pam[i].maxval;
+    }
+    return maxMaxvalSoFar;
+}
+
+
+
+static bool
+maxvalsAreEqual(const struct pam * const pam,
+                unsigned int       const pamCt) {
+
+    bool equalSoFar;
+    unsigned int i;
+
+    assert(pamCt >= 1);
+
+    equalSoFar = true;
+
+    for (i = 1; i < pamCt; ++i) {
+        if (pam[i].maxval != pam[0].maxval)
+            equalSoFar = false;
+    }
+    return equalSoFar;
+}
+
+
+
 static void
-computeOutputType(struct pam *  const outpamP,
-                  struct pam    const inpam1,
-                  struct pam    const inpam2,
-                  enum function const function) {
+computeOutputType(struct pam *       const outpamP,
+                  const struct pam * const inpam,   /* array */
+                  unsigned int       const operandCt,
+                  enum Function      const function) {
+
+    assert(operandCt >= 1);
 
     outpamP->size        = sizeof(struct pam);
     outpamP->len         = PAM_STRUCT_SIZE(tuple_type);
     outpamP->file        = stdout;
     outpamP->plainformat = FALSE;
-    outpamP->height      = inpam1.height;
-    outpamP->width       = inpam1.width;
-    outpamP->depth       = MAX(inpam1.depth, inpam2.depth);
+    outpamP->height      = inpam[0].height;
+    outpamP->width       = inpam[0].width;
+    outpamP->depth       = maxDepth(inpam, operandCt);
 
-    if (function == FN_COMPARE)
-        outpamP->format = outFmtForCompare(inpam1.format, inpam2.format);
-    else
-        outpamP->format = MAX(inpam1.format, inpam2.format);
+
+    if (function == FN_COMPARE) {
+        assert(operandCt == 2);
+        outpamP->format = outFmtForCompare(inpam[0].format, inpam[1].format);
+    } else
+        outpamP->format = maxFormat(inpam, operandCt);
 
     switch (functionCategory(function)) {
     case CATEGORY_FRACTIONAL_ARITH:
         if (function == FN_COMPARE)
             outpamP->maxval = 2;
+        else if (function == FN_EQUAL)
+            outpamP->maxval = 1;
         else
-            outpamP->maxval = MAX(inpam1.maxval, inpam2.maxval);
+            outpamP->maxval = maxMaxval(inpam, operandCt);
         break;
     case CATEGORY_BITSTRING:
-        if (inpam2.maxval != inpam1.maxval)
+        if (!maxvalsAreEqual(inpam, operandCt))
             pm_error("For a bit string operation, the maxvals of the "
-                     "left and right image must be the same.  You have "
-                     "left=%u and right=%u", 
-                     (unsigned)inpam1.maxval, (unsigned)inpam2.maxval);
+                     "operand images must be the same.  Yours differ");
 
-        if (pm_bitstomaxval(pm_maxvaltobits(inpam1.maxval)) != inpam1.maxval)
+        if (pm_bitstomaxval(pm_maxvaltobits(inpam[0].maxval)) !=
+            inpam[0].maxval)
             pm_error("For a bit string operation, the maxvals of the inputs "
                      "must be a full binary count, i.e. a power of two "
                      "minus one such as 0xff.  You have 0x%x",
-                     (unsigned)inpam1.maxval);
+                     (unsigned)inpam[0].maxval);
 
-        outpamP->maxval = inpam1.maxval;
+        outpamP->maxval = inpam[0].maxval;
         break;
     case CATEGORY_SHIFT:
-        if (pm_bitstomaxval(pm_maxvaltobits(inpam1.maxval)) != inpam1.maxval)
+        if (pm_bitstomaxval(pm_maxvaltobits(inpam[0].maxval)) !=
+            inpam[0].maxval)
             pm_error("For a bit shift operation, the maxval of the left "
                      "input image "
                      "must be a full binary count, i.e. a power of two "
                      "minus one such as 0xff.  You have 0x%x",
-                     (unsigned)inpam1.maxval);
-        outpamP->maxval = inpam1.maxval;
+                     (unsigned)inpam[0].maxval);
+        outpamP->maxval = inpam[0].maxval;
     }
     outpamP->bytes_per_sample = (pm_maxvaltobits(outpamP->maxval)+7)/8;
-    strcpy(outpamP->tuple_type, inpam1.tuple_type);
+
+    if (outpamP->maxval > 1 &&
+        strneq(inpam[0].tuple_type, "BLACKANDWHITE", 13)) {
+
+        strcpy(outpamP->tuple_type, "");
+    } else
+        strcpy(outpamP->tuple_type, inpam[0].tuple_type);
 }
 
 
@@ -311,6 +485,21 @@ samplenSum(samplen      const operands[],
 
 
 static samplen
+samplenProduct(samplen      const operands[],
+               unsigned int const operandCt) {
+
+    unsigned int i;
+    double product;
+
+    for (i = 1, product = operands[0]; i < operandCt; ++i)
+        product *= operands[i];
+
+    return product;
+}
+
+
+
+static samplen
 samplenMin(samplen      const operands[],
            unsigned int const operandCt) {
 
@@ -358,24 +547,28 @@ samplenMean(samplen      const operands[],
 
 
 static samplen
-samplenProduct(samplen      const operands[],
-               unsigned int const operandCt) {
+samplenEqual(samplen      const operands[],
+             unsigned int const operandCt,
+             double       const closeness) {
 
     unsigned int i;
-    double product;
+    bool allEqual;
 
-    for (i = 1, product = operands[0]; i < operandCt; ++i)
-        product *= operands[i];
+    for (i = 1, allEqual = true; i < operandCt; ++i) {
+        if (fabs(operands[i]- operands[0]) > closeness)
+            allEqual = false;
+    }
 
-    return product;
+    return allEqual ? 1.0 : 0.0;
 }
 
 
 
 static samplen
-applyNormalizedFunction(enum function const function,
+applyNormalizedFunction(enum Function const function,
                         samplen       const operands[],
-                        unsigned int  const operandCt) {
+                        unsigned int  const operandCt,
+                        double        const closeness) {
 
     samplen result;
 
@@ -394,7 +587,7 @@ applyNormalizedFunction(enum function const function,
         operands[0] / operands[1] : 1.;
         break;
     case FN_DIFFERENCE:
-        result = operands[0] > operands[1] ? 
+        result = operands[0] > operands[1] ?
             operands[0] - operands[1] : operands[1] - operands[0];
         break;
     case FN_MINIMUM:
@@ -406,8 +599,11 @@ applyNormalizedFunction(enum function const function,
     case FN_MEAN:
         result = samplenMean(operands, operandCt);
         break;
+    case FN_EQUAL:
+        result = samplenEqual(operands, operandCt, closeness);
+        break;
     case FN_COMPARE:
-        result = 
+        result =
             operands[0] > operands[1] ?
             1. : operands[0] < operands[1] ?
             0. : .5;
@@ -423,16 +619,11 @@ applyNormalizedFunction(enum function const function,
 
 
 static void
-doNormalizedArith(struct pam *  const inpam1P,
-                  struct pam *  const inpam2P,
-                  struct pam *  const outpamP,
-                  enum function const function) {
-
-    /* Some of the logic in this subroutine is designed for future
-       expansion into non-dyadic computations.  But for now, all
-       computations have exactly two operands.
-    */
-    unsigned int const operandCt = 2;
+doNormalizedArith(const struct pam * const inpam,  /* array */
+                  unsigned int       const operandCt,
+                  const struct pam * const outpamP,
+                  enum Function      const function,
+                  double             const closeness) {
 
     tuplen ** tuplerown;
         /* tuplerown[0] is the current row in the first operand image */
@@ -444,38 +635,43 @@ doNormalizedArith(struct pam *  const inpam1P,
            computation
         */
     unsigned int * plane;
-        /* plane[0] is the plane number in the first operand image for 
+        /* plane[0] is the plane number in the first operand image for
            the current one-sample computation.  plane[1] is the plane number
            in the second operand image, etc.
          */
+    unsigned int i;
 
     MALLOCARRAY_NOFAIL(operands, operandCt);
     MALLOCARRAY_NOFAIL(plane, operandCt);
     MALLOCARRAY_NOFAIL(tuplerown, operandCt);
 
-    tuplerown[0] = pnm_allocpamrown(inpam1P);
-    tuplerown[1] = pnm_allocpamrown(inpam2P);
+    for (i = 0; i < operandCt; ++i)
+        tuplerown[i] = pnm_allocpamrown(&inpam[i]);
     tuplerownOut = pnm_allocpamrown(outpamP);
 
     for (row = 0; row < outpamP->height; ++row) {
+        unsigned int i;
         unsigned int col;
-        pnm_readpamrown(inpam1P, tuplerown[0]);
-        pnm_readpamrown(inpam2P, tuplerown[1]);
-        
+
+        for (i = 0; i < operandCt; ++i)
+            pnm_readpamrown(&inpam[i], tuplerown[i]);
+
         for (col = 0; col < outpamP->width; ++col) {
             unsigned int outplane;
-            
+
             for (outplane = 0; outplane < outpamP->depth; ++outplane) {
+                unsigned int i;
                 unsigned int op;
 
-                plane[0] = MIN(outplane, inpam1P->depth-1);
-                plane[1] = MIN(outplane, inpam2P->depth-1);
+                for (i = 0; i < operandCt; ++i)
+                    plane[i] = MIN(outplane, inpam[i].depth-1);
 
                 for (op = 0; op < operandCt; ++op)
                     operands[op] = tuplerown[op][col][plane[op]];
 
-                tuplerownOut[col][outplane] = 
-                    applyNormalizedFunction(function, operands, operandCt); 
+                tuplerownOut[col][outplane] =
+                    applyNormalizedFunction(function, operands, operandCt,
+                                            closeness);
                 assert(tuplerownOut[col][outplane] >= 0.);
                 assert(tuplerownOut[col][outplane] <= 1.);
             }
@@ -483,8 +679,8 @@ doNormalizedArith(struct pam *  const inpam1P,
         pnm_writepamrown(outpamP, tuplerownOut);
     }
 
-    pnm_freepamrown(tuplerown[0]);
-    pnm_freepamrown(tuplerown[1]);
+    for (i = 0; i < operandCt; ++i)
+        pnm_freepamrown(tuplerown[i]);
     free(tuplerown);
     pnm_freepamrown(tuplerownOut);
     free(plane);
@@ -512,6 +708,22 @@ sampleSum(sample       const operands[],
 
 
 static sample
+sampleProduct(sample       const operands[],
+              unsigned int const operandCt,
+              sample       const maxval) {
+
+    unsigned int i;
+    double product;
+
+    for (i = 0, product = 1.0; i < operandCt; ++i) {
+        product *= ((double)operands[i]/maxval);
+    }
+    return (sample)(product * maxval + 0.5);
+}
+
+
+
+static sample
 sampleMin(sample       const operands[],
           unsigned int const operandCt) {
 
@@ -561,17 +773,18 @@ sampleMean(sample       const operands[],
 
 
 static sample
-sampleProduct(sample       const operands[],
-              unsigned int const operandCt,
-              sample       const maxval) {
+sampleEqual(sample       const operands[],
+            unsigned int const operandCt,
+            sample       const maxval) {
 
     unsigned int i;
-    double product;
+    bool allEqual;
 
-    for (i = 0, product = 1.0; i < operandCt; ++i) {
-        product *= ((double)operands[i]/maxval);
+    for (i = 1, allEqual = true; i < operandCt; ++i) {
+        if (operands[i] != operands[0])
+            allEqual = false;
     }
-    return (sample)(product * maxval + 0.5);
+    return allEqual ? maxval : 0;
 }
 
 
@@ -654,7 +867,7 @@ sampleXor(sample       const operands[],
 
 
 static sample
-applyUnNormalizedFunction(enum function const function,
+applyUnNormalizedFunction(enum Function const function,
                           sample        const operands[],
                           unsigned int  const operandCt,
                           sample        const maxval) {
@@ -694,6 +907,9 @@ applyUnNormalizedFunction(enum function const function,
     case FN_MEAN:
         result = sampleMean(operands, operandCt);
         break;
+    case FN_EQUAL:
+        result = sampleEqual(operands, operandCt, maxval);
+        break;
     case FN_COMPARE:
         result = operands[0] > operands[1] ?
             2 : operands[0] < operands[1] ? 0 : 1;
@@ -738,21 +954,15 @@ applyUnNormalizedFunction(enum function const function,
 
 
 static void
-doUnNormalizedArith(struct pam *  const inpam1P,
-                    struct pam *  const inpam2P,
-                    struct pam *  const outpamP,
-                    enum function const function) {
+doUnNormalizedArith(const struct pam * const inpam,  /* array */
+                    unsigned int       const operandCt,
+                    const struct pam * const outpamP,
+                    enum Function      const function) {
 /*----------------------------------------------------------------------------
    Take advantage of the fact that both inputs and the output use the same
    maxval to do the computation without time-consuming normalization of
    sample values.
 -----------------------------------------------------------------------------*/
-    /* Some of the logic in this subroutine is designed for future
-       expansion into non-dyadic computations.  But for now, all
-       computations have exactly two operands.
-    */
-    unsigned int const operandCt = 2;
-
     sample const maxval = outpamP->maxval;
 
     tuple ** tuplerow;
@@ -765,42 +975,42 @@ doUnNormalizedArith(struct pam *  const inpam1P,
            computation
         */
     unsigned int * plane;
-        /* plane[0] is the plane number in the first operand image for 
+        /* plane[0] is the plane number in the first operand image for
            the current one-sample computation.  plane[1] is the plane number
            in the second operand image, etc.
          */
-
-    /* Input conditions: */
-    assert(inpam1P->maxval == maxval);
-    assert(inpam2P->maxval == maxval);
-    assert(outpamP->maxval == maxval);
+    unsigned int i;
 
     MALLOCARRAY_NOFAIL(operands, operandCt);
     MALLOCARRAY_NOFAIL(plane, operandCt);
     MALLOCARRAY_NOFAIL(tuplerow, operandCt);
 
-    tuplerow[0]   = pnm_allocpamrow(inpam1P);
-    tuplerow[1]   = pnm_allocpamrow(inpam2P);
+    for (i = 0; i < operandCt; ++i)
+        tuplerow[i]   = pnm_allocpamrow(&inpam[i]);
+
     tuplerowOut = pnm_allocpamrow(outpamP);
 
     for (row = 0; row < outpamP->height; ++row) {
+        unsigned int i;
         unsigned int col;
-        pnm_readpamrow(inpam1P, tuplerow[0]);
-        pnm_readpamrow(inpam2P, tuplerow[1]);
-        
+
+        for (i = 0; i < operandCt; ++i)
+            pnm_readpamrow(&inpam[i], tuplerow[i]);
+
         for (col = 0; col < outpamP->width; ++col) {
             unsigned int outplane;
-            
+
             for (outplane = 0; outplane < outpamP->depth; ++outplane) {
                 unsigned int op;
+                unsigned int i;
 
-                plane[0] = MIN(outplane, inpam1P->depth-1);
-                plane[1] = MIN(outplane, inpam2P->depth-1);
+                for (i = 0; i < operandCt; ++i)
+                    plane[i] = MIN(outplane, inpam[i].depth-1);
 
                 for (op = 0; op < operandCt; ++op)
                     operands[op] = tuplerow[op][col][plane[op]];
 
-                tuplerowOut[col][outplane] = 
+                tuplerowOut[col][outplane] =
                     applyUnNormalizedFunction(function, operands, operandCt,
                                               maxval);
 
@@ -811,8 +1021,8 @@ doUnNormalizedArith(struct pam *  const inpam1P,
         pnm_writepamrow(outpamP, tuplerowOut);
     }
 
-    pnm_freepamrow(tuplerow[0]);
-    pnm_freepamrow(tuplerow[1]);
+    for (i = 0; i < operandCt; ++i)
+        pnm_freepamrow(tuplerow[i]);
     free(tuplerow);
     pnm_freepamrow(tuplerowOut);
     free(plane);
@@ -825,57 +1035,60 @@ int
 main(int argc, const char *argv[]) {
 
     struct CmdlineInfo cmdline;
-    struct pam inpam1;
-    struct pam inpam2;
-    struct pam outpam;
-    FILE * if1P;
-    FILE * if2P;
-    
+    struct pam * inpam;  /* malloc'ed array */
+    FILE **      ifP;    /* malloc'ed array */
+    struct pam   outpam;
+
     pm_proginit(&argc, argv);
 
     parseCommandLine(argc, argv, &cmdline);
 
-    if (cmdline.operandCt != 2)
-        /* Code for > 2 operands not written yet */
-        pm_error("You specified %u operands.  We understand only 2.",
-                 cmdline.operandCt);
+    MALLOCARRAY(inpam, cmdline.operandCt);
+    MALLOCARRAY(ifP,   cmdline.operandCt);
 
-    if1P = pm_openr(cmdline.operandFileNames[0]);
-    if2P = pm_openr(cmdline.operandFileNames[1]);
-
-    pnm_readpaminit(if1P, &inpam1, PAM_STRUCT_SIZE(tuple_type));
-    pnm_readpaminit(if2P, &inpam2, PAM_STRUCT_SIZE(tuple_type));
+    if (!inpam || !ifP)
+        pm_error("Failed to allocate arrays for %u operands",
+                 cmdline.operandCt);
 
-    if (inpam1.width != inpam2.width || inpam1.height != inpam2.height)
-        pm_error("The two images must be the same width and height.  "
-                 "The first is %ux%ux%u, but the second is %ux%ux%u",
-                 inpam1.width, inpam1.height, inpam1.depth,
-                 inpam2.width, inpam2.height, inpam2.depth);
+    {
+        unsigned int i;
+        for (i = 0; i < cmdline.operandCt; ++i) {
+            ifP[i] = pm_openr(cmdline.operandFileNames[i]);
 
-    if (inpam1.depth != 1 && inpam2.depth != 1 && inpam1.depth != inpam2.depth)
-        pm_error("The two images must have the same depth or one of them "
-                 "must have depth 1.  The first has depth %u and the second "
-                 "has depth %u", inpam1.depth, inpam2.depth);
+            pnm_readpaminit(ifP[i], &inpam[i], PAM_STRUCT_SIZE(tuple_type));
+        }
+    }
+    validateConsistentDimensions(inpam, cmdline.operandCt);
 
-    computeOutputType(&outpam, inpam1, inpam2, cmdline.function);
+    computeOutputType(&outpam, inpam, cmdline.operandCt, cmdline.function);
 
     pnm_writepaminit(&outpam);
 
-    switch (functionCategory(cmdline.function)) {    
+    switch (functionCategory(cmdline.function)) {
     case CATEGORY_FRACTIONAL_ARITH:
-        if (inpam1.maxval == inpam2.maxval && inpam2.maxval == outpam.maxval)
-            doUnNormalizedArith(&inpam1, &inpam2, &outpam, cmdline.function);
+        if (maxvalsAreEqual(inpam, cmdline.operandCt) &&
+            inpam[0].maxval == outpam.maxval)
+            doUnNormalizedArith(inpam, cmdline.operandCt, &outpam,
+                                cmdline.function);
         else
-            doNormalizedArith(&inpam1, &inpam2, &outpam, cmdline.function);
+            doNormalizedArith(inpam, cmdline.operandCt, &outpam,
+                              cmdline.function, cmdline.closeness);
         break;
     case CATEGORY_BITSTRING:
     case CATEGORY_SHIFT:
-        doUnNormalizedArith(&inpam1, &inpam2, &outpam, cmdline.function);
+        doUnNormalizedArith(inpam, cmdline.operandCt, &outpam,
+                            cmdline.function);
         break;
     }
 
-    pm_close(if1P);
-    pm_close(if2P);
-    
+    {
+        unsigned int i;
+        for (i = 0; i < cmdline.operandCt; ++i)
+            pm_close(ifP[i]);
+    }
+
     return 0;
 }
+
+
+
diff --git a/other/pambayer.c b/other/pambayer.c
index 9cffc8f0..599dc11e 100644
--- a/other/pambayer.c
+++ b/other/pambayer.c
@@ -18,7 +18,7 @@
   USA
 
   Copyright Alexandre Becoulet <diaxen AT free DOT fr>
-  
+
   Completely rewritten for Netpbm by Bryan Henderson August 2005.
 */
 
@@ -32,24 +32,24 @@
 #include "nstring.h"
 
 
-enum bayerType {
+enum BayerType {
     BAYER1,
     BAYER2,
     BAYER3,
     BAYER4
 };
 
-struct cmdlineInfo {
-    const char * inputFilespec;
-    enum bayerType bayerType;
-    unsigned int nointerpolate;
+struct CmdlineInfo {
+    const char *   inputFilespec;
+    enum BayerType bayerType;
+    unsigned int   nointerpolate;
 };
 
 
 
 static void
 parseCommandLine(int argc, const char ** argv,
-                 struct cmdlineInfo * 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.
@@ -66,10 +66,10 @@ parseCommandLine(int argc, const char ** argv,
     MALLOCARRAY_NOFAIL(option_def, 100);
 
     option_def_index = 0;   /* incremented by OPTENT3 */
-    OPTENT3(0, "type",     OPT_UINT, &type,
-            &typeSpec, 0);
+    OPTENT3(0, "type",          OPT_UINT, &type,
+            &typeSpec,                      0);
     OPTENT3(0, "nointerpolate", OPT_FLAG, NULL,
-            &cmdlineP->nointerpolate, 0);
+            &cmdlineP->nointerpolate,       0);
 
     opt.opt_table = option_def;
     opt.short_allowed = FALSE;  /* We have no short (old-fashioned) options */
@@ -104,47 +104,60 @@ parseCommandLine(int argc, const char ** argv,
 
 static void
 clearTuples(const struct pam * const pamP,
-	    tuple **           const outtuples) {
+            tuple **           const outtuples) {
+/*----------------------------------------------------------------------------
+  Make tuples at the edge that may not get set to anything by the normal
+  computation of the bayer pattern black.
+-----------------------------------------------------------------------------*/
+    if (pamP->height <= 4 || pamP->width <= 4) {
+        unsigned int row;
+
+        for (row = 0; row < pamP->height; ++row) {
+            unsigned int col;
+            for (col = 0; col < pamP->width; ++col) {
+                unsigned int plane;
+                for (plane = 0; plane < pamP->depth; ++plane)
+                    outtuples[row][col][plane] = 0;
+            }
+        }
+    } else {
+        unsigned int col;
+        unsigned int row;
 
-    unsigned int row;
-    unsigned int col;
-    unsigned int plane;
+        for (col = 0; col < pamP->width; ++col) {
+            unsigned int plane;
 
-    if(pamP->height <= 4 || pamP->width <= 4) {
-        for(row=0; row < pamP->height; ++row)
-          for(col=0; col < pamP->width; ++col)
-            for (plane=0; plane < pamP->depth; ++plane)
-              outtuples[row][col][plane] = 0;
-    }
-    else {
-        for(col = 0; col < pamP->width; ++col)
             for (plane = 0; plane < pamP->depth; ++plane) {
                 outtuples[0][col][plane] = 0;
                 outtuples[1][col][plane] = 0;
                 outtuples[pamP->height-2][col][plane] = 0;
                 outtuples[pamP->height-1][col][plane] = 0;
-          }
+            }
 
-        for(row = 2; row < pamP->height - 2; ++row)
-            for (plane = 0; plane < pamP->depth; ++plane) {
-                outtuples[row][0][plane] = 0;
-                outtuples[row][1][plane] = 0;
-                outtuples[row][pamP->width-2][plane] = 0;
-                outtuples[row][pamP->width-1][plane] = 0;
+            for (row = 2; row < pamP->height - 2; ++row) {
+                unsigned int plane;
+
+                for (plane = 0; plane < pamP->depth; ++plane) {
+                    outtuples[row][0][plane] = 0;
+                    outtuples[row][1][plane] = 0;
+                    outtuples[row][pamP->width-2][plane] = 0;
+                    outtuples[row][pamP->width-1][plane] = 0;
+                }
             }
+        }
     }
 }
 
 
 
 static void
-calc_4(const struct pam * const pamP,
-       tuple **           const intuples,
-       tuple **           const outtuples,
-       unsigned int       const plane,
-       bool               const noInterpolation,
-       unsigned int       const xoffset,
-       unsigned int       const yoffset) {
+calc4(const struct pam * const pamP,
+      tuple **           const intuples,
+      tuple **           const outtuples,
+      unsigned int       const plane,
+      bool               const noInterpolation,
+      unsigned int       const xoffset,
+      unsigned int       const yoffset) {
 /*----------------------------------------------------------------------------
     X . X
     . . .
@@ -158,7 +171,7 @@ calc_4(const struct pam * const pamP,
   (even/odd is with respect to ('xoffset', 'yoffset')).
 -----------------------------------------------------------------------------*/
     unsigned int row;
-    
+
     /* Do the even rows -- the even column pixels get copied from the input,
        while the odd column pixels get the mean of adjacent even ones
     */
@@ -177,7 +190,7 @@ calc_4(const struct pam * const pamP,
     for (row = yoffset; row + 2 < pamP->height; row += 2) {
         unsigned int col;
         for (col = xoffset; col < pamP->width; ++col) {
-            outtuples[row + 1][col][plane] = 
+            outtuples[row + 1][col][plane] =
                 noInterpolation ?
                 0 :
                 (outtuples[row][col][plane] +
@@ -189,13 +202,13 @@ calc_4(const struct pam * const pamP,
 
 
 static void
-calc_5(const struct pam * const pamP,
-       tuple **           const intuples,
-       tuple **           const outtuples,
-       unsigned int       const plane,
-       bool               const noInterpolation,
-       unsigned int       const xoffset,
-       unsigned int       const yoffset) {
+calc5(const struct pam * const pamP,
+      tuple **           const intuples,
+      tuple **           const outtuples,
+      unsigned int       const plane,
+      bool               const noInterpolation,
+      unsigned int       const xoffset,
+      unsigned int       const yoffset) {
 /*----------------------------------------------------------------------------
    . X .
    X . X
@@ -231,7 +244,7 @@ calc_5(const struct pam * const pamP,
 
 
 
-struct compAction {
+struct CompAction {
     unsigned int xoffset;
     unsigned int yoffset;
     void (*calc)(const struct pam * const pamP,
@@ -244,54 +257,68 @@ struct compAction {
 };
 
 
+struct BayerPattern {
+
+    struct CompAction compAction[3];
+        /* compAction[n] tells how to compute Plane 'n' */
+};
+
 
-static struct compAction const comp_1[3] = {
+
+static struct BayerPattern const bayer1 = {
 /*----------------------------------------------------------------------------
   G B G B
   R G R G
   G B G B
   R G R G
 -----------------------------------------------------------------------------*/
-
-    { 0, 1, calc_4 },
-    { 0, 1, calc_5 },
-    { 1, 0, calc_4 }
+    {  /* compAction */
+        { 0, 1, calc4 },
+        { 0, 1, calc5 },
+        { 1, 0, calc4 }
+    }
 };
 
-static struct compAction const comp_2[3] = {
+static struct BayerPattern const bayer2 = {
 /*----------------------------------------------------------------------------
   R G R G
   G B G B
   R G R G
   G B G B
 -----------------------------------------------------------------------------*/
-    { 0, 0, calc_4 },
-    { 0, 0, calc_5 },
-    { 1, 1, calc_4 }
+    {  /* compAction */
+        { 0, 0, calc4 },
+        { 0, 0, calc5 },
+        { 1, 1, calc4 }
+    }
 };
 
-static struct compAction const comp_3[3] = {
+static struct BayerPattern const bayer3 = {
 /*----------------------------------------------------------------------------
   B G B G
   G R G R
   B G B G
   G R G R
 -----------------------------------------------------------------------------*/
-    { 1, 1, calc_4 },
-    { 0, 0, calc_5 },
-    { 0, 0, calc_4 }
+    {  /* compAction */
+        { 1, 1, calc4 },
+        { 0, 0, calc5 },
+        { 0, 0, calc4 }
+    }
 };
 
-static struct compAction const comp_4[3] = {
+static struct BayerPattern const bayer4 = {
 /*----------------------------------------------------------------------------
   G R G R
   B G B G
   G R G R
   B G B G
 -----------------------------------------------------------------------------*/
-    { 1, 0, calc_4 },
-    { 0, 1, calc_5 },
-    { 0, 1, calc_4 }
+    {  /* compAction */
+        { 1, 0, calc4 },
+        { 0, 1, calc5 },
+        { 0, 1, calc4 }
+    }
 };
 
 
@@ -315,56 +342,83 @@ makeOutputPam(const struct pam * const inpamP,
 
 
 
-static const struct compAction *
-actionTableForType(enum bayerType const bayerType) {
+struct XyOffset {
+/*----------------------------------------------------------------------------
+   A two-dimensional offset within a matrix.
+-----------------------------------------------------------------------------*/
+    unsigned int row;
+    unsigned int col;
+};
+
+
+
+static const struct CompAction *
+actionTableForType(enum BayerType const bayerType) {
 
-    const struct compAction * retval;
+    const struct CompAction * retval;
 
     switch (bayerType) {
-    case BAYER1: retval = comp_1; break;
-    case BAYER2: retval = comp_2; break;
-    case BAYER3: retval = comp_3; break;
-    case BAYER4: retval = comp_4; break;
+    case BAYER1: retval = bayer1.compAction; break;
+    case BAYER2: retval = bayer2.compAction; break;
+    case BAYER3: retval = bayer3.compAction; break;
+    case BAYER4: retval = bayer4.compAction; break;
     }
     return retval;
 }
 
 
 
-int
+static void
+calcImage(struct pam *   const inpamP,
+          tuple **       const intuples,
+          struct pam *   const outpamP,
+          tuple **       const outtuples,
+          enum BayerType const bayerType,
+          bool           const wantNoInterpolate) {
+
+    const struct CompAction * const compActionTable =
+        actionTableForType(bayerType);
+
+    unsigned int plane;
+
+    clearTuples(outpamP, outtuples);
+
+    for (plane = 0; plane < 3; ++plane) {
+
+        struct CompAction const compAction = compActionTable[plane];
+
+        compAction.calc(inpamP, intuples, outtuples, plane,
+                        wantNoInterpolate,
+                        compAction.xoffset, compAction.yoffset);
+    }
+}
+
+
+            int
 main(int argc, const char **argv) {
 
-    struct cmdlineInfo cmdline;
+    struct CmdlineInfo cmdline;
     FILE * ifP;
     struct pam inpam;
     struct pam outpam;
     tuple ** intuples;
     tuple ** outtuples;
-    const struct compAction * compActionTable;
-    unsigned int plane;
 
     pm_proginit(&argc, argv);
-    
+
     parseCommandLine(argc, argv, &cmdline);
-    
+
     ifP = pm_openr(cmdline.inputFilespec);
-    
-    intuples = pnm_readpam(ifP, &inpam, PAM_STRUCT_SIZE(tuple_type));
 
-    compActionTable = actionTableForType(cmdline.bayerType);
+    intuples = pnm_readpam(ifP, &inpam, PAM_STRUCT_SIZE(tuple_type));
 
     makeOutputPam(&inpam, &outpam);
 
     outtuples = pnm_allocpamarray(&outpam);
-    clearTuples(&outpam, outtuples);
 
-    for (plane = 0; plane < 3; ++plane) {
-        struct compAction const compAction = compActionTable[plane];
+    calcImage(&inpam, intuples, &outpam, outtuples,cmdline.bayerType,
+              !!cmdline.nointerpolate);
 
-        compAction.calc(&inpam, intuples, outtuples, plane,
-                        cmdline.nointerpolate,
-                        compAction.xoffset, compAction.yoffset);
-    }
     pnm_writepam(&outpam, outtuples);
 
     pnm_freepamarray(outtuples, &outpam);
@@ -372,3 +426,6 @@ main(int argc, const char **argv) {
 
     return 0;
 }
+
+
+
diff --git a/other/pamdepth.c b/other/pamdepth.c
index 96613d07..46601864 100644
--- a/other/pamdepth.c
+++ b/other/pamdepth.c
@@ -17,7 +17,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,7 +30,7 @@ struct cmdlineInfo {
 
 static void
 parseCommandLine(int argc, const char ** argv,
-                 struct cmdlineInfo *cmdlineP) {
+                 struct CmdlineInfo * const cmdlineP) {
 /*----------------------------------------------------------------------------
    Note that the file spec strings we return are stored in the storage that
    was passed to us as the argv array.
@@ -45,7 +45,7 @@ parseCommandLine(int argc, const char ** argv,
     MALLOCARRAY_NOFAIL(option_def, 100);
 
     option_def_index = 0;   /* incremented by OPTENTRY */
-    OPTENT3(0, "verbose",  OPT_STRING, NULL, 
+    OPTENT3(0, "verbose",  OPT_STRING, NULL,
             &cmdlineP->verbose, 0);
 
     opt.opt_table = option_def;
@@ -64,8 +64,8 @@ parseCommandLine(int argc, const char ** argv,
             pm_error("New maxval must be at least 1.  You specified %d",
                      intval);
         else if (intval > PNM_OVERALLMAXVAL)
-            pm_error("newmaxval (%d) is too large.\n"
-                     "The maximum allowed by the PNM formats is %d.",
+            pm_error("newmaxval (%d) is too large.  "
+                     "The maximum allowed by the PNM formats is %u.",
                      intval, PNM_OVERALLMAXVAL);
         else
             cmdlineP->newMaxval = intval;
@@ -101,7 +101,7 @@ createSampleMap(sample   const oldMaxval,
 static void
 transformRaster(struct pam * const inpamP,
                 struct pam * const outpamP) {
-                
+
     tuple * tuplerow;
     unsigned int row;
     sample * sampleMap;  /* malloc'ed */
@@ -136,7 +136,7 @@ transformRaster(struct pam * const inpamP,
 int
 main(int argc, const char * argv[]) {
 
-    struct cmdlineInfo cmdline;
+    struct CmdlineInfo cmdline;
     FILE * ifP;
     struct pam inpam;
     struct pam outpam;
@@ -153,15 +153,15 @@ main(int argc, const char * argv[]) {
         pnm_readpaminit(ifP, &inpam, PAM_STRUCT_SIZE(tuple_type));
 
         outpam = inpam;  /* initial value */
-        
+
         outpam.file = stdout;
         outpam.maxval = cmdline.newMaxval;
-        
+
         if (PNM_FORMAT_TYPE(inpam.format) == PBM_TYPE) {
             outpam.format = PGM_TYPE;
         } else
             outpam.format = inpam.format;
-        
+
         if (streq(inpam.tuple_type, PAM_PBM_TUPLETYPE)) {
             pm_message("promoting from black and white to grayscale");
             strcpy(outpam.tuple_type, PAM_PGM_TUPLETYPE);
@@ -182,3 +182,6 @@ main(int argc, const char * argv[]) {
 
     return 0;
 }
+
+
+
diff --git a/other/pamendian.c b/other/pamendian.c
index c24d9629..8826b06f 100644
--- a/other/pamendian.c
+++ b/other/pamendian.c
@@ -1,20 +1,21 @@
-/******************************************************************************
+/*=============================================================================
                               pamendian
-*******************************************************************************
+===============================================================================
 
   Reverse the endianness of multi-byte samples in a Netpbm stream.
   I.e. convert between the true format and the little endian variation of
   it.
-******************************************************************************/
-  
+=============================================================================*/
+
 #include "pam.h"
 
 
 static sample
-reverseSample(sample const insample, unsigned int const bytesPerSample) {
+reverseSample(sample       const insample,
+              unsigned int const bytesPerSample) {
 /*----------------------------------------------------------------------------
   Return a sample whose value is the least significant
-  'bytes_per_sample' bytes, in reverse order.
+  'bytesPerSample' bytes, in reverse order.
 -----------------------------------------------------------------------------*/
     unsigned int bytePos;
     sample shiftedInsample;
@@ -30,14 +31,14 @@ reverseSample(sample const insample, unsigned int const bytesPerSample) {
 
 
 
-int main(int argc, char *argv[]) {
+int main(int argc, const char ** argv) {
 
     struct pam inpam, outpam;
     tuple * intuplerow;
     tuple * outtuplerow;
     unsigned int row;
 
-    pnm_init(&argc, argv);
+    pm_proginit(&argc, argv);
 
     pnm_readpaminit(stdin, &inpam, PAM_STRUCT_SIZE(tuple_type));
 
@@ -46,27 +47,39 @@ int main(int argc, char *argv[]) {
 
     pnm_writepaminit(&outpam);
 
-    intuplerow = pnm_allocpamrow(&inpam);      
-    outtuplerow = pnm_allocpamrow(&outpam);
-
+    /* We read the samples as if the maxval is 65535 so pnm_readpamrow, which
+       assumes big-endian samples, doesn't choke on a little-endian sample,
+       finding it to exceed the maxval.  (The pure way to do this would be not
+       to use libnetpbm row reading and writing facilities on little-endian
+       pseudo-Netpbm images, but this program isn't important enough to
+       justify that effort).
+    */
     inpam.maxval = 65535;
 
-    for (row = 0; row < inpam.height; row++) {
+    intuplerow  = pnm_allocpamrow(&inpam);
+    outtuplerow = pnm_allocpamrow(&outpam);
+
+    for (row = 0; row < inpam.height; ++row) {
         unsigned int col;
+
         pnm_readpamrow(&inpam, intuplerow);
-        for (col = 0; col < inpam.width; col++) {
+        for (col = 0; col < inpam.width; ++col) {
             unsigned int plane;
-            for (plane = 0; plane < inpam.depth; plane++) 
-                outtuplerow[col][plane] = 
-                    reverseSample(intuplerow[col][plane], 
+
+            for (plane = 0; plane < inpam.depth; ++plane) {
+                outtuplerow[col][plane] =
+                    reverseSample(intuplerow[col][plane],
                                   inpam.bytes_per_sample);
+            }
         }
         pnm_writepamrow(&outpam, outtuplerow);
     }
 
-    pnm_freepamrow(outtuplerow);        
-    pnm_freepamrow(intuplerow);        
+    pnm_freepamrow(outtuplerow);
+    pnm_freepamrow(intuplerow);
 
-    exit(0);
+    return 0;
 }
 
+
+
diff --git a/other/pamexec.c b/other/pamexec.c
index c3a1ee78..cbde81ac 100644
--- a/other/pamexec.c
+++ b/other/pamexec.c
@@ -15,8 +15,11 @@
 #define _BSD_SOURCE 1      /* Make sure strdup() is in string.h */
 #define _XOPEN_SOURCE 500  /* Make sure strdup() is in string.h */
 
+#include <stdbool.h>
 #include <string.h>
 #include <stdio.h>
+#include <signal.h>
+#include <setjmp.h>
 #include <errno.h>
 
 #include "pm_c_util.h"
@@ -25,7 +28,27 @@
 #include "mallocvar.h"
 #include "pam.h"
 
-struct cmdlineInfo {
+
+/* About SIGPIPE:
+
+   Unix has a strange function where by default, if you write into a pipe when
+   the reading side of the pipe has been closed, it generates a signal (of
+   class SIGPIPE), but if you tell the OS to ignore signals of class SIGPIPE,
+   then instead of generating that signal, the system call to write to the
+   pipe just fails.
+
+   Pamexec writes to a pipe when it feeds an image to the user's program's
+   Standard Input.  Should the user's program close its end of the pipe (such
+   as by exiting) before reading the whole image, that would, if we did
+   nothing to deal with it, cause Pamexec to receive a SIGPIPE signal, which
+   would make the OS terminate Pamexec.
+
+   We don't want that, so we tell the OS to ignore SIGPIPE signals, so that
+   instead our attempt to write to the pipe just fails and we fail with a
+   meaningful error message.
+*/
+
+struct CmdlineInfo {
     /* All the information the user supplied in the command line,
        in a form easy for the program to use.
     */
@@ -40,7 +63,7 @@ struct cmdlineInfo {
 
 static void
 parseCommandLine(int argc, const char ** argv,
-                 struct cmdlineInfo * const cmdlineP) {
+                 struct CmdlineInfo * const cmdlineP) {
 /*----------------------------------------------------------------------------
    Note that the pointers we place into *cmdlineP are sometimes to storage
    in the argv array.
@@ -59,13 +82,13 @@ parseCommandLine(int argc, const char ** argv,
     OPTENT3(0,   "check",   OPT_FLAG,   NULL,         &cmdlineP->check, 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 */
+    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 < 1) 
+    if (argc-1 < 1)
         pm_error("You must specify at least one argument: the shell command "
                  "to execute");
     else {
@@ -92,7 +115,7 @@ pipeOneImage(FILE * const infileP,
     struct pam inpam;
     struct pam outpam;
     enum pm_check_code checkRetval;
-    
+
     unsigned int row;
     tuple * tuplerow;
 
@@ -107,9 +130,25 @@ pipeOneImage(FILE * const infileP,
 
     tuplerow = pnm_allocpamrow(&inpam);
 
-    for (row = 0; row < inpam.height; ++row) {
-        pnm_readpamrow(&inpam, tuplerow);
-        pnm_writepamrow(&outpam, tuplerow);
+    {
+        jmp_buf jmpbuf;
+        int rc;
+        rc = setjmp(jmpbuf);
+        if (rc == 0) {
+            pm_setjmpbuf(&jmpbuf);
+
+            for (row = 0; row < inpam.height; ++row) {
+                pnm_readpamrow(&inpam, tuplerow);
+                pnm_writepamrow(&outpam, tuplerow);
+            }
+        } else {
+            pm_setjmpbuf(NULL);
+            pm_error("Failed to read image and pipe it to program's "
+                     "Standard Input.  If previous messages indicate "
+                     "a broken pipe error, that means the program closed "
+                     "its Standard Error (possibly by exiting) before "
+                     "it had read the entire image>");
+        }
     }
 
     pnm_freepamrow(tuplerow);
@@ -133,7 +172,7 @@ doOneImage(FILE *        const ifP,
     ofP = popen(command, "w");
 
     if (ofP == NULL)
-        pm_asprintf(errorP, 
+        pm_asprintf(errorP,
                     "Failed to start shell to run command '%s'.  "
                     "errno = %d (%s)",
                     command, errno, strerror(errno));
@@ -141,7 +180,7 @@ doOneImage(FILE *        const ifP,
         int rc;
 
         pipeOneImage(ifP, ofP);
-            
+
         rc = pclose(ofP);
 
         if (check && rc != 0)
@@ -157,7 +196,7 @@ doOneImage(FILE *        const ifP,
 int
 main(int argc, const char *argv[]) {
 
-    struct cmdlineInfo cmdline;
+    struct CmdlineInfo cmdline;
 
     FILE *       ifP;         /* Input file pointer */
     int          eof;         /* No more images in input */
@@ -169,12 +208,17 @@ main(int argc, const char *argv[]) {
     pm_proginit(&argc, argv);
 
     parseCommandLine(argc, argv, &cmdline);
-    
+
+    /* Make write to closed pipe fail rather than generate a signal.
+       See comments at top of program.
+    */
+    signal(SIGPIPE, SIG_IGN);
+
     ifP = pm_openr(cmdline.inputFileName);
 
-    for (eof = FALSE, imageSeq = 0; !eof; ++imageSeq) {
+    for (eof = false, imageSeq = 0; !eof; ++imageSeq) {
         const char * error;
-        
+
         doOneImage(ifP, cmdline.command, cmdline.check, &error);
 
         if (error) {
@@ -185,9 +229,8 @@ main(int argc, const char *argv[]) {
         pnm_nextimage(ifP, &eof);
     }
     pm_close(ifP);
-    
+
     return 0;
 }
 
 
-
diff --git a/other/pamfix.c b/other/pamfix.c
index e40f51f4..4f4eedfc 100644
--- a/other/pamfix.c
+++ b/other/pamfix.c
@@ -109,7 +109,7 @@ analyzeRaster(const struct pam * const pamP,
               bool               const verbose) {
 /*----------------------------------------------------------------------------
    Go through the raster at which the stream described by *tweakedPamP is
-   presently positioned and count how many rows can be successfully read
+   currently positioned and count how many rows can be successfully read
    (including validating the samples against pamP->maxval) and determine the
    highest sample value in those rows.
 
diff --git a/other/pamunlookup.c b/other/pamunlookup.c
index 77c2807b..defa7b1f 100644
--- a/other/pamunlookup.c
+++ b/other/pamunlookup.c
@@ -19,7 +19,6 @@
 #include "pm_c_util.h"
 #include "mallocvar.h"
 #include "shhopt.h"
-#include "pm_system.h"
 #include "nstring.h"
 #include "pam.h"
 #include "pammap.h"
@@ -49,13 +48,13 @@ parseCommandLine(int argc, const char ** const argv,
     optStruct3 opt;
 
     unsigned int option_def_index;
-    
+
     unsigned int lookupfileSpec;
 
     MALLOCARRAY_NOFAIL(option_def, 100);
 
     option_def_index = 0;   /* incremented by OPTENTRY */
-    OPTENT3(0, "lookupfile",     OPT_STRING, &cmdlineP->lookupfile,  
+    OPTENT3(0, "lookupfile",     OPT_STRING, &cmdlineP->lookupfile,
             &lookupfileSpec, 0);
 
     opt.opt_table = option_def;
@@ -74,17 +73,17 @@ parseCommandLine(int argc, const char ** const argv,
         cmdlineP->inputFileName = argv[1];
 
     free(option_def);
-}        
+}
 
 
 
 static void
-getLookup(const char * const lookupFileName, 
+getLookup(const char * const lookupFileName,
           tuple ***    const lookupP,
           struct pam * const lookuppamP) {
 /*----------------------------------------------------------------------------
    Get the lookup image (the one that maps integers to tuples, e.g. a
-   color index / color map / palette) from the file named 
+   color index / color map / palette) from the file named
    'lookupFileName'.
 
    Return the image as *lookupP and *lookuppamP.
@@ -95,14 +94,14 @@ getLookup(const char * const lookupFileName,
     tuple ** inputLookup;
 
     lookupfileP = pm_openr(lookupFileName);
-    inputLookup = pnm_readpam(lookupfileP, 
+    inputLookup = pnm_readpam(lookupfileP,
                               &inputLookuppam, PAM_STRUCT_SIZE(tuple_type));
 
     pm_close(lookupfileP);
-    
+
     if (inputLookuppam.height != 1)
         pm_error("The lookup table image must be one row.  "
-                 "Yours is %u rows.", 
+                 "Yours is %u rows.",
                  inputLookuppam.height);
 
     *lookupP = inputLookup;
@@ -129,7 +128,7 @@ makeReverseLookupHash(struct pam * const lookuppamP,
 
     for (col = 0; col < lookuppamP->width; ++col) {
         tuple const thisValue = lookup[0][col];
-        
+
         int found;
         int priorValue;
 
@@ -181,7 +180,7 @@ doUnlookup(struct pam * const inpamP,
 
     for (row = 0; row < inpamP->height; ++row) {
         unsigned int col;
-        
+
         pnm_readpamrow(inpamP, inrow);
 
         for (col = 0; col < inpamP->width; ++col) {
@@ -214,7 +213,7 @@ main(int argc, const char ** const argv) {
     tuple ** lookup;
 
     tuplehash lookupHash;
-    
+
     pm_proginit(&argc, argv);
 
     parseCommandLine(argc, argv, &cmdline);
@@ -230,7 +229,7 @@ main(int argc, const char ** const argv) {
                  "has depth %u.  They must be the same",
                  lookuppam.depth, inpam.depth);
     if (!streq(inpam.tuple_type, lookuppam.tuple_type))
-        pm_error("The lookup image has tupel type '%s', "
+        pm_error("The lookup image has tuple type '%s', "
                  "but the input image "
                  "has tuple type '%s'.  They must be the same",
                  lookuppam.tuple_type, inpam.tuple_type);
@@ -243,7 +242,7 @@ main(int argc, const char ** const argv) {
 
     pnm_destroytuplehash(lookupHash);
     pnm_freepamarray(lookup, &lookuppam);
-    
+
     return 0;
 }
 
diff --git a/other/pamx/Makefile b/other/pamx/Makefile
index e4892540..7fe7b9db 100644
--- a/other/pamx/Makefile
+++ b/other/pamx/Makefile
@@ -9,7 +9,10 @@ include $(BUILDDIR)/config.mk
 
 EXTERN_INCLUDE =
 
-ifeq ($(shell $(PKG_CONFIG) x11 --modversion --silence-errors),)
+TEST_PKGCONFIG_X11 := \
+  if $(PKG_CONFIG) x11 --exists; then echo exists; fi
+
+ifeq ($(shell $(TEST_PKGCONFIG_X11)),)
   # Pkg-config has never heard of X11, or doesn't even exist
 
   ifneq ($(X11LIB),NONE)
@@ -46,7 +49,8 @@ all: $(BINARIES)
 
 include $(SRCDIR)/common.mk
 
-ifeq ($(shell $(PKG_CONFIG) x11 --libs),)
+ifeq ($(shell $(TEST_PKGCONFIG_X11)),)
+  # Pkg-config has never heard of X11, or doesn't even exist
   X11_LIBOPTS = $(shell $(LIBOPT) $(LIBOPTR) $(X11LIB))
 else
   X11_LIBOPTS = $(shell $(PKG_CONFIG) x11 --libs)
diff --git a/other/pamx/window.c b/other/pamx/window.c
index d44008c9..61089ea5 100644
--- a/other/pamx/window.c
+++ b/other/pamx/window.c
@@ -1,8 +1,8 @@
 /*
    Functions to allocate and deallocate structures and structure data
- 
+
    By Jim Frost 1989.10.03, Bryan Henderson 2006.03.25.
- 
+
    See COPYRIGHT file for copyright information.
 */
 
@@ -80,10 +80,10 @@ static void
 setDeleteWindow(viewer * const viewerP) {
 
     Atom const protoAtom = XInternAtom(viewerP->dispP, "WM_PROTOCOLS", False);
-    
+
     viewerP->deleteAtom = XInternAtom(viewerP->dispP,
                                       "WM_DELETE_WINDOW", False);
-    
+
     if ((protoAtom != None) && (viewerP->deleteAtom != None))
         XChangeProperty(viewerP->dispP, viewerP->viewportWin,
                         protoAtom, XA_ATOM, 32, PropModeReplace,
@@ -139,7 +139,7 @@ createViewer(viewer **     const viewerPP,
              bool          const fullscreen) {
 
     viewer * viewerP;
-    
+
     XSetWindowAttributes  swa_view;
 
     MALLOCVAR_NOFAIL(viewerP);
@@ -164,7 +164,7 @@ createViewer(viewer **     const viewerPP,
         ButtonPressMask | Button1MotionMask | KeyPressMask |
         StructureNotifyMask | EnterWindowMask | LeaveWindowMask;
     swa_view.save_under       = FALSE;
-    
+
     viewerP->viewportWin =
         XCreateWindow(dispP, RootWindow(dispP, scrn),
                       viewerP->xpos, viewerP->ypos,
@@ -176,7 +176,7 @@ createViewer(viewer **     const viewerPP,
                       &swa_view);
 
     setXloadimageClassHint(viewerP->dispP, viewerP->viewportWin);
-    
+
     setDeleteWindow(viewerP);
 
     viewerP->blank = TRUE;
@@ -192,7 +192,7 @@ determineRepaintStrategy(viewer  *    const viewerP,
                          bool         const verbose,
                          XImageInfo * const ximageinfoP,
                          Pixmap *     const pixmapP) {
-                        
+
     /* Decide how we're going to handle repaints.  We have three modes:
        use backing-store, use background pixmap, and use exposures.
        If the server allows backing-store, we enable it and use it.
@@ -200,7 +200,7 @@ determineRepaintStrategy(viewer  *    const viewerP,
        server does not have backing-store, we try to send the image to
        a pixmap and use that as backing-store.  If that fails, we use
        exposures to blit the image (which is ugly but it works).
-       
+
        'use_pixmap' forces background pixmap mode, which may improve
        performance.
     */
@@ -303,7 +303,7 @@ createImageWindow(viewer *      const viewerP,
 
 static void
 destroyImageWindow(viewer * const viewerP) {
-    
+
     if (viewerP->imageWin) {
         if (viewerP->imageColormap &&
             (viewerP->imageColormap != DefaultColormap(viewerP->dispP, viewerP->scrn)))
@@ -311,7 +311,7 @@ destroyImageWindow(viewer * const viewerP) {
         XDestroyWindow(viewerP->dispP, viewerP->imageWin);
     }
 }
-                       
+
 
 
 static void
@@ -348,12 +348,32 @@ placeImage(viewer * const viewerP,
            int      const height,
            int *    const rxP,     /* input and output */
            int *    const ryP) {   /* input and output */
+/*----------------------------------------------------------------------------
+   Move the X window of *viewerP on the right place on the display to contain
+   an image whose dimensions are 'width' by 'height'.
+
+   As input *rxP and *ryP are the current location of the window (offset from
+   the top left corner of the display of the top left corner of the window).
+   We update them with the new location.
+
+   If the image is narrower than the display, center it.
+
+   Otherwise, either left-justify it or right-justify it depending upon the
+   current position of the window:
+
+      If in the current position, the right edge of the image is within the
+      display, right-justify.
 
+      Otherwise left-justify.
+
+
+   Do the analogous thing for vertical placement.
+-----------------------------------------------------------------------------*/
     int pixx, pixy;
-    
+
     pixx = *rxP;
     pixy = *ryP;
-    
+
     if (viewerP->width > width)
         pixx = (viewerP->width - width) / 2;
     else {
@@ -434,10 +454,10 @@ destroyViewer(viewer * const viewerP) {
     if (viewerP->imageWin)
         XDestroyWindow(viewerP->dispP, viewerP->imageWin);
     viewerP->imageWin = 0;
-    
+
     if (viewerP->viewportWin)
         XDestroyWindow(viewerP->dispP, viewerP->viewportWin);
-    
+
     viewerP->viewportWin = 0;
 
     XFreeCursor(viewerP->dispP, viewerP->cursor);
@@ -459,12 +479,12 @@ setViewportColormap(viewer *  const viewerP,
 
     if (cmap_atom == None)
         cmap_atom = XInternAtom(viewerP->dispP, "WM_COLORMAP_WINDOWS", False);
-    
+
     /* If the visual we're using is the same as the default visual (used by
        the viewport window) then we can set the viewport window to use the
        image's colormap.  This keeps most window managers happy.
     */
-    
+
     if (visualP == DefaultVisual(viewerP->dispP, viewerP->scrn)) {
         swa.colormap = viewerP->imageColormap;
         XChangeWindowAttributes(viewerP->dispP, viewerP->viewportWin,
@@ -507,11 +527,11 @@ iconName(const char * const s) {
         t = strchr(buf, ' ');
         if (t)
             *t = '\0';
-    
+
         /* Strip off leading path.  if you don't use unix-style paths,
            You might want to change this.
         */
-    
+
         t= strrchr(buf, '/');
         if (t) {
             char * p;
@@ -557,7 +577,7 @@ visualClassFromName(const char * const name) {
     unsigned int a;
     int class;
     bool found;
-    
+
     for (a = 0, found = FALSE; VisualClassName[a].name; ++a) {
         if (strcaseeq(VisualClassName[a].name, name)) {
             /* Check for uniqueness.  We special-case StaticGray
@@ -674,10 +694,10 @@ bestVisual(Display *      const disp,
         pm_message("bestVisual: didn't find any depths?!?");
         depth = DefaultDepth(disp, scrn);
     }
-    
+
     /* given this depth, find the best possible visual
      */
-    
+
     default_visualP = DefaultVisual(disp, scrn);
     switch (imageP->type) {
     case ITRUE: {
@@ -697,7 +717,7 @@ bestVisual(Display *      const disp,
                 visualP = bestVisualOfClassAndDepth(disp, scrn, TrueColor,
                                                     depth);
         }
-        
+
         if (!visualP || ((depth <= 8) &&
                          bestVisualOfClassAndDepth(disp, scrn, PseudoColor,
                                                    depth)))
@@ -711,7 +731,7 @@ bestVisual(Display *      const disp,
         if (!visualP)
             visualP = bestVisualOfClassAndDepth(disp, scrn, StaticGray, depth);
     } break;
-        
+
     case IRGB: {
         /* if it's an RGB image, we want PseudoColor if we can get it */
 
@@ -814,7 +834,7 @@ getImageDispDimensions(viewer *       const viewerP,
         *widthP  = viewerP->width;
         *heightP = viewerP->height;
     } else {
-        unsigned int const displayWidth = 
+        unsigned int const displayWidth =
             DisplayWidth(viewerP->dispP, viewerP->scrn);
         unsigned int const displayHeight =
             DisplayHeight(viewerP->dispP, viewerP->scrn);
@@ -878,9 +898,9 @@ getVisualAndDepth(Image *        const imageP,
 static void
 setNormalSizeHints(viewer *     const viewerP,
                    Image *      const imageP) {
-    
+
     XSizeHints sh;
-    
+
     sh.width  = viewerP->width;
     sh.height = viewerP->height;
     if (viewerP->fullscreen) {
@@ -984,7 +1004,7 @@ run(viewer *     const viewerP,
             KeySym ks;
             XComposeStatus status;
             Status rc;
-            
+
             rc = XLookupString(&event.key, buf, 128, &ks, &status);
             if (rc == 1) {
                 char const ret = buf[0];
@@ -1093,7 +1113,7 @@ static void
 resizeViewer(viewer *     const viewerP,
              unsigned int const newWidth,
              unsigned int const newHeight) {
-    
+
     if (viewerP->width != newWidth || viewerP->height != newHeight) {
         XResizeWindow(viewerP->dispP, viewerP->viewportWin,
                       newWidth, newHeight);
@@ -1116,7 +1136,7 @@ displayInViewer(viewer *       const viewerP,
                 const char *   const title,
                 bool           const verbose,
                 int *          const retvalP) {
-    
+
     XImageInfo *     ximageInfoP;
     Visual *         visual;
     unsigned int     depth;
@@ -1134,7 +1154,7 @@ displayInViewer(viewer *       const viewerP,
     getVisualAndDepth(imageP, viewerP->dispP, viewerP->scrn,
                       fit, visualSpec, visualClass,
                       &visual, &depth);
-    
+
     if (verbose && (visual != DefaultVisual(viewerP->dispP, viewerP->scrn)))
         pm_message("Using %s visual", nameOfVisualClass(visual->class));
 
@@ -1211,3 +1231,6 @@ displayInViewer(viewer *       const viewerP,
 
     *retvalP = retvalueFromExitReason(exitReason);
 }
+
+
+
diff --git a/other/pnmcolormap.c b/other/pnmcolormap.c
index f687f037..fbe85d4e 100644
--- a/other/pnmcolormap.c
+++ b/other/pnmcolormap.c
@@ -1,7 +1,6 @@
-/******************************************************************************
-                               pnmcolormap.c
-*******************************************************************************
-
+/*=============================================================================
+                               pnmcolormap
+===============================================================================
   Create a colormap file (a PPM image containing one pixel of each of a set
   of colors).  Base the set of colors on an input image.
 
@@ -20,9 +19,8 @@
   copyright notice and this permission notice appear in supporting
   documentation.  This software is provided "as is" without express or
   implied warranty.
-
-******************************************************************************/
-
+=============================================================================*/
+#include <assert.h>
 #include <math.h>
 
 #include "pm_config.h"
@@ -33,39 +31,73 @@
 #include "pam.h"
 #include "pammap.h"
 
-enum methodForLargest {LARGE_NORM, LARGE_LUM};
+enum MethodForLargest {LARGE_NORM, LARGE_LUM};
 
-enum methodForRep {REP_CENTER_BOX, REP_AVERAGE_COLORS, REP_AVERAGE_PIXELS};
+enum MethodForRep {REP_CENTER_BOX, REP_AVERAGE_COLORS, REP_AVERAGE_PIXELS};
+
+enum MethodForSplit {SPLIT_MAX_PIXELS, SPLIT_MAX_COLORS, SPLIT_MAX_SPREAD};
+
+struct Box {
+/*----------------------------------------------------------------------------
+   A box contains an extent of a color frequency table, i.e. the colors
+   with some consecutive index values in the color frequency table.
+-----------------------------------------------------------------------------*/
+    unsigned int startIndex;
+        /* First index in the extent */
+    unsigned int colorCt;
+        /* Size of the extent (Number of colors in it -- at least 1) */
+    unsigned int sum;
+        /* Number of pixels of all colors in the extent */
+    unsigned int maxdim;
+        /* Which dimension has the largest spread.  RGB plane number. */
+        /* Meaningless if box contains only 1 color */
+    sample       spread;
+        /* spread in dimension 'maxdim' */
+        /* Meaningless if box contains only 1 color */
+};
 
-typedef struct box* boxVector;
-struct box {
-    int ind;
-    int colors;
-    int sum;
+struct BoxVector {
+    tupletable2 colorFreqTable;
+        /* The colors and their frequencies (number of pixels in the image of
+           that color), ordered into consecutive boxes, as defined by 'box'.
+        */
+    unsigned int colorDepth;
+        /* Number of planes in the tuples of 'colorFreqTable' */
+    struct Box * box;  /* malloc'ed array */
+        /* An array of boxes that contain consecutive extents of
+           'colorFreqTable'.  The list covers the entire table.
+        */
+    unsigned int boxCt;
+        /* Number of boxes in the above list */
+    unsigned int capacity;
+        /* Number of boxes the array is capable of containing */
 };
 
-struct cmdlineInfo {
+struct CmdlineInfo {
     /* All the information the user supplied in the command line,
        in a form easy for the program to use.
     */
-    const char *inputFilespec;  /* Filespec of input file */
+    const char * inputFileNm;  /* Name of input file */
     unsigned int allcolors;  /* boolean: select all colors from the input */
     unsigned int newcolors;
         /* Number of colors argument; meaningless if allcolors true */
-    enum methodForLargest methodForLargest;
+    enum MethodForLargest methodForLargest;
         /* -spreadintensity/-spreadluminosity options */
-    enum methodForRep methodForRep;
+    enum MethodForRep methodForRep;
         /* -center/-meancolor/-meanpixel options */
+    enum MethodForSplit methodForSplit;
+        /* -splitpixelct/-splitcolorct/-splitspread options */
     unsigned int sort;
     unsigned int square;
     unsigned int verbose;
+    unsigned int debug;
 };
 
 
 
 static void
-parseCommandLine (int argc, char ** argv,
-                  struct cmdlineInfo *cmdlineP) {
+parseCommandLine (int argc, const char ** argv,
+                  struct CmdlineInfo * const cmdlineP) {
 /*----------------------------------------------------------------------------
    parse program command line described in Unix standard form by argc
    and argv.  Return the information in the options as *cmdlineP.
@@ -85,6 +117,7 @@ parseCommandLine (int argc, char ** argv,
 
     unsigned int spreadbrightness, spreadluminosity;
     unsigned int center, meancolor, meanpixel;
+    unsigned int splitpixelct, splitcolorct, splitspread;
 
     MALLOCARRAY_NOFAIL(option_def, 100);
 
@@ -99,18 +132,26 @@ parseCommandLine (int argc, char ** argv,
             NULL,                       &meancolor,        0);
     OPTENT3(0,   "meanpixel",        OPT_FLAG,
             NULL,                       &meanpixel,        0);
+    OPTENT3(0,   "splitpixelct",     OPT_FLAG,
+            NULL,                       &splitpixelct,     0);
+    OPTENT3(0,   "splitcolorct",     OPT_FLAG,
+            NULL,                       &splitcolorct,     0);
+    OPTENT3(0,   "splitspread",      OPT_FLAG,
+            NULL,                       &splitspread,      0);
     OPTENT3(0, "sort",     OPT_FLAG,   NULL,
             &cmdlineP->sort,       0 );
     OPTENT3(0, "square",   OPT_FLAG,   NULL,
-            &cmdlineP->square,       0 );
+            &cmdlineP->square,     0 );
     OPTENT3(0, "verbose",   OPT_FLAG,   NULL,
-            &cmdlineP->verbose,       0 );
+            &cmdlineP->verbose,    0 );
+    OPTENT3(0, "debug",     OPT_FLAG,   NULL,
+            &cmdlineP->debug,      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, argv, opt, sizeof(opt), 0 );
+    pm_optParseOptions3( &argc, (char **)argv, opt, sizeof(opt), 0 );
         /* Uses and sets argc, argv, and some of *cmdline_p and others. */
 
 
@@ -132,15 +173,24 @@ parseCommandLine (int argc, char ** argv,
     else
         cmdlineP->methodForRep = REP_CENTER_BOX;
 
+    if (splitpixelct)
+        cmdlineP->methodForSplit = SPLIT_MAX_PIXELS;
+    else if (splitcolorct)
+        cmdlineP->methodForSplit = SPLIT_MAX_COLORS;
+    else if (splitspread)
+        cmdlineP->methodForSplit = SPLIT_MAX_SPREAD;
+    else
+        cmdlineP->methodForSplit = SPLIT_MAX_PIXELS;
+
     if (argc-1 > 2)
         pm_error("Program takes at most two arguments: number of colors "
                  "and input file specification.  "
                  "You specified %d arguments.", argc-1);
     else {
         if (argc-1 < 2)
-            cmdlineP->inputFilespec = "-";
+            cmdlineP->inputFileNm = "-";
         else
-            cmdlineP->inputFilespec = argv[2];
+            cmdlineP->inputFileNm = argv[2];
 
         if (argc-1 < 1)
             pm_error("You must specify the number of colors in the "
@@ -185,8 +235,13 @@ compareplane(const void * const arg1,
     const struct tupleint * const * const comparandPP  = arg1;
     const struct tupleint * const * const comparatorPP = arg2;
 
-    return (*comparandPP)->tuple[compareplanePlane] -
-        (*comparatorPP)->tuple[compareplanePlane];
+    sample const comparandSample  = (*comparandPP) ->tuple[compareplanePlane];
+    sample const comparatorSample = (*comparatorPP)->tuple[compareplanePlane];
+
+    return
+        comparandSample < comparatorSample ? -1 :
+        comparandSample > comparatorSample ? 1 :
+        0;
 }
 
 
@@ -196,84 +251,105 @@ static qsort_comparison_fn sumcompare;
 #endif
 
 static int
-sumcompare(const void * const b1, const void * const b2) {
-    return(((boxVector)b2)->sum - ((boxVector)b1)->sum);
+sumcompare(const void * const arg1,
+           const void * const arg2) {
+
+    struct Box * const comparandP  = (struct Box *)arg1;
+    struct Box * const comparatorP = (struct Box *)arg2;
+
+    return
+        comparatorP->sum < comparandP->sum ? -1 :
+        comparatorP->sum > comparandP->sum ? 1 :
+        0;
 }
 
 
 
+#ifndef LITERAL_FN_DEF_MATCH
+static qsort_comparison_fn colcompare;
+#endif
 
-/*
-** Here is the fun part, the median-cut colormap generator.  This is based
-** on Paul Heckbert's paper "Color Image Quantization for Frame Buffer
-** Display", SIGGRAPH '82 Proceedings, page 297.
-*/
+static int
+colcompare(const void * const arg1,
+           const void * const arg2) {
 
-static tupletable2
-newColorMap(unsigned int const newcolors,
-            unsigned int const depth) {
+    struct Box * const comparandP  = (struct Box *)arg1;
+    struct Box * const comparatorP = (struct Box *)arg2;
 
-    tupletable2 colormap;
-    unsigned int i;
-    struct pam pam;
+    return
+        comparatorP->colorCt < comparandP->colorCt ? -1 :
+        comparatorP->colorCt > comparandP->colorCt ? 1 :
+        0;
+}
 
-    pam.depth = depth;
 
-    colormap.table = pnm_alloctupletable(&pam, newcolors);
 
-    for (i = 0; i < newcolors; ++i) {
-        unsigned int plane;
-        for (plane = 0; plane < depth; ++plane)
-            colormap.table[i]->tuple[plane] = 0;
-    }
-    colormap.size = newcolors;
+#ifndef LITERAL_FN_DEF_MATCH
+static qsort_comparison_fn spreadcompare;
+#endif
 
-    return colormap;
+static int
+spreadcompare(const void * const arg1,
+              const void * const arg2) {
+
+    struct Box * const comparandP  = (struct Box *)arg1;
+    struct Box * const comparatorP = (struct Box *)arg2;
+
+    return
+        comparatorP->spread < comparandP->spread ? -1 :
+        comparatorP->spread > comparandP->spread ? 1 :
+        0;
 }
 
 
 
-static boxVector
-newBoxVector(int const colors, int const sum, int const newcolors) {
+static void
+sortBoxes(struct BoxVector *  const boxVectorP,
+          enum MethodForSplit const methodForSplit) {
 
-    boxVector bv;
-    MALLOCARRAY(bv, newcolors);
-    if (bv == NULL)
-        pm_error("out of memory allocating box vector table");
+    qsort_comparison_fn * comparisonFn;
 
-    /* Set up the initial box. */
-    bv[0].ind = 0;
-    bv[0].colors = colors;
-    bv[0].sum = sum;
+    switch (methodForSplit){
+    case SPLIT_MAX_PIXELS: comparisonFn = &sumcompare;    break;
+    case SPLIT_MAX_COLORS: comparisonFn = &colcompare;    break;
+    case SPLIT_MAX_SPREAD: comparisonFn = &spreadcompare; break;
+    }
 
-    return bv;
+    qsort((char*) &boxVectorP->box[0], boxVectorP->boxCt, sizeof(struct Box),
+          comparisonFn);
 }
 
 
 
+/*
+** Here is the fun part, the median-cut colormap generator.  This is based
+** on Paul Heckbert's paper "Color Image Quantization for Frame Buffer
+** Display", SIGGRAPH '82 Proceedings, page 297.
+*/
+
 static void
-findBoxBoundaries(tupletable2  const colorfreqtable,
+findBoxBoundaries(tupletable2  const colorFreqTable,
                   unsigned int const depth,
                   unsigned int const boxStart,
                   unsigned int const boxSize,
                   sample             minval[],
                   sample             maxval[]) {
 /*----------------------------------------------------------------------------
-  Go through the box finding the minimum and maximum of each
-  component - the boundaries of the box.
+  Go through the box finding the minimum and maximum of each component - the
+  boundaries of the box.
 -----------------------------------------------------------------------------*/
     unsigned int plane;
     unsigned int i;
 
     for (plane = 0; plane < depth; ++plane) {
-        minval[plane] = colorfreqtable.table[boxStart]->tuple[plane];
+        minval[plane] = colorFreqTable.table[boxStart]->tuple[plane];
         maxval[plane] = minval[plane];
     }
 
     for (i = 1; i < boxSize; ++i) {
         unsigned int plane;
         for (plane = 0; plane < depth; ++plane) {
-            sample const v = colorfreqtable.table[boxStart + i]->tuple[plane];
+            sample const v = colorFreqTable.table[boxStart + i]->tuple[plane];
             if (v < minval[plane]) minval[plane] = v;
             if (v > maxval[plane]) maxval[plane] = v;
         }
@@ -282,57 +358,161 @@ findBoxBoundaries(tupletable2  const colorfreqtable,
 
 
 
-static unsigned int
-largestByNorm(sample minval[], sample maxval[], unsigned int const depth) {
-
-    unsigned int largestDimension;
+static void
+findPlaneWithLargestSpreadByNorm(sample         const minval[],
+                                 sample         const maxval[],
+                                 unsigned int   const depth,
+                                 unsigned int * const planeP,
+                                 sample *       const spreadP) {
+
+    unsigned int planeWithLargest;
+    sample       largestSpreadSoFar;
     unsigned int plane;
 
-    sample largestSpreadSoFar = 0;
-    largestDimension = 0;
-    for (plane = 0; plane < depth; ++plane) {
+    for (plane = 0, largestSpreadSoFar = 0; plane < depth; ++plane) {
+
         sample const spread = maxval[plane]-minval[plane];
         if (spread > largestSpreadSoFar) {
-            largestDimension = plane;
             largestSpreadSoFar = spread;
+            planeWithLargest   = plane;
         }
     }
-    return largestDimension;
+    *planeP  = planeWithLargest;
+    *spreadP = largestSpreadSoFar;
 }
 
 
 
-static unsigned int
-largestByLuminosity(sample minval[], sample maxval[],
-                    unsigned int const depth) {
+static void
+findPlaneWithLargestSpreadByLuminosity(sample         const minval[],
+                                       sample         const maxval[],
+                                       unsigned int   const depth,
+                                       unsigned int * const planeP,
+                                       sample *       const spreadP) {
 /*----------------------------------------------------------------------------
    This subroutine presumes that the tuple type is either
    BLACKANDWHITE, GRAYSCALE, or RGB (which implies pamP->depth is 1 or 3).
    To save time, we don't actually check it.
 -----------------------------------------------------------------------------*/
-    unsigned int retval;
-
-    if (depth == 1)
-        retval = 0;
-    else {
+    if (depth == 1){
+        *planeP  = 0;
+        *spreadP = 0;
+    } else {
         /* An RGB tuple */
-        unsigned int largestDimension;
+        unsigned int planeWithLargest;
+        sample       largestSpreadSoFar;
         unsigned int plane;
-        double largestSpreadSoFar;
 
-        largestSpreadSoFar = 0.0;
+        assert(depth >= 3);
 
-        for (plane = 0; plane < 3; ++plane) {
+        for (plane = 0, largestSpreadSoFar = 0; plane < 3; ++plane) {
             double const spread =
                 pnm_lumin_factor[plane] * (maxval[plane]-minval[plane]);
             if (spread > largestSpreadSoFar) {
-                largestDimension = plane;
                 largestSpreadSoFar = spread;
+                planeWithLargest   = plane;
             }
         }
-        retval = largestDimension;
+        *planeP  = planeWithLargest;
+        *spreadP = largestSpreadSoFar;
+    }
+}
+
+
+
+static void
+computeBoxSpread(const struct Box *    const boxP,
+                 tupletable2           const colorFreqTable,
+                 unsigned int          const depth,
+                 enum MethodForLargest const methodForLargest,
+                 unsigned int *        const planeWithLargestP,
+                 sample *              const spreadP
+                 ) {
+/*----------------------------------------------------------------------------
+  Find the spread in the dimension in which it is greatest.
+
+  Return as *planeWithLargestP the number of that plane and as *spreadP the
+  spread in that plane.
+-----------------------------------------------------------------------------*/
+    sample * minval;  /* malloc'ed array */
+    sample * maxval;  /* malloc'ed array */
+
+    MALLOCARRAY_NOFAIL(minval, depth);
+    MALLOCARRAY_NOFAIL(maxval, depth);
+
+    findBoxBoundaries(colorFreqTable, depth, boxP->startIndex, boxP->colorCt,
+                      minval, maxval);
+
+    switch (methodForLargest) {
+    case LARGE_NORM:
+        findPlaneWithLargestSpreadByNorm(minval, maxval, depth,
+                                         planeWithLargestP, spreadP);
+        break;
+    case LARGE_LUM:
+        findPlaneWithLargestSpreadByLuminosity(minval, maxval, depth,
+                                               planeWithLargestP, spreadP);
+        break;
     }
-    return retval;
+    free(minval); free(maxval);
+}
+
+
+
+static unsigned int
+freqTotal(tupletable2 const freqTable) {
+
+    unsigned int i;
+    unsigned int sum;
+
+    for (i = 0, sum = 0; i < freqTable.size; ++i)
+        sum += freqTable.table[i]->value;
+
+    return sum;
+}
+
+
+
+static struct BoxVector
+newBoxVector(tupletable2           const colorFreqTable,
+             unsigned int          const capacity,
+             unsigned int          const depth,
+             enum MethodForLargest const methodForLargest) {
+
+    unsigned int const colorCt = colorFreqTable.size;
+    unsigned int const sum     = freqTotal(colorFreqTable);
+
+    struct BoxVector boxVector;
+
+    boxVector.colorFreqTable = colorFreqTable;
+    boxVector.colorDepth     = depth;
+
+    MALLOCARRAY(boxVector.box, capacity);
+
+    if (!boxVector.box)
+        pm_error("out of memory allocating box vector table");
+
+    /* Set up the initial box. */
+    boxVector.box[0].startIndex = 0;
+    boxVector.box[0].colorCt    = colorCt;
+    boxVector.box[0].sum        = sum;
+
+    computeBoxSpread(&boxVector.box[0], colorFreqTable, depth,
+                     methodForLargest,
+                     &boxVector.box[0].maxdim,
+                     &boxVector.box[0].spread);
+
+    boxVector.boxCt    = 1;
+    boxVector.capacity = capacity;
+
+    return boxVector;
+}
+
+
+
+static void
+destroyBoxVector(struct BoxVector const boxVector) {
+
+    free(boxVector.box);
 }
 
 
@@ -340,7 +520,7 @@ largestByLuminosity(sample minval[], sample maxval[],
 static void
 centerBox(int          const boxStart,
           int          const boxSize,
-          tupletable2  const colorfreqtable,
+          tupletable2  const colorFreqTable,
           unsigned int const depth,
           tuple        const newTuple) {
 
@@ -350,10 +530,10 @@ centerBox(int          const boxStart,
         int minval, maxval;
         unsigned int i;
 
-        minval = maxval = colorfreqtable.table[boxStart]->tuple[plane];
+        minval = maxval = colorFreqTable.table[boxStart]->tuple[plane];
 
         for (i = 1; i < boxSize; ++i) {
-            int const v = colorfreqtable.table[boxStart + i]->tuple[plane];
+            int const v = colorFreqTable.table[boxStart + i]->tuple[plane];
             minval = MIN( minval, v);
             maxval = MAX( maxval, v);
         }
@@ -363,10 +543,34 @@ centerBox(int          const boxStart,
 
 
 
+static tupletable2
+newColorMap(unsigned int const colorCt,
+            unsigned int const depth) {
+
+    tupletable2 colormap;
+    unsigned int i;
+    struct pam pam;
+
+    pam.depth = depth;
+
+    colormap.table = pnm_alloctupletable(&pam, colorCt);
+
+    for (i = 0; i < colorCt; ++i) {
+        unsigned int plane;
+        for (plane = 0; plane < depth; ++plane)
+            colormap.table[i]->tuple[plane] = 0;
+    }
+    colormap.size = colorCt;
+
+    return colormap;
+}
+
+
+
 static void
 averageColors(int          const boxStart,
               int          const boxSize,
-              tupletable2  const colorfreqtable,
+              tupletable2  const colorFreqTable,
               unsigned int const depth,
               tuple        const newTuple) {
 
@@ -379,7 +583,7 @@ averageColors(int          const boxStart,
         sum = 0;
 
         for (i = 0; i < boxSize; ++i)
-            sum += colorfreqtable.table[boxStart+i]->tuple[plane];
+            sum += colorFreqTable.table[boxStart+i]->tuple[plane];
 
         newTuple[plane] = ROUNDDIV(sum, boxSize);
     }
@@ -390,7 +594,7 @@ averageColors(int          const boxStart,
 static void
 averagePixels(int          const boxStart,
               int          const boxSize,
-              tupletable2  const colorfreqtable,
+              tupletable2  const colorFreqTable,
               unsigned int const depth,
               tuple        const newTuple) {
 
@@ -402,7 +606,7 @@ averagePixels(int          const boxStart,
     /* Count the tuples in question */
     n = 0;  /* initial value */
     for (i = 0; i < boxSize; ++i)
-        n += colorfreqtable.table[boxStart + i]->value;
+        n += colorFreqTable.table[boxStart + i]->value;
 
 
     for (plane = 0; plane < depth; ++plane) {
@@ -412,8 +616,8 @@ averagePixels(int          const boxStart,
         sum = 0;
 
         for (i = 0; i < boxSize; ++i)
-            sum += colorfreqtable.table[boxStart+i]->tuple[plane]
-                * colorfreqtable.table[boxStart+i]->value;
+            sum += colorFreqTable.table[boxStart+i]->tuple[plane]
+                * colorFreqTable.table[boxStart+i]->value;
 
         newTuple[plane] = ROUNDDIV(sum, n);
     }
@@ -422,12 +626,9 @@ averagePixels(int          const boxStart,
 
 
 static tupletable2
-colormapFromBv(unsigned int      const newcolors,
-               boxVector         const bv,
-               unsigned int      const boxes,
-               tupletable2       const colorfreqtable,
-               unsigned int      const depth,
-               enum methodForRep const methodForRep) {
+colormapFromBv(unsigned int      const colorCt,
+               struct BoxVector  const boxVector,
+               enum MethodForRep const methodForRep) {
     /*
     ** Ok, we've got enough boxes.  Now choose a representative color for
     ** each box.  There are a number of possible ways to make this choice.
@@ -437,23 +638,29 @@ colormapFromBv(unsigned int      const newcolors,
     ** method is to average all the pixels in the box.
     */
     tupletable2 colormap;
-    unsigned int bi;
+    unsigned int boxIdx;
 
-    colormap = newColorMap(newcolors, depth);
+    colormap = newColorMap(colorCt, boxVector.colorDepth);
 
-    for (bi = 0; bi < boxes; ++bi) {
+    for (boxIdx = 0; boxIdx < boxVector.boxCt; ++boxIdx) {
         switch (methodForRep) {
         case REP_CENTER_BOX:
-            centerBox(bv[bi].ind, bv[bi].colors, colorfreqtable, depth,
-                      colormap.table[bi]->tuple);
+            centerBox(boxVector.box[boxIdx].startIndex,
+                      boxVector.box[boxIdx].colorCt,
+                      boxVector.colorFreqTable, boxVector.colorDepth,
+                      colormap.table[boxIdx]->tuple);
             break;
         case REP_AVERAGE_COLORS:
-            averageColors(bv[bi].ind, bv[bi].colors, colorfreqtable, depth,
-                          colormap.table[bi]->tuple);
+            averageColors(boxVector.box[boxIdx].startIndex,
+                          boxVector.box[boxIdx].colorCt,
+                          boxVector.colorFreqTable, boxVector.colorDepth,
+                          colormap.table[boxIdx]->tuple);
             break;
         case REP_AVERAGE_PIXELS:
-            averagePixels(bv[bi].ind, bv[bi].colors, colorfreqtable, depth,
-                          colormap.table[bi]->tuple);
+            averagePixels(boxVector.box[boxIdx].startIndex,
+                          boxVector.box[boxIdx].colorCt,
+                          boxVector.colorFreqTable, boxVector.colorDepth,
+                          colormap.table[boxIdx]->tuple);
             break;
         default:
             pm_error("Internal error: invalid value of methodForRep: %d",
@@ -465,153 +672,170 @@ colormapFromBv(unsigned int      const newcolors,
 
 
 
-static unsigned int
-freqTotal(tupletable2 const freqtable) {
-
-    unsigned int i;
-    unsigned int sum;
-
-    sum = 0;
-
-    for (i = 0; i < freqtable.size; ++i)
-        sum += freqtable.table[i]->value;
-
-    return sum;
-}
-
-
 static void
-splitBox(boxVector             const bv,
-         unsigned int *        const boxesP,
-         unsigned int          const bi,
-         tupletable2           const colorfreqtable,
-         unsigned int          const depth,
-         enum methodForLargest const methodForLargest) {
+splitBox(struct BoxVector *    const boxVectorP,
+         unsigned int          const boxIdx,
+         enum MethodForLargest const methodForLargest,
+         enum MethodForSplit   const methodForSplit) {
 /*----------------------------------------------------------------------------
-   Split Box 'bi' in the box vector bv (so that bv contains one more box
-   than it did as input).  Split it so that each new box represents about
-   half of the pixels in the distribution given by 'colorfreqtable' for
-   the colors in the original box, but with distinct colors in each of the
-   two new boxes.
+   Split Box 'boxIdx' in the box vector 'boxVector' (so that 'boxVector'
+   contains one more box than it did as input).  Split it so that each new box
+   represents about half of the pixels in the image for the colors in the
+   original box, but with distinct colors in each of the two new boxes.
 
    Assume the box contains at least two colors.
 -----------------------------------------------------------------------------*/
-    unsigned int const boxStart = bv[bi].ind;
-    unsigned int const boxSize  = bv[bi].colors;
-    unsigned int const sm       = bv[bi].sum;
-
-    sample * minval;  /* malloc'ed array */
-    sample * maxval;  /* malloc'ed array */
+    unsigned int const boxStart = boxVectorP->box[boxIdx].startIndex;
+    unsigned int const boxSize  = boxVectorP->box[boxIdx].colorCt;
+    unsigned int const sum      = boxVectorP->box[boxIdx].sum;
 
-    unsigned int largestDimension;
-        /* number of the plane with the largest spread */
     unsigned int medianIndex;
-    int lowersum;
+    unsigned int lowerSum;
         /* Number of pixels whose value is "less than" the median */
 
-    MALLOCARRAY_NOFAIL(minval, depth);
-    MALLOCARRAY_NOFAIL(maxval, depth);
-
-    findBoxBoundaries(colorfreqtable, depth, boxStart, boxSize,
-                      minval, maxval);
-
-    /* Find the largest dimension, and sort by that component.  I have
-       included two methods for determining the "largest" dimension;
-       first by simply comparing the range in RGB space, and second by
-       transforming into luminosities before the comparison.
-    */
-    switch (methodForLargest) {
-    case LARGE_NORM:
-        largestDimension = largestByNorm(minval, maxval, depth);
-        break;
-    case LARGE_LUM:
-        largestDimension = largestByLuminosity(minval, maxval, depth);
-        break;
-    }
 
-    /* TODO: I think this sort should go after creating a box,
-       not before splitting.  Because you need the sort to use
-       the REP_CENTER_BOX method of choosing a color to
-       represent the final boxes
+    /* Perhaps this sort should go after creating a box, not before splitting.
+       Because you need the sort to use the REP_CENTER_BOX method of choosing
+       a color to represent the final boxes
     */
 
     /* Set the gross global variable 'compareplanePlane' as a
        parameter to compareplane(), which is called by qsort().
     */
-    compareplanePlane = largestDimension;
-    qsort((char*) &colorfreqtable.table[boxStart], boxSize,
-          sizeof(colorfreqtable.table[boxStart]),
+    compareplanePlane = boxVectorP->box[boxIdx].maxdim;
+    qsort((char*) &boxVectorP->colorFreqTable.table[boxStart], boxSize,
+          sizeof(boxVectorP->colorFreqTable.table[boxStart]),
           compareplane);
 
     {
-        /* Now find the median based on the counts, so that about half
-           the pixels (not colors, pixels) are in each subdivision.  */
-
+        /* Find the median based on the counts, so that about half the pixels
+           (not colors, pixels) are in each subdivision.
+        */
         unsigned int i;
 
-        lowersum = colorfreqtable.table[boxStart]->value; /* initial value */
-        for (i = 1; i < boxSize - 1 && lowersum < sm/2; ++i) {
-            lowersum += colorfreqtable.table[boxStart + i]->value;
+        lowerSum = boxVectorP->colorFreqTable.table[boxStart]->value;
+            /* initial value */
+        for (i = 1; i < boxSize - 1 && lowerSum < sum/2; ++i) {
+            lowerSum += boxVectorP->colorFreqTable.table[boxStart + i]->value;
         }
         medianIndex = i;
     }
     /* Split the box, and sort to bring the biggest boxes to the top.  */
+    {
+        struct Box * const oldBoxP = &boxVectorP->box[boxIdx];
+
+        oldBoxP->colorCt = medianIndex;
+        oldBoxP->sum     = lowerSum;
+        computeBoxSpread(oldBoxP, boxVectorP->colorFreqTable,
+                         boxVectorP->colorDepth, methodForLargest,
+                         &oldBoxP->maxdim, &oldBoxP->spread);
+    }
+    {
+        struct Box * const newBoxP = &boxVectorP->box[boxVectorP->boxCt];
+
+        newBoxP->startIndex = boxStart + medianIndex;
+        newBoxP->colorCt    = boxSize - medianIndex;
+        newBoxP->sum        = sum - lowerSum;
+        computeBoxSpread(newBoxP, boxVectorP->colorFreqTable,
+                         boxVectorP->colorDepth, methodForLargest,
+                         &newBoxP->maxdim, &newBoxP->spread);
+        ++boxVectorP->boxCt;
+    }
 
-    bv[bi].colors = medianIndex;
-    bv[bi].sum = lowersum;
-    bv[*boxesP].ind = boxStart + medianIndex;
-    bv[*boxesP].colors = boxSize - medianIndex;
-    bv[*boxesP].sum = sm - lowersum;
-    ++(*boxesP);
-    qsort((char*) bv, *boxesP, sizeof(struct box), sumcompare);
+    sortBoxes(boxVectorP, methodForSplit);
+}
 
-    free(minval); free(maxval);
+
+
+static void
+reportBoxVector(struct BoxVector const boxVector) {
+
+    unsigned int i;
+
+    pm_message("All colors of image, sorted into %u boxes:", boxVector.boxCt);
+
+    for (i = 0; i < boxVector.boxCt; ++i) {
+        const struct Box * const boxP = &boxVector.box[i];
+
+        unsigned int j;
+
+        pm_message("Box %u, %u colors starting with index %u (%u pixels):",
+                   i, boxP->colorCt, boxP->startIndex, boxP->sum);
+        if (boxP->colorCt > 1)
+            pm_message("Largest spread is %lu, in plane %u",
+                       boxP->spread, boxP->maxdim);
+
+        for (j = 0; j < boxP->colorCt; ++j) {
+            unsigned int colorIdx = boxP->startIndex + j;
+
+            assert(colorIdx < boxVector.colorFreqTable.size);
+
+            tuple const color =
+                boxVector.colorFreqTable.table[colorIdx]->tuple;
+
+            pm_message("(%lu, %lu, %lu)",
+                       color[PAM_RED_PLANE],
+                       color[PAM_GRN_PLANE],
+                       color[PAM_BLU_PLANE]);
+        }
+    }
 }
 
 
 
 static void
-mediancut(tupletable2           const colorfreqtable,
+mediancut(tupletable2           const colorFreqTable,
           unsigned int          const depth,
-          int                   const newcolors,
-          enum methodForLargest const methodForLargest,
-          enum methodForRep     const methodForRep,
+          unsigned int          const newColorCt,
+          enum MethodForLargest const methodForLargest,
+          enum MethodForRep     const methodForRep,
+          enum MethodForSplit   const methodForSplit,
+          bool                  const wantBvReport,
           tupletable2 *         const colormapP) {
 /*----------------------------------------------------------------------------
-   Compute a set of only 'newcolors' colors that best represent an
+   Compute a set of only 'newColorCt' colors that best represent an
    image whose pixels are summarized by the histogram
-   'colorfreqtable'.  Each tuple in that table has depth 'depth'.
-   colorfreqtable.table[i] tells the number of pixels in the subject image
-   have a particular color.
+   'colorFreqTable'.  Each tuple in that table has depth 'depth'.
+   colorFreqTable.table[i] tells the number of pixels in the subject image
+   that have a particular color.
 
-   As a side effect, sort 'colorfreqtable'.
+   As a side effect, sort 'colorFreqTable'.
 -----------------------------------------------------------------------------*/
-    boxVector bv;
-    unsigned int bi;
-    unsigned int boxes;
+    struct BoxVector boxVector;
     bool multicolorBoxesExist;
         /* There is at least one box that contains at least 2 colors; ergo,
            there is more splitting we can do.
         */
 
-    bv = newBoxVector(colorfreqtable.size, freqTotal(colorfreqtable),
-                      newcolors);
-    boxes = 1;
-    multicolorBoxesExist = (colorfreqtable.size > 1);
+    boxVector = newBoxVector(colorFreqTable, newColorCt, depth,
+                             methodForLargest);
 
-    /* Main loop: split boxes until we have enough. */
-    while (boxes < newcolors && multicolorBoxesExist) {
-        /* Find the first splittable box. */
-        for (bi = 0; bi < boxes && bv[bi].colors < 2; ++bi);
-        if (bi >= boxes)
+    multicolorBoxesExist = (colorFreqTable.size > 1);
+
+    /* Split boxes until we have enough. */
+    while (boxVector.boxCt < newColorCt && multicolorBoxesExist) {
+        unsigned int boxIdx;
+
+        for (boxIdx = 0;
+             boxIdx < boxVector.boxCt && boxVector.box[boxIdx].colorCt < 2;
+             ++boxIdx);
+            /* Find the first splittable box. */
+
+        if (boxIdx >= boxVector.boxCt)
             multicolorBoxesExist = FALSE;
         else
-            splitBox(bv, &boxes, bi, colorfreqtable, depth, methodForLargest);
+            splitBox(&boxVector, boxIdx, methodForLargest, methodForSplit);
+                /* Side effect: sorts the extent of 'colorfreqTable' that is
+                   in the box
+                */
     }
-    *colormapP = colormapFromBv(newcolors, bv, boxes, colorfreqtable, depth,
-                                methodForRep);
 
-    free(bv);
+    if (wantBvReport)
+        reportBoxVector(boxVector);
+
+    *colormapP = colormapFromBv(newColorCt, boxVector, methodForRep);
+
+    destroyBoxVector(boxVector);
 }
 
 
@@ -673,7 +897,7 @@ static void
 computeHistogram(FILE *         const ifP,
                  int *          const formatP,
                  struct pam *   const freqPamP,
-                 tupletable2 *  const colorfreqtableP) {
+                 tupletable2 *  const colorFreqTableP) {
 /*----------------------------------------------------------------------------
   Make a histogram of the colors in the image stream in the file '*ifP'.
 
@@ -715,13 +939,13 @@ computeHistogram(FILE *         const ifP,
 
         pnm_nextimage(ifP, &eof);
     }
-    colorfreqtableP->table =
+    colorFreqTableP->table =
         pnm_tuplehashtotable(&firstPam, tuplehash, colorCount);
-    colorfreqtableP->size = colorCount;
+    colorFreqTableP->size = colorCount;
 
     pnm_destroytuplehash(tuplehash);
 
-    pm_message("%u colors found", colorfreqtableP->size);
+    pm_message("%u colors found", colorFreqTableP->size);
 
     freqPamP->size   = sizeof(*freqPamP);
     freqPamP->len    = PAM_STRUCT_SIZE(tuple_type);
@@ -738,9 +962,11 @@ computeHistogram(FILE *         const ifP,
 static void
 computeColorMapFromInput(FILE *                const ifP,
                          bool                  const allColors,
-                         int                   const reqColors,
-                         enum methodForLargest const methodForLargest,
-                         enum methodForRep     const methodForRep,
+                         unsigned int          const reqColors,
+                         enum MethodForLargest const methodForLargest,
+                         enum MethodForRep     const methodForRep,
+                         enum MethodForSplit   const methodForSplit,
+                         bool                  const wantBvReport,
                          int *                 const formatP,
                          struct pam *          const freqPamP,
                          tupletable2 *         const colormapP) {
@@ -763,23 +989,26 @@ computeColorMapFromInput(FILE *                const ifP,
    *formatP and *freqPamP.  (This information is not really
    relevant to our colormap mission; just a fringe benefit).
 -----------------------------------------------------------------------------*/
-    tupletable2 colorfreqtable;
+    tupletable2 colorFreqTable;
+        /* Table of all colors in the image, with the number of pixels of
+           each color.
+        */
 
-    computeHistogram(ifP, formatP, freqPamP, &colorfreqtable);
+    computeHistogram(ifP, formatP, freqPamP, &colorFreqTable);
 
     if (allColors) {
-        *colormapP = colorfreqtable;
+        *colormapP = colorFreqTable;
     } else {
-        if (colorfreqtable.size <= reqColors) {
-            pm_message("Image already has few enough colors (<=%d).  "
+        if (colorFreqTable.size <= reqColors) {
+            pm_message("Image already has few enough colors (<=%u).  "
                        "Keeping same colors.", reqColors);
-            *colormapP = colorfreqtable;
+            *colormapP = colorFreqTable;
         } else {
-            pm_message("choosing %d colors...", reqColors);
-            mediancut(colorfreqtable, freqPamP->depth,
+            pm_message("choosing %u colors...", reqColors);
+            mediancut(colorFreqTable, freqPamP->depth,
                       reqColors, methodForLargest, methodForRep,
-                      colormapP);
-            pnm_freetupletable2(freqPamP, colorfreqtable);
+                      methodForSplit, wantBvReport, colormapP);
+            pnm_freetupletable2(freqPamP, colorFreqTable);
         }
     }
 }
@@ -933,9 +1162,9 @@ colormapToImage(int                const format,
 
 
 int
-main(int argc, char * argv[] ) {
+main(int argc, const char * argv[] ) {
 
-    struct cmdlineInfo cmdline;
+    struct CmdlineInfo cmdline;
     FILE * ifP;
     int format;
     struct pam colormapPam;
@@ -943,16 +1172,18 @@ main(int argc, char * argv[] ) {
     tuple ** colormapRaster;
     tupletable2 colormap;
 
-    pnm_init(&argc, argv);
+    pm_proginit(&argc, argv);
 
     parseCommandLine(argc, argv, &cmdline);
 
-    ifP = pm_openr(cmdline.inputFilespec);
+    ifP = pm_openr(cmdline.inputFileNm);
 
     computeColorMapFromInput(ifP,
                              cmdline.allcolors, cmdline.newcolors,
                              cmdline.methodForLargest,
                              cmdline.methodForRep,
+                             cmdline.methodForSplit,
+                             cmdline.debug,
                              &format, &colormapPam, &colormap);
 
     pm_close(ifP);