Promote Development to Advanced, Release 10.94.00

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

15 files changed, 1243 insertions, 381 deletions

diff --git a/editor/Makefile b/editor/Makefile
index 88409dad..395deaf4 100644
--- a/editor/Makefile
+++ b/editor/Makefile
@@ -19,7 +19,8 @@ SUBDIRS = pamflip specialty
 PORTBINARIES = pamaddnoise pamaltsat pambackground pambrighten pamcomp pamcut \
 	       pamdice pamditherbw pamedge \
 	       pamenlarge \
-	       pamfunc pamhue pamlevels pammasksharpen pammixmulti \
+	       pamfunc pamhomography pamhue pamlevels \
+	       pammasksharpen pammixmulti \
 	       pamperspective pamrecolor pamrubber \
 	       pamscale pamsistoaglyph pamstretch pamthreshold pamundice \
 	       pamwipeout \
diff --git a/editor/pamaddnoise.c b/editor/pamaddnoise.c
index 20c68d99..9ca80394 100644
--- a/editor/pamaddnoise.c
+++ b/editor/pamaddnoise.c
@@ -33,18 +33,20 @@
 
 #include "pm_c_util.h"
 #include "mallocvar.h"
+#include "rand.h"
 #include "shhopt.h"
 #include "pm_gamma.h"
 #include "pam.h"
 
 static double const EPSILON = 1.0e-5;
+static double const SALT_RATIO = 0.5;
 
 
 
 static double
-rand1() {
+rand1(struct pm_randSt * const randStP) {
 
-    return (double)rand()/RAND_MAX;
+    return (double)pm_rand(randStP)/RAND_MAX;
 }
 
 
@@ -67,6 +69,7 @@ struct CmdlineInfo {
 
     enum NoiseType noiseType;
 
+    unsigned int seedSpec;
     unsigned int seed;
 
     float lambda;
@@ -119,7 +122,7 @@ parseCommandLine(int argc, const char ** const argv,
 
     unsigned int option_def_index;
 
-    unsigned int typeSpec, seedSpec, lambdaSpec, lsigmaSpec, mgsigmaSpec,
+    unsigned int typeSpec, lambdaSpec, lsigmaSpec, mgsigmaSpec,
         sigma1Spec, sigma2Spec, toleranceSpec;
 
     const char * type;
@@ -128,21 +131,21 @@ parseCommandLine(int argc, const char ** const argv,
 
     option_def_index = 0;   /* incremented by OPTENT3 */
     OPTENT3(0,   "type",            OPT_STRING,   &type,
-            &typeSpec,         0);
+            &typeSpec,           0);
     OPTENT3(0,   "seed",            OPT_UINT,     &cmdlineP->seed,
-            &seedSpec,         0);
+            &cmdlineP->seedSpec, 0);
     OPTENT3(0,   "lambda",          OPT_FLOAT,    &cmdlineP->lambda,
-            &lambdaSpec,       0);
+            &lambdaSpec,         0);
     OPTENT3(0,   "lsigma",          OPT_FLOAT,    &cmdlineP->lsigma,
-            &lsigmaSpec,       0);
+            &lsigmaSpec,         0);
     OPTENT3(0,   "mgsigma",         OPT_FLOAT,    &cmdlineP->mgsigma,
-            &mgsigmaSpec,      0);
+            &mgsigmaSpec,        0);
     OPTENT3(0,   "sigma1",          OPT_FLOAT,    &cmdlineP->sigma1,
-            &sigma1Spec,       0);
+            &sigma1Spec,         0);
     OPTENT3(0,   "sigma2",          OPT_FLOAT,    &cmdlineP->sigma2,
-            &sigma2Spec,       0);
+            &sigma2Spec,         0);
     OPTENT3(0,   "tolerance",       OPT_FLOAT,    &cmdlineP->tolerance,
-            &toleranceSpec,    0);
+            &toleranceSpec,      0);
 
     opt.opt_table = option_def;
     opt.short_allowed = FALSE;  /* We have no short (old-fashioned) options */
@@ -193,7 +196,7 @@ parseCommandLine(int argc, const char ** const argv,
     if (!toleranceSpec)
         cmdlineP->tolerance = 0.10;
 
-    if (!seedSpec)
+    if (!cmdlineP->seedSpec)
         cmdlineP->seed = pm_randseed();
 
     if (argc-1 > 1)
@@ -211,11 +214,12 @@ parseCommandLine(int argc, const char ** const argv,
 
 
 static void
-addGaussianNoise(sample   const maxval,
-                 sample   const origSample,
-                 sample * const newSampleP,
-                 float    const sigma1,
-                 float    const sigma2) {
+addGaussianNoise(sample             const maxval,
+                 sample             const origSample,
+                 sample *           const newSampleP,
+                 float              const sigma1,
+                 float              const sigma2,
+                 struct pm_randSt * const randStP) {
 /*----------------------------------------------------------------------------
    Add Gaussian noise.
 
@@ -225,11 +229,11 @@ addGaussianNoise(sample   const maxval,
     double x1, x2, xn, yn;
     double rawNewSample;
 
-    x1 = rand1();
+    x1 = rand1(randStP);
 
     if (x1 == 0.0)
         x1 = 1.0;
-    x2 = rand1();
+    x2 = rand1(randStP);
     xn = sqrt(-2.0 * log(x1)) * cos(2.0 * M_PI * x2);
     yn = sqrt(-2.0 * log(x1)) * sin(2.0 * M_PI * x2);
 
@@ -242,40 +246,42 @@ addGaussianNoise(sample   const maxval,
 
 
 static void
-addImpulseNoise(sample   const maxval,
-                sample   const origSample,
-                sample * const newSampleP,
-                float    const tolerance) {
+addImpulseNoise(sample             const maxval,
+                sample             const origSample,
+                sample *           const newSampleP,
+                float              const tolerance,
+                double             const saltRatio,
+                struct pm_randSt * const randStP) {
 /*----------------------------------------------------------------------------
    Add impulse (salt and pepper) noise
 -----------------------------------------------------------------------------*/
 
-    double const low_tol  = tolerance / 2.0;
-    double const high_tol = 1.0 - (tolerance / 2.0);
-    double const sap = rand1();
+    double const pepperRatio = 1.0 - saltRatio;
+    double const loTolerance = tolerance * pepperRatio;
+    double const hiTolerance = 1.0 - tolerance * saltRatio;
+    double const sap         = rand1(randStP);
 
-    if (sap < low_tol)
-        *newSampleP = 0;
-    else if ( sap >= high_tol )
-        *newSampleP = maxval;
-    else
-        *newSampleP = origSample;
+    *newSampleP =
+        sap < loTolerance ? 0 :
+        sap >= hiTolerance? maxval :
+        origSample;
 }
 
 
 
 static void
-addLaplacianNoise(sample   const maxval,
-                  double   const infinity,
-                  sample   const origSample,
-                  sample * const newSampleP,
-                  float    const lsigma) {
+addLaplacianNoise(sample             const maxval,
+                  double             const infinity,
+                  sample             const origSample,
+                  sample *           const newSampleP,
+                  float              const lsigma,
+                  struct pm_randSt * const randStP) {
 /*----------------------------------------------------------------------------
    Add Laplacian noise
 
    From Pitas' book.
 -----------------------------------------------------------------------------*/
-    double const u = rand1();
+    double const u = rand1(randStP);
 
     double rawNewSample;
 
@@ -297,11 +303,12 @@ addLaplacianNoise(sample   const maxval,
 
 
 static void
-addMultiplicativeGaussianNoise(sample   const maxval,
-                               double   const infinity,
-                               sample   const origSample,
-                               sample * const newSampleP,
-                               float    const mgsigma) {
+addMultiplicativeGaussianNoise(sample             const maxval,
+                               double             const infinity,
+                               sample             const origSample,
+                               sample *           const newSampleP,
+                               float              const mgsigma,
+                               struct pm_randSt * const randStP) {
 /*----------------------------------------------------------------------------
    Add multiplicative Gaussian noise
 
@@ -311,14 +318,14 @@ addMultiplicativeGaussianNoise(sample   const maxval,
     double rawNewSample;
 
     {
-        double const uniform = rand1();
+        double const uniform = rand1(randStP);
         if (uniform <= EPSILON)
             rayleigh = infinity;
         else
             rayleigh = sqrt(-2.0 * log( uniform));
     }
     {
-        double const uniform = rand1();
+        double const uniform = rand1(randStP);
         gauss = rayleigh * cos(2.0 * M_PI * uniform);
     }
     rawNewSample = origSample + (origSample * mgsigma * gauss);
@@ -363,10 +370,11 @@ poissonPmf(double       const lambda,
 
 
 static void
-addPoissonNoise(struct pam * const pamP,
-                sample       const origSample,
-                sample *     const newSampleP,
-                float        const lambdaOfMaxval) {
+addPoissonNoise(struct pam *       const pamP,
+                sample             const origSample,
+                sample *           const newSampleP,
+                float              const lambdaOfMaxval,
+                struct pm_randSt * const randStP) {
 /*----------------------------------------------------------------------------
    Add Poisson noise
 -----------------------------------------------------------------------------*/
@@ -376,7 +384,7 @@ addPoissonNoise(struct pam * const pamP,
 
     double const lambda  = origSampleIntensity * lambdaOfMaxval;
 
-    double const u = rand1();
+    double const u = rand1(randStP);
 
     /* We now apply the inverse CDF (cumulative distribution function) of the
        Poisson distribution to uniform random variable 'u' to get a Poisson
@@ -416,12 +424,14 @@ main(int argc, const char ** argv) {
     const tuple * newtuplerow;
     unsigned int row;
     double infinity;
+    struct pm_randSt randSt;
 
     pm_proginit(&argc, argv);
 
     parseCommandLine(argc, argv, &cmdline);
 
-    srand(cmdline.seed);
+    pm_randinit(&randSt);
+    pm_srand2(&randSt, cmdline.seedSpec, cmdline.seed);
 
     ifP = pm_openr(cmdline.inputFileName);
 
@@ -448,35 +458,40 @@ main(int argc, const char ** argv) {
                     addGaussianNoise(inpam.maxval,
                                      tuplerow[col][plane],
                                      &newtuplerow[col][plane],
-                                     cmdline.sigma1, cmdline.sigma2);
+                                     cmdline.sigma1, cmdline.sigma2,
+                                     &randSt);
                     break;
 
                 case NOISETYPE_IMPULSE:
                     addImpulseNoise(inpam.maxval,
                                     tuplerow[col][plane],
                                     &newtuplerow[col][plane],
-                                    cmdline.tolerance);
+                                    cmdline.tolerance, SALT_RATIO,
+                                    &randSt);
                    break;
 
                 case NOISETYPE_LAPLACIAN:
                     addLaplacianNoise(inpam.maxval, infinity,
                                       tuplerow[col][plane],
                                       &newtuplerow[col][plane],
-                                      cmdline.lsigma);
+                                      cmdline.lsigma,
+                                      &randSt);
                     break;
 
                 case NOISETYPE_MULTIPLICATIVE_GAUSSIAN:
                     addMultiplicativeGaussianNoise(inpam.maxval, infinity,
                                                    tuplerow[col][plane],
                                                    &newtuplerow[col][plane],
-                                                   cmdline.mgsigma);
+                                                   cmdline.mgsigma,
+                                                   &randSt);
                     break;
 
                 case NOISETYPE_POISSON:
                     addPoissonNoise(&inpam,
                                     tuplerow[col][plane],
                                     &newtuplerow[col][plane],
-                                    cmdline.lambda);
+                                    cmdline.lambda,
+                                    &randSt);
                     break;
 
                 }
@@ -484,6 +499,7 @@ main(int argc, const char ** argv) {
         }
         pnm_writepamrow(&outpam, newtuplerow);
     }
+    pm_randterm(&randSt);
     pnm_freepamrow(newtuplerow);
     pnm_freepamrow(tuplerow);
 
diff --git a/editor/pamditherbw.c b/editor/pamditherbw.c
index ae91a26f..694b2c21 100644
--- a/editor/pamditherbw.c
+++ b/editor/pamditherbw.c
@@ -14,12 +14,14 @@
 #include <string.h>
 
 #include "pm_c_util.h"
-#include "pam.h"
-#include "dithers.h"
+#include "rand.h"
 #include "mallocvar.h"
 #include "shhopt.h"
+#include "pam.h"
+#include "dithers.h"
 #include "pm_gamma.h"
 
+
 enum halftone {QT_FS,
                QT_ATKINSON,
                QT_THRESH,
@@ -588,7 +590,9 @@ fsDestroy(struct converter * const converterP) {
 
 static struct converter
 createFsConverter(struct pam * const graypamP,
-                  float        const threshFraction) {
+                  float        const threshFraction,
+                  bool         const randomseedSpec,
+                  unsigned int const randomseed) {
 
     struct fsState * stateP;
     struct converter converter;
@@ -605,9 +609,17 @@ createFsConverter(struct pam * const graypamP,
 
     {
         /* (random errors in [-1/8 .. 1/8]) */
+
         unsigned int col;
+        struct pm_randSt randSt;
+
+        pm_randinit(&randSt);
+        pm_srand2(&randSt, randomseedSpec, randomseed);
+
         for (col = 0; col < graypamP->width + 2; ++col)
-            stateP->thiserr[col] = ((float)rand()/RAND_MAX - 0.5) / 4;
+            stateP->thiserr[col] = (pm_drand(&randSt) - 0.5) / 4;
+
+        pm_randterm(&randSt);
     }
 
     stateP->halfWhite = threshFraction;
@@ -725,7 +737,9 @@ atkinsonDestroy(struct converter * const converterP) {
 
 static struct converter
 createAtkinsonConverter(struct pam * const graypamP,
-                        float        const threshFraction) {
+                        float        const threshFraction,
+                        bool         const randomseedSpec,
+                        unsigned int const randomseed) {
 
     struct atkinsonState * stateP;
     struct converter converter;
@@ -743,11 +757,18 @@ createAtkinsonConverter(struct pam * const graypamP,
     {
         /* (random errors in [-1/8 .. 1/8]) */
         unsigned int col;
+        struct pm_randSt randSt;
+
+        pm_randinit(&randSt);
+        pm_srand2(&randSt, randomseedSpec, randomseed);
+
         for (col = 0; col < graypamP->width + 2; ++col) {
-            stateP->error[0][col] = ((float)rand()/RAND_MAX - 0.5) / 4;
+            stateP->error[0][col] = (pm_drand(&randSt) - 0.5) / 4;
             stateP->error[1][col] = 0.0;
             stateP->error[2][col] = 0.0;
         }
+
+        pm_randterm(&randSt);
     }
 
     stateP->halfWhite = threshFraction;
@@ -934,8 +955,6 @@ main(int argc, char *argv[]) {
 
     parseCommandLine(argc, argv, &cmdline);
 
-    srand(cmdline.randomseedSpec ? cmdline.randomseed : pm_randseed());
-
     ifP = pm_openr(cmdline.inputFilespec);
 
     if (cmdline.halftone == QT_HILBERT)
@@ -956,10 +975,14 @@ main(int argc, char *argv[]) {
 
         switch (cmdline.halftone) {
         case QT_FS:
-            converter = createFsConverter(&graypam, cmdline.threshval);
+            converter = createFsConverter(&graypam, cmdline.threshval,
+                                          cmdline.randomseedSpec,
+                                          cmdline.randomseed);
             break;
         case QT_ATKINSON:
-            converter = createAtkinsonConverter(&graypam, cmdline.threshval);
+            converter = createAtkinsonConverter(&graypam, cmdline.threshval,
+                                                cmdline.randomseedSpec,
+                                                cmdline.randomseed);
             break;
         case QT_THRESH:
             converter = createThreshConverter(&graypam, cmdline.threshval);
diff --git a/editor/pamhomography.c b/editor/pamhomography.c
new file mode 100644
index 00000000..59b59ed7
--- /dev/null
+++ b/editor/pamhomography.c
@@ -0,0 +1,677 @@
+/* ----------------------------------------------------------------------
+ *
+ * Map one quadrilateral to another
+ * by Scott Pakin <scott+pbm@pakin.org>
+ *
+ * ----------------------------------------------------------------------
+ *
+ * Copyright (C) 2020 Scott Pakin <scott+pbm@pakin.org>
+ *
+ * All rights reserved.
+ *
+ * Redistribution and use in source and binary forms, with or without
+ * modification, are permitted provided that the following conditions
+ * are met:
+ *
+ * 1. Redistributions of source code must retain the above copyright
+ *    notice, this list of conditions and the following disclaimer.
+ * 2. Redistributions in binary form must reproduce the above copyright
+ *    notice, this list of conditions and the following disclaimer in the
+ *    documentation and/or other materials provided with the distribution.
+ * 3. The name of the author may not be used to endorse or promote products
+ *    derived from this software without specific prior written permission.
+ *
+ * THIS SOFTWARE IS PROVIDED BY THE AUTHOR "AS IS" AND ANY EXPRESS OR
+ * IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES
+ * OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
+ * IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY DIRECT, INDIRECT,
+ * INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT
+ * NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
+ * DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
+ * THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
+ * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF
+ * THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+ *
+ * ----------------------------------------------------------------------
+ */
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include <ctype.h>
+#include <errno.h>
+#include <math.h>
+#include "pm_c_util.h"
+#include "shhopt.h"
+#include "mallocvar.h"
+#include "pam.h"
+
+#define MIN4(A, B, C, D) MIN(MIN(A, B), MIN(C, D))
+#define MAX4(A, B, C, D) MAX(MAX(A, B), MAX(C, D))
+
+/* A point on the image plane.  It may or may not lie within the
+   bounds of the image itself. */
+typedef struct Point {
+    int x;
+    int y;
+} Point;
+
+/* A quadrilateral on the image plane */
+typedef struct Quad {
+    Point ul;
+    Point ur;
+    Point lr;
+    Point ll;
+} Quad;
+
+/* A user-specified mapping from one quadrilateral to another */
+typedef struct QuadMap {
+    Quad from;
+    Quad to;
+} QuadMap;
+
+struct CmdlineInfo {
+    /* All the information the user supplied in the command line,
+       in a form easy for the program to use.
+    */
+    const char * inputFilespec;  /* "-" if stdin */
+    QuadMap      qmap;           /* Source and target quadrilaterals */
+    Quad         bbox;           /* Bounding box for the target image */
+    const char * fillColor;      /* Fill color for unused coordinates */
+};
+
+
+
+static unsigned int
+parseCoords(const char * const str,
+            int *        const coords) {
+/*----------------------------------------------------------------------------
+  Parse a list of up to 16 integers.  The function returns the number
+  of integers encountered.
+-----------------------------------------------------------------------------*/
+
+    const char * p;
+    char * pnext;
+    unsigned int i;
+
+    for (i = 0, p = str; i < 16; ++i, p = pnext) {
+        long int val;
+
+        /* Skip punctuation, except "+" and "-", and white space. */
+        while (*p != '\0' && *p != '+' && *p != '-' &&
+               (isspace(*p) || ispunct(*p)))
+            ++p;
+
+        /* Parse the next integer. */
+        errno = 0;  /* strtol() sets errno on error. */
+        val = strtol(p, &pnext, 10);
+        if (errno == ERANGE)
+            return i;  /* Integer lies out of long int range */
+        if (errno != 0 || pnext == p)
+            return i;  /* Too few integers */
+        coords[i] = (int)val;
+        if ((long int)coords[i] != val)
+            return i;  /* Integer lies out of int range */
+    }
+    return i;
+}
+
+
+
+static void
+parseViewString(const char * const str,
+                Quad *       const quad) {
+/*----------------------------------------------------------------------------
+  Parse a list of four integers in the order {ulx, uly, lrx, lry} into a
+  quadrilateral.  The function aborts on error.
+-----------------------------------------------------------------------------*/
+
+    int coords[16];
+
+    if (parseCoords(str, coords) != 4)
+        pm_error("failed to parse \"%s\" as a list of four integers", str);
+    quad->ul.x = quad->ll.x = coords[0];
+    quad->ul.y = quad->ur.y = coords[1];
+    quad->lr.x = quad->ur.x = coords[2];
+    quad->lr.y = quad->ll.y = coords[3];
+}
+
+
+
+static void
+parseQuadString(const char * const str,
+                Quad *       const quad) {
+/*----------------------------------------------------------------------------
+  Parse a list of eight integers in the order {ulx, uly, urx, ury,
+  lrx, lry, llx, lly} into a quadrilateral.  The function aborts on
+  error.
+-----------------------------------------------------------------------------*/
+
+    int coords[16];
+
+    if (parseCoords(str, coords) != 8)
+        pm_error("failed to parse \"%s\" as a list of eight integers", str);
+    quad->ul.x = coords[0];
+    quad->ul.y = coords[1];
+    quad->ur.x = coords[2];
+    quad->ur.y = coords[3];
+    quad->lr.x = coords[4];
+    quad->lr.y = coords[5];
+    quad->ll.x = coords[6];
+    quad->ll.y = coords[7];
+}
+
+
+
+static void
+readMapFile(const char * const fname,
+            QuadMap *    const qmap) {
+/*----------------------------------------------------------------------------
+  Read from a file either 16 numbers in the order {ulx1, uly1, urx1, ury1,
+  lrx1, lry1, llx1, lly1, ulx2, uly2, urx2, ury2, lrx2, lry2, llx2, lly2}
+  or 8 numbers in the order {ulx2, uly2, urx2, ury2, lrx2, lry2, llx2,
+  lly2}.  This function aborts on error.
+-----------------------------------------------------------------------------*/
+
+    FILE * fp;
+    char * str;      /* Entire file contents */
+    int coords[16];  /* File as a list of up to 16 coordinates */
+    char * c;
+    long int nread;
+
+    /* Read the entire file. */
+    fp = pm_openr(fname);
+    str = pm_read_unknown_size(fp, &nread);
+    REALLOCARRAY_NOFAIL(str, nread + 1);
+    str[nread] = '\0';
+    pm_close(fp);
+
+    /* Replace newlines and tabs with spaces to prettify error reporting. */
+    for (c = str; *c != '\0'; ++c)
+        if (isspace(*c))
+            *c = ' ';
+
+    /* Read either {from, to} or just a {to} quadrilateral. */
+    switch (parseCoords(str, coords)) {
+    case 16:
+        /* 16 integers: assign both the "from" and the "to" quadrilateral. */
+        qmap->from.ul.x = coords[0];
+        qmap->from.ul.y = coords[1];
+        qmap->from.ur.x = coords[2];
+        qmap->from.ur.y = coords[3];
+        qmap->from.lr.x = coords[4];
+        qmap->from.lr.y = coords[5];
+        qmap->from.ll.x = coords[6];
+        qmap->from.ll.y = coords[7];
+        qmap->to.ul.x = coords[8];
+        qmap->to.ul.y = coords[9];
+        qmap->to.ur.x = coords[10];
+        qmap->to.ur.y = coords[11];
+        qmap->to.lr.x = coords[12];
+        qmap->to.lr.y = coords[13];
+        qmap->to.ll.x = coords[14];
+        qmap->to.ll.y = coords[15];
+        break;
+    case 8:
+        /* 8 integers: assign only the "to" quadrilateral. */
+        memset((void *)&qmap->from, 0, sizeof(Quad));
+        qmap->to.ul.x = coords[0];
+        qmap->to.ul.y = coords[1];
+        qmap->to.ur.x = coords[2];
+        qmap->to.ur.y = coords[3];
+        qmap->to.lr.x = coords[4];
+        qmap->to.lr.y = coords[5];
+        qmap->to.ll.x = coords[6];
+        qmap->to.ll.y = coords[7];
+        break;
+    default:
+        /* Any other number of integers: issue an error message. */
+        pm_error("failed to parse \"%s\" as a list of either 8 or 16 integers",
+                 str);
+        break;
+    }
+
+    free(str);
+}
+
+
+
+static void
+parseCommandLine(int                  argc,
+                 const char **        const 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.
+
+   If command line is internally inconsistent (invalid options, etc.),
+   issue error message to stderr and abort program.
+-----------------------------------------------------------------------------*/
+
+    optEntry *option_def;
+        /* Instructions to pm_optParseOptions3 on how to parse our options.
+         */
+    optStruct3 opt;
+
+    unsigned int option_def_index;
+
+    unsigned int mapFileSpec = 0, fromSpec = 0, toSpec = 0;
+    unsigned int viewSpec = 0, fillColorSpec = 0;
+    char *mapFile, *from, *to, *view;
+
+    MALLOCARRAY_NOFAIL(option_def, 100);
+
+    option_def_index = 0;   /* incremented by OPTENT3 */
+    OPTENT3(0, "mapfile",         OPT_STRING, &mapFile,
+            &mapFileSpec,             0);
+    OPTENT3(0, "from",            OPT_STRING, &from,
+            &fromSpec,                0);
+    OPTENT3(0, "to",              OPT_STRING, &to,
+            &toSpec,                  0);
+    OPTENT3(0, "view",            OPT_STRING, &view,
+            &viewSpec,                0);
+    OPTENT3(0, "fill",            OPT_STRING, &cmdlineP->fillColor,
+            &fillColorSpec,           0);
+
+    opt.opt_table = option_def;
+    opt.short_allowed = FALSE;  /* We have no short (old-fashioned) options */
+    opt.allowNegNum = FALSE;  /* We have no parms that are negative numbers */
+
+    pm_optParseOptions3(&argc, (char **)argv, opt, sizeof(opt), 0);
+        /* Uses and sets argc, argv, and local variables. */
+
+    if (!fillColorSpec)
+        cmdlineP->fillColor = NULL;
+
+    memset((void *)&cmdlineP->qmap, 0, sizeof(QuadMap));
+    if (mapFileSpec)
+        readMapFile(mapFile, &cmdlineP->qmap);
+    if (fromSpec)
+        parseQuadString(from, &cmdlineP->qmap.from);
+    if (toSpec)
+        parseQuadString(to, &cmdlineP->qmap.to);
+    if (!mapFileSpec && !fromSpec && !toSpec && !viewSpec)
+        pm_error("You must specify at least one of "
+                 "-mapfile, -qin, -qout, and -view");
+    if (viewSpec)
+        parseViewString(view, &cmdlineP->bbox);
+    else
+        memset((void *)&cmdlineP->bbox, 0, sizeof(Quad));
+
+    if (argc < 2)
+        cmdlineP->inputFilespec = "-";
+    else if (argc == 2)
+        cmdlineP->inputFilespec = argv[1];
+    else
+        pm_error("Too many non-option arguments: %u.  "
+                 "Only argument is input file name", argc - 1);
+
+    free((void *) option_def);
+}
+
+
+
+static tuple
+parseFillColor(const struct pam *         const pamP,
+               const struct CmdlineInfo * const cmdlineP) {
+/*----------------------------------------------------------------------------
+  Parse the fill color into the correct format for the given PAM metadata.
+-----------------------------------------------------------------------------*/
+
+    tuple rgb;
+    tuple fillColor;
+
+    if (!cmdlineP->fillColor) {
+        pnm_createBlackTuple(pamP, &fillColor);
+        return fillColor;
+    }
+
+    rgb = pnm_parsecolor(cmdlineP->fillColor, pamP->maxval);
+    fillColor = pnm_allocpamtuple(pamP);
+    switch (pamP->depth) {
+    case 1:
+        /* Grayscale */
+        fillColor[0] = (rgb[PAM_RED_PLANE]*299 +
+                        rgb[PAM_GRN_PLANE]*587 +
+                        rgb[PAM_BLU_PLANE]*114)/1000;
+        break;
+    case 2:
+        /* Grayscale + alpha */
+        fillColor[0] = (rgb[PAM_RED_PLANE]*299 +
+                        rgb[PAM_GRN_PLANE]*587 +
+                        rgb[PAM_BLU_PLANE]*114)/1000;
+        fillColor[PAM_GRAY_TRN_PLANE] = pamP->maxval;
+        break;
+    case 3:
+        /* RGB */
+        pnm_assigntuple(pamP, fillColor, rgb);
+        break;
+    case 4:
+        /* RGB + alpha */
+        pnm_assigntuple(pamP, fillColor, rgb);
+        fillColor[PAM_TRN_PLANE] = pamP->maxval;
+        break;
+    default:
+        pm_error("unexpected image depth %d", pamP->depth);
+        break;
+    }
+
+    return fillColor;
+}
+
+
+
+static tuple **
+initOutputImage(const struct pam *         const pamP,
+                const struct CmdlineInfo * const cmdlineP) {
+/*----------------------------------------------------------------------------
+  Allocate and initialize the output image.
+-----------------------------------------------------------------------------*/
+
+    tuple fillColor;    /* Fill color to use for unused coordinates */
+    tuple ** outImg;    /* Output image */
+    unsigned int row;
+
+    outImg = pnm_allocpamarray(pamP);
+
+    fillColor = parseFillColor(pamP, cmdlineP);
+    for (row = 0; row < pamP->height; ++row) {
+        unsigned int col;
+
+        for (col = 0; col < pamP->width; ++col) {
+            pnm_assigntuple(pamP, outImg[row][col], fillColor);
+        }
+    }
+
+    free((void *) fillColor);
+    return outImg;
+}
+
+
+
+static void
+computeSteps(const Quad * const qfrom,
+             const Quad * const qto,
+             double *     const ustep,
+             double *     const vstep) {
+/*----------------------------------------------------------------------------
+  Compute increments for u and v as these range from 0.0 to 1.0.
+-----------------------------------------------------------------------------*/
+
+    double fx0, fx1, fxd;
+    double tx0, tx1, txd;
+    double fy0, fy1, fyd;
+    double ty0, ty1, tyd;
+
+    /* Compute ustep as the inverse of the maximum possible x delta across
+       either the "from" or "to" quadrilateral. */
+    fx0 = MIN4((double)qfrom->ur.x,
+               (double)qfrom->ul.x,
+               (double)qfrom->lr.x,
+               (double)qfrom->ll.x);
+    fx1 = MAX4((double)qfrom->ur.x,
+               (double)qfrom->ul.x,
+               (double)qfrom->lr.x,
+               (double)qfrom->ll.x);
+    fxd = fx1 - fx0;
+    tx0 = MIN4((double)qto->ur.x,
+               (double)qto->ul.x,
+               (double)qto->lr.x,
+               (double)qto->ll.x);
+    tx1 = MAX4((double)qto->ur.x,
+               (double)qto->ul.x,
+               (double)qto->lr.x,
+               (double)qto->ll.x);
+    txd = tx1 - tx0;
+    if (fxd == 0.0 && txd == 0.0)
+        *ustep = 1.0;  /* Arbitrary nonzero step */
+    *ustep = 0.5/MAX(fxd, txd);
+        /* Divide into 0.5 instead of 1.0 for additional smoothing. */
+
+    /* Compute vstep as the inverse of the maximum possible y delta across
+       either the "from" or "to" quadrilateral
+  . */
+    fy0 = MIN4((double)qfrom->ur.y,
+               (double)qfrom->ul.y,
+               (double)qfrom->lr.y,
+               (double)qfrom->ll.y);
+    fy1 = MAX4((double)qfrom->ur.y,
+               (double)qfrom->ul.y,
+               (double)qfrom->lr.y,
+               (double)qfrom->ll.y);
+    fyd = fy1 - fy0;
+    ty0 = MIN4((double)qto->ur.y,
+               (double)qto->ul.y,
+               (double)qto->lr.y,
+               (double)qto->ll.y);
+    ty1 = MAX4((double)qto->ur.y,
+               (double)qto->ul.y,
+               (double)qto->lr.y,
+               (double)qto->ll.y);
+    tyd = ty1 - ty0;
+    if (fyd == 0.0 && tyd == 0.0)
+        *vstep = 1.0;  /* Arbitrary nonzero step */
+    *vstep = 0.5/MAX(fyd, tyd);
+        /* Divide into 0.5 instead of 1.0 for additional smoothing. */
+}
+
+
+
+static Quad *
+prepareQuadrilateral(const struct pam * const pamP,
+                     const Quad *       const qdata) {
+/*----------------------------------------------------------------------------
+  If a quadrilateral has all zero points, replace it with a quadrilateral
+  of the full size of the image.  The caller should free the result.
+-----------------------------------------------------------------------------*/
+
+    Quad * qcopy;
+
+    MALLOCVAR_NOFAIL(qcopy);
+
+    if (qdata->ul.x == 0 && qdata->ul.y == 0 &&
+        qdata->ur.x == 0 && qdata->ur.y == 0 &&
+        qdata->ll.x == 0 && qdata->ll.y == 0 &&
+        qdata->lr.x == 0 && qdata->lr.y == 0) {
+        /* Set the quadrilateral to the image's bounding box. */
+        memset((void *)qcopy, 0, sizeof(Quad));
+        qcopy->ur.x = pamP->width - 1;
+        qcopy->lr.x = pamP->width - 1;
+        qcopy->lr.y = pamP->height - 1;
+        qcopy->ll.y = pamP->height - 1;
+    } else {
+        /* Use the quadrilateral as specified. */
+        memcpy(qcopy, qdata, sizeof(Quad));
+    }
+
+    return qcopy;
+}
+
+
+static void
+coordsAtPercent(const Quad * const quad,
+                double       const u,
+                double       const v,
+                int *        const x,
+                int *        const y) {
+/*----------------------------------------------------------------------------
+  Return the (x, y) coordinates that lie at (u%, v%) from the upper left to
+  the lower right of a given quadrilateral.
+-----------------------------------------------------------------------------*/
+
+    *x = (int) nearbyint((1.0 - u)*(1.0 - v)*quad->ul.x +
+                         u*(1.0 - v)*quad->ur.x +
+                         u*v*quad->lr.x +
+                         (1.0 - u)*v*quad->ll.x);
+    *y = (int) nearbyint((1.0 - u)*(1.0 - v)*quad->ul.y +
+                         u*(1.0 - v)*quad->ur.y +
+                         u*v*quad->lr.y +
+                         (1.0 - u)*v*quad->ll.y);
+}
+
+
+
+static void
+computeBoundingBox(const Quad * const q,
+                   Quad *       const bbox) {
+/*----------------------------------------------------------------------------
+  Compute the bounding box of a given quadrilateral.
+-----------------------------------------------------------------------------*/
+
+    bbox->ul.x = bbox->ll.x = MIN4(q->ul.x, q->ur.x, q->lr.x, q->ll.x);
+    bbox->ul.y = bbox->ur.y = MIN4(q->ul.y, q->ur.y, q->lr.y, q->ll.y);
+    bbox->ur.x = bbox->lr.x = MAX4(q->ul.x, q->ur.x, q->lr.x, q->ll.x);
+    bbox->ll.y = bbox->lr.y = MAX4(q->ul.y, q->ur.y, q->lr.y, q->ll.y);
+}
+
+
+
+static void
+mapQuadrilaterals(const struct pam * const inPamP,
+                  const struct pam * const outPamP,
+                  const Quad *       const qfrom,
+                  const Quad *       const qto,
+                  tuple **           const inImg,
+                  tuple **           const outImg,
+                  int                const xofs,
+                  int                const yofs) {
+/*----------------------------------------------------------------------------
+  Map the quadrilateral in the source image to the quadrilateral in the
+  target image.  This is the function that implemens pamhomography's
+  primary functionality.
+-----------------------------------------------------------------------------*/
+
+    sample ** channel;
+        /* Aggregated values for a single channel */
+    unsigned long ** tally;
+        /* Number of values at each coordinate in the above */
+    double ustep, vstep;
+        /* Steps to use when iterating from 0.0 to 1.0 */
+    double u, v;
+    unsigned int plane, row, col;
+
+    MALLOCARRAY2_NOFAIL(channel, outPamP->height, outPamP->width);
+    MALLOCARRAY2_NOFAIL(tally, outPamP->height, outPamP->width);
+
+    computeSteps(qfrom, qto, &ustep, &vstep);
+
+    for (plane = 0; plane < outPamP->depth; ++plane) {
+        /* Reset the channel colors and tally for each plane, */
+        for (row = 0; row < outPamP->height; ++row)
+            for (col = 0; col < outPamP->width; ++col) {
+                channel[row][col] = 0;
+                tally[row][col] = 0;
+            }
+
+        /* Iterate from 0% to 100% in the y dimension. */
+        for (v = 0.0; v <= 1.0; v += vstep) {
+            /* Iterate from 0% to 100% in the x dimension. */
+            for (u = 0.0; u <= 1.0; u += ustep) {
+                int x0, y0;  /* "From" coordinate */
+                int x1, y1;  /* "To" coordinate */
+
+                /* Map (u%, v%) of one quadrilateral to (u%, v%) of the
+                   other quadrilateral. */
+                coordsAtPercent(qfrom, u, v, &x0, &y0);
+                coordsAtPercent(qto, u, v, &x1, &y1);
+
+                /* Copy the source image's (x0, y0) to the destination
+                   image's (x1, y1) in the current plane. */
+                x1 += xofs;
+                y1 += yofs;
+                if (x0 >= 0 && y0 >= 0 &&
+                    x0 < inPamP->width && y0 < inPamP->height &&
+                    x1 >= 0 && y1 >= 0 &&
+                    x1 < outPamP->width && y1 < outPamP->height) {
+                    channel[y1][x1] += inImg[y0][x0][plane];
+                    tally[y1][x1]++;
+                }
+            }
+        }
+
+        /* Assign the current plane in the output image the average color
+           at each point. */
+        for (row = 0; row < outPamP->height; ++row)
+            for (col = 0; col < outPamP->width; ++col)
+                if (tally[row][col] != 0)
+                    outImg[row][col][plane] =
+                        (channel[row][col] + tally[row][col]/2) /
+                        tally[row][col];
+    }
+
+    pm_freearray2((void ** const)tally);
+    pm_freearray2((void ** const)channel);
+    free((void *)qto);
+    free((void *)qfrom);
+}
+
+
+
+static void
+processFile(FILE *                     const ifP,
+            const struct CmdlineInfo * const cmdlineP) {
+/*----------------------------------------------------------------------------
+  Read the input image, create the output image, and map a quadrilateral in
+  the former to a quadrilateral in the latter.
+-----------------------------------------------------------------------------*/
+
+    struct pam inPam;    /* PAM metadata for the input file */
+    struct pam outPam;   /* PAM metadata for the output file */
+    tuple ** inImg;      /* Input image */
+    tuple ** outImg;     /* Output image */
+    Quad *qfrom, *qto;   /* Source and target quadrilaterals */
+    Quad bbox;           /* Bounding box around the transformed input image */
+
+    inImg = pnm_readpam(ifP, &inPam, PAM_STRUCT_SIZE(tuple_type));
+
+    /* Extract quadrilaterals and populate them with the image bounds
+       if necessary. */
+    qfrom = prepareQuadrilateral(&inPam, &cmdlineP->qmap.from);
+    qto = prepareQuadrilateral(&inPam, &cmdlineP->qmap.to);
+
+    /* Allocate storage for the target image. */
+    if (cmdlineP->bbox.ul.x == 0 && cmdlineP->bbox.ul.y == 0 &&
+        cmdlineP->bbox.lr.x == 0 && cmdlineP->bbox.lr.y == 0)
+        /* User did not specify a target bounding box.  Compute optimal
+           dimensions. */
+        computeBoundingBox(qto, &bbox);
+    else
+        /* User specified a target bounding box.  Use it. */
+        bbox = cmdlineP->bbox;
+    outPam = inPam;
+    outPam.file = stdout;
+    outPam.width = bbox.lr.x - bbox.ul.x + 1;
+    outPam.height = bbox.lr.y - bbox.ul.y + 1;
+    outImg = initOutputImage(&outPam, cmdlineP);
+
+    mapQuadrilaterals(&inPam, &outPam,
+                      qfrom, qto,
+                      inImg, outImg,
+                      -bbox.ul.x, -bbox.ul.y);
+
+    pnm_writepam(&outPam, outImg);
+
+    pnm_freepamarray(outImg, &outPam);
+    pnm_freepamarray(inImg, &inPam);
+}
+
+
+
+int
+main(int argc, const char *argv[]) {
+
+    struct CmdlineInfo cmdline;      /* Parsed command line */
+    FILE * ifP;
+
+    pm_proginit(&argc, argv);
+
+    parseCommandLine(argc, argv, &cmdline);
+
+    ifP = pm_openr(cmdline.inputFilespec);
+
+    processFile(ifP, &cmdline);
+
+    pm_close(ifP);
+
+    return 0;
+}
diff --git a/editor/pammixmulti.c b/editor/pammixmulti.c
index f5012d7a..2b45d807 100644
--- a/editor/pammixmulti.c
+++ b/editor/pammixmulti.c
@@ -13,6 +13,7 @@
 #include "shhopt.h"
 #include "mallocvar.h"
 #include "nstring.h"
+#include "rand.h"
 
 typedef enum {
     BLEND_AVERAGE,   /* Take the average color of all pixels */
@@ -36,6 +37,8 @@ struct ProgramState {
         /* Standard deviation when selecting images via a mask */
     unsigned long ** imageWeights;
         /* Per-image weights as a function of grayscale level */
+    struct pm_randSt randSt;
+        /* Random number generator parameters and internal state */
 };
 
 
@@ -134,9 +137,9 @@ parseCommandLine(int argc, const char ** argv,
 }
 
 static void
-openInputFiles(unsigned int          const inFileCt,
-               const char **         const inFileName,
-               struct ProgramState * const stateP) {
+initInput(unsigned int          const inFileCt,
+          const char **         const inFileName,
+          struct ProgramState * const stateP) {
 /*----------------------------------------------------------------------------
   Open all of the input files.
 
@@ -180,6 +183,25 @@ openInputFiles(unsigned int          const inFileCt,
 }
 
 
+
+static void
+termInput(struct ProgramState * const stateP) {
+/*----------------------------------------------------------------------------
+  Deallocate all of the resources we allocated.
+-----------------------------------------------------------------------------*/
+    unsigned int i;
+
+    for (i = 0; i < stateP->inFileCt; ++i) {
+        pnm_freepamrow(stateP->inTupleRows[i]);
+        pm_close(stateP->inPam[i].file);
+    }
+
+    free(stateP->inTupleRows);
+    free(stateP->inPam);
+}
+
+
+
 static void
 initMask(const char *          const maskFileName,
          struct ProgramState * const stateP) {
@@ -233,6 +255,15 @@ initOutput(FILE *                const ofP,
 }
 
 
+static void
+termOutput(struct ProgramState * const stateP) {
+
+    free(stateP->outTupleRow);
+
+    pm_close(stateP->outPam.file);
+}
+
+
 
 static void
 blendTuplesRandom(struct ProgramState * const stateP,
@@ -243,8 +274,8 @@ blendTuplesRandom(struct ProgramState * const stateP,
   from a random input image.
 -----------------------------------------------------------------------------*/
     unsigned int const depth = stateP->inPam[0].depth;
-    unsigned int const img = (unsigned int) (rand() % stateP->inFileCt);
-
+    unsigned int const img = (unsigned int) (pm_rand(&stateP->randSt) %
+                                             stateP->inFileCt);
     unsigned int samp;
 
     for (samp = 0; samp < depth; ++samp)
@@ -276,23 +307,26 @@ blendTuplesAverage(struct ProgramState * const stateP,
 
 
 
+#if 0
 static void
-randomNormal2(double * const r1P,
-              double * const r2P) {
+randomNormal2(double *           const r1P,
+              double *           const r2P,
+              struct pm_randSt * const randStP) {
 /*----------------------------------------------------------------------------
   Return two normally distributed random numbers.
 -----------------------------------------------------------------------------*/
     double u1, u2;
 
     do {
-        u1 = drand48();
-        u2 = drand48();
+        u1 = drand48(randStP);
+        u2 = drand48(randStP);
     }
     while (u1 <= DBL_EPSILON);
 
     *r1P = sqrt(-2.0*log(u1)) * cos(2.0*M_PI*u2);
     *r2P = sqrt(-2.0*log(u1)) * sin(2.0*M_PI*u2);
 }
+#endif
 
 
 
@@ -332,7 +366,8 @@ precomputeImageWeights(struct ProgramState * const stateP,
             double r[2];
             unsigned int k;
 
-            randomNormal2(&r[0], &r[1]);
+            pm_gaussrand2(&stateP->randSt, &r[0], &r[1]);
+
             for (k = 0; k < 2; ++k) {
                 int const img =
                     r[k] * sigma + pctGray * stateP->inFileCt * 0.999999;
@@ -453,26 +488,6 @@ blendImages(BlendType             const blend,
 
 
 
-static void
-termState(struct ProgramState * const stateP) {
-/*----------------------------------------------------------------------------
-  Deallocate all of the resources we allocated.
------------------------------------------------------------------------------*/
-    unsigned int i;
-
-    for (i = 0; i < stateP->inFileCt; ++i) {
-        pnm_freepamrow(stateP->inTupleRows[i]);
-        pm_close(stateP->inPam[i].file);
-    }
-
-    free(stateP->outTupleRow);
-    free(stateP->inTupleRows);
-    free(stateP->inPam);
-    pm_close(stateP->outPam.file);
-}
-
-
-
 int
 main(int argc, const char * argv[]) {
 
@@ -483,13 +498,14 @@ main(int argc, const char * argv[]) {
 
     parseCommandLine(argc, argv, &cmdline);
 
-    srand(cmdline.randomseedSpec ? cmdline.randomseed : pm_randseed());
-
-    openInputFiles(cmdline.inFileNameCt, cmdline.inFileName, &state);
+    initInput(cmdline.inFileNameCt, cmdline.inFileName, &state);
 
     if (cmdline.blend == BLEND_MASK)
         initMask(cmdline.maskfile, &state);
 
+    pm_randinit(&state.randSt);
+    pm_srand2(&state.randSt, cmdline.randomseedSpec, cmdline.randomseed);
+
     initOutput(stdout, &state);
 
     if (cmdline.blend == BLEND_MASK)
@@ -497,10 +513,14 @@ main(int argc, const char * argv[]) {
 
     blendImages(cmdline.blend, &state);
 
+    termOutput(&state);
+
+    pm_randterm(&state.randSt);
+
     if (cmdline.blend == BLEND_MASK)
         termMask(&state);
 
-    termState(&state);
+    termInput(&state);
 
     freeCmdline(&cmdline);
 
diff --git a/editor/pamrecolor.c b/editor/pamrecolor.c
index 8c5bce12..86c1965c 100644
--- a/editor/pamrecolor.c
+++ b/editor/pamrecolor.c
@@ -1,3 +1,4 @@
+
 /* ----------------------------------------------------------------------
  *
  * Replace every pixel in an image with one of equal luminance
@@ -32,6 +33,7 @@
 
 #include "mallocvar.h"
 #include "nstring.h"
+#include "rand.h"
 #include "shhopt.h"
 #include "pam.h"
 
@@ -42,7 +44,7 @@
 #define CLAMPxy(N, A, B) MAX(MIN((float)(N), (float)(B)), (float)(A))
 
 
-struct rgbfrac {
+struct Rgbfrac {
     /* This structure represents red, green, and blue, each expressed
        as a fraction from 0.0 to 1.0.
     */
@@ -51,19 +53,19 @@ struct rgbfrac {
     float bfrac;
 };
 
-struct cmdlineInfo {
+struct CmdlineInfo {
     /* This structure represents all of the information the user
        supplied in the command line but in a form that's easy for the
        program to use.
     */
     const char *    inputFileName;  /* '-' if stdin */
     const char *    colorfile;      /* NULL if unspecified */
-    struct rgbfrac  color2gray;
+    struct Rgbfrac  color2gray;
         /* colorspace/rmult/gmult/bmult options.  Negative numbers if
            unspecified.
         */
     unsigned int    targetcolorSpec;
-    struct rgbfrac  targetcolor;
+    struct Rgbfrac  targetcolor;
     unsigned int    randomseed;
     unsigned int    randomseedSpec;
 };
@@ -71,7 +73,7 @@ struct cmdlineInfo {
 
 
 static float
-rgb2gray(struct rgbfrac * const color2grayP,
+rgb2gray(struct Rgbfrac * const color2grayP,
          float            const red,
          float            const grn,
          float            const blu) {
@@ -122,7 +124,7 @@ getColorRow(struct pam  * const pamP,
 
 static void
 convertRowToGray(struct pam     * const pamP,
-                 struct rgbfrac * const color2gray,
+                 struct Rgbfrac * const color2gray,
                  tuplen         * const tupleRow,
                  samplen        * const grayRow) {
 /*----------------------------------------------------------------------
@@ -160,7 +162,7 @@ convertRowToGray(struct pam     * const pamP,
 static void
 explicitlyColorRow(struct pam *   const pamP,
                    tuplen *       const rowData,
-                   struct rgbfrac const tint) {
+                   struct Rgbfrac const tint) {
 
     unsigned int col;
 
@@ -175,17 +177,25 @@ explicitlyColorRow(struct pam *   const pamP,
 
 static void
 randomlyColorRow(struct pam *   const pamP,
-                 tuplen *       const rowData) {
+                 tuplen *       const rowData,
+                 bool           const randomseedSpec,
+                 unsigned int   const randomseed) {
 /*----------------------------------------------------------------------
   Assign each tuple in a row a random color.
 ------------------------------------------------------------------------*/
     unsigned int col;
+    struct pm_randSt randSt;
+
+    pm_randinit(&randSt);
+    pm_srand2(&randSt, randomseedSpec, randomseed);
 
     for (col = 0; col < pamP->width; ++col) {
-        rowData[col][PAM_RED_PLANE] = rand() / (float)RAND_MAX;
-        rowData[col][PAM_GRN_PLANE] = rand() / (float)RAND_MAX;
-        rowData[col][PAM_BLU_PLANE] = rand() / (float)RAND_MAX;
+        rowData[col][PAM_RED_PLANE] = pm_drand(&randSt);
+        rowData[col][PAM_GRN_PLANE] = pm_drand(&randSt);
+        rowData[col][PAM_BLU_PLANE] = pm_drand(&randSt);
     }
+
+    pm_randterm(&randSt);
 }
 
 
@@ -193,7 +203,7 @@ randomlyColorRow(struct pam *   const pamP,
 static void
 recolorRow(struct pam     * const inPamP,
            tuplen         * const inRow,
-           struct rgbfrac * const color2grayP,
+           struct Rgbfrac * const color2grayP,
            tuplen         * const colorRow,
            struct pam     * const outPamP,
            tuplen         * const outRow) {
@@ -293,10 +303,10 @@ recolorRow(struct pam     * const inPamP,
 
 
 
-static struct rgbfrac
+static struct Rgbfrac
 color2GrayFromCsName(const char * const csName) {
 
-    struct rgbfrac retval;
+    struct Rgbfrac retval;
 
     /* Thanks to
        http://www.brucelindbloom.com/index.html?Eqn_RGB_XYZ_Matrix.html
@@ -352,7 +362,7 @@ color2GrayFromCsName(const char * const csName) {
 
 static void
 parseCommandLine(int argc, const char ** const argv,
-                 struct cmdlineInfo * const cmdlineP ) {
+                 struct CmdlineInfo * const cmdlineP ) {
 
     optEntry     * option_def;
         /* Instructions to OptParseOptions3 on how to parse our options */
@@ -445,7 +455,7 @@ parseCommandLine(int argc, const char ** const argv,
 
 int
 main(int argc, const char *argv[]) {
-    struct cmdlineInfo cmdline;          /* Command-line parameters */
+    struct CmdlineInfo cmdline;          /* Command-line parameters */
     struct pam         inPam;
     struct pam         outPam;
     struct pam         colorPam;
@@ -462,8 +472,6 @@ main(int argc, const char *argv[]) {
 
     parseCommandLine(argc, argv, &cmdline);
 
-    srand(cmdline.randomseedSpec ? cmdline.randomseed : pm_randseed());
-
     ifP = pm_openr(cmdline.inputFileName);
     inPam.comment_p = &comments;
     pnm_readpaminit(ifP, &inPam, PAM_STRUCT_SIZE(comment_p));
@@ -474,7 +482,6 @@ main(int argc, const char *argv[]) {
     outPam.depth = 4 - (inPam.depth % 2);
     outPam.allocation_depth = outPam.depth;
     strcpy(outPam.tuple_type, PAM_PPM_TUPLETYPE);
-    pnm_writepaminit(&outPam);
 
     if (cmdline.colorfile) {
         colorfP = pm_openr(cmdline.colorfile);
@@ -492,6 +499,8 @@ main(int argc, const char *argv[]) {
 
     colorRowBuffer = pnm_allocpamrown(&outPam);
 
+    pnm_writepaminit(&outPam);
+
     for (row = 0; row < inPam.height; ++row) {
         tuplen * colorRow;
 
@@ -505,7 +514,8 @@ main(int argc, const char *argv[]) {
             if (cmdline.targetcolorSpec)
                 explicitlyColorRow(&colorPam, colorRow, cmdline.targetcolor);
             else
-                randomlyColorRow(&colorPam, colorRow);
+                randomlyColorRow(&colorPam, colorRow,
+                                 cmdline.randomseedSpec, cmdline.randomseed);
         }
         recolorRow(&inPam, inRow,
                    &cmdline.color2gray, colorRow,
diff --git a/editor/pamrubber.c b/editor/pamrubber.c
index 602701ec..fda31203 100644
--- a/editor/pamrubber.c
+++ b/editor/pamrubber.c
@@ -1,20 +1,18 @@
-/*----------------------------------------------------------------------------*/
-
-/* pamrubber.c - transform images using Rubber Sheeting algorithm
-**               see: http://www.schaik.com/netpbm/rubber/
-**
-** Copyright (C) 2011 by Willem van Schaik (willem@schaik.com)
-**
-** Permission to use, copy, modify, and distribute this software and its
-** documentation for any purpose and without fee is hereby granted, provided
-** that the above copyright notice appear in all copies and that both that
-** copyright notice and this permission notice appear in supporting
-** documentation.  This software is provided "as is" without express or
-** implied warranty.
-*/
-
-/*----------------------------------------------------------------------------*/
-
+/*=============================================================================
+                              pamrubber
+===============================================================================
+  Transform images using Rubber Sheeting algorithm
+  See: http://www.schaik.com/netpbm/rubber/
+
+  Copyright (C) 2011 by Willem van Schaik (willem@schaik.com)
+
+  Permission to use, copy, modify, and distribute this software and its
+  documentation for any purpose and without fee is hereby granted, provided
+  that the above copyright notice appear in all copies and that both that
+  copyright notice and this permission notice appear in supporting
+  documentation.  This software is provided "as is" without express or
+  implied warranty.
+=============================================================================*/
 #include <stdio.h>
 #include <stdlib.h>
 #include <unistd.h>
@@ -25,12 +23,12 @@
 
 #include "pm_c_util.h"
 #include "mallocvar.h"
+#include "rand.h"
 #include "shhopt.h"
 #include "pam.h"
 #include "pamdraw.h"
 
 
-
 typedef struct {
   double x;
   double y;
@@ -54,7 +52,7 @@ typedef struct {
     point br;  /* bottom right */
 } quadrilateral;
 
-struct cmdlineInfo {
+struct CmdlineInfo {
     unsigned int nCP;
     point        oldCP[4];
     point        newCP[4];
@@ -64,14 +62,14 @@ struct cmdlineInfo {
     unsigned int frame;
     unsigned int linear;
     unsigned int verbose;
-    unsigned int randseedSpec;
-    unsigned int randseed;
+    unsigned int randomseedSpec;
+    unsigned int randomseed;
 };
 
 
 static void
 parseCmdline(int argc, const char ** argv,
-             struct cmdlineInfo * const cmdlineP) {
+             struct CmdlineInfo * const cmdlineP) {
 
 /* parse all parameters from the command line */
 
@@ -92,10 +90,10 @@ parseCmdline(int argc, const char ** argv,
     OPTENT3(0, "frame",    OPT_FLAG, NULL, &cmdlineP->frame,    0);
     OPTENT3(0, "linear",   OPT_FLAG, NULL, &cmdlineP->linear,   0);
     OPTENT3(0, "verbose",  OPT_FLAG, NULL, &cmdlineP->verbose,  0);
-    OPTENT3(0, "randseed", OPT_UINT, &cmdlineP->randseed,
-            &cmdlineP->randseedSpec, 0);
-    OPTENT3(0, "randomseed", OPT_UINT, &cmdlineP->randseed,
-            &cmdlineP->randseedSpec, 0);
+    OPTENT3(0, "randseed", OPT_UINT, &cmdlineP->randomseed,
+            &cmdlineP->randomseedSpec, 0);
+    OPTENT3(0, "randomseed", OPT_UINT, &cmdlineP->randomseed,
+            &cmdlineP->randomseedSpec, 0);
 
     opt.opt_table = option_def;
     opt.short_allowed = FALSE;  /* we have no short (old-fashioned) options */
@@ -339,28 +337,29 @@ windtriangle(triangle * const tP,
 
 
 static double
-tiny(void) {
+tiny(struct pm_randSt * const randStP) {
 
-    if (rand() % 2)
-        return +1E-6 * (double) ((rand() % 90) + 9);
+    if (pm_rand(randStP) % 2)
+        return +1E-6 * (double) ((pm_rand(randStP) % 90) + 9);
     else
-        return -1E-6 * (double) ((rand() % 90) + 9);
+        return -1E-6 * (double) ((pm_rand(randStP) % 90) + 9);
 }
 
 
 
 static void
 angle(point * const p1P,
-      point * const p2P) {
+      point * const p2P,
+      struct pm_randSt * const randStP) {
 /*----------------------------------------------------------------------------
    Move *p2P slightly if necessary to make sure the line (*p1P, *p2P)
    is not horizontal or vertical.
 -----------------------------------------------------------------------------*/
     if (p1P->x == p2P->x) { /* vertical line */
-        p2P->x += tiny();
+        p2P->x += tiny(randStP);
     }
     if (p1P->y == p2P->y) { /* horizontal line */
-        p2P->y += tiny();
+        p2P->y += tiny(randStP);
     }
 }
 
@@ -720,8 +719,10 @@ static void drawClippedTriangle(const struct pam * const pamP,
 
 
 static void
-prepTrig(int const wd,
-         int const ht) {
+prepTrig(int          const wd,
+         int          const ht,
+         bool         const randomseedSpec,
+         unsigned int const randomseed) {
 
 /* create triangles using control points */
 
@@ -731,16 +732,20 @@ prepTrig(int const wd,
     point c2p1, c2p2, c2p3, c2p4;
     line l1, l2;
     point p0;
+    struct pm_randSt randSt;
+
+    pm_randinit(&randSt);
+    pm_srand2(&randSt, randomseedSpec, randomseed);
 
-    rtl1 = makepoint(0.0 + tiny(),               0.0 + tiny());
-    rtr1 = makepoint((double) wd - 1.0 + tiny(), 0.0 + tiny());
-    rbl1 = makepoint(0.0 + tiny(),               (double) ht - 1.0 + tiny());
-    rbr1 = makepoint((double) wd - 1.0 + tiny(), (double) ht - 1.0 + tiny());
+    rtl1 = makepoint(0.0 + tiny(&randSt),               0.0 + tiny(&randSt));
+    rtr1 = makepoint((double) wd - 1.0 + tiny(&randSt), 0.0 + tiny(&randSt));
+    rbl1 = makepoint(0.0 + tiny(&randSt),               (double) ht - 1.0 + tiny(&randSt));
+    rbr1 = makepoint((double) wd - 1.0 + tiny(&randSt), (double) ht - 1.0 + tiny(&randSt));
 
-    rtl2 = makepoint(0.0 + tiny(),               0.0 + tiny());
-    rtr2 = makepoint((double) wd - 1.0 + tiny(), 0.0 + tiny());
-    rbl2 = makepoint(0.0 + tiny(),               (double) ht - 1.0 + tiny());
-    rbr2 = makepoint((double) wd - 1.0 + tiny(), (double) ht - 1.0 + tiny());
+    rtl2 = makepoint(0.0 + tiny(&randSt),               0.0 + tiny(&randSt));
+    rtr2 = makepoint((double) wd - 1.0 + tiny(&randSt), 0.0 + tiny(&randSt));
+    rbl2 = makepoint(0.0 + tiny(&randSt),               (double) ht - 1.0 + tiny(&randSt));
+    rbr2 = makepoint((double) wd - 1.0 + tiny(&randSt), (double) ht - 1.0 + tiny(&randSt));
 
     if (nCP == 1) {
         c1p1 = oldCP[0];
@@ -772,8 +777,8 @@ prepTrig(int const wd,
         c2p2 = newCP[1];
 
         /* check for hor/ver edges */
-        angle (&c1p1, &c1p2);
-        angle (&c2p1, &c2p2);
+        angle (&c1p1, &c1p2, &randSt);
+        angle (&c2p1, &c2p2, &randSt);
 
         /* connect two control points to corners to get 6 triangles */
         /* left side */
@@ -811,13 +816,13 @@ prepTrig(int const wd,
         /* Move vertices slightly if necessary to make sure no edge is
            horizontal or vertical.
         */
-        angle(&c1p1, &c1p2);
-        angle(&c1p2, &c1p3);
-        angle(&c1p3, &c1p1);
+        angle(&c1p1, &c1p2, &randSt);
+        angle(&c1p2, &c1p3, &randSt);
+        angle(&c1p3, &c1p1, &randSt);
 
-        angle(&c2p1, &c2p2);
-        angle(&c2p2, &c2p3);
-        angle(&c2p3, &c2p1);
+        angle(&c2p1, &c2p2, &randSt);
+        angle(&c2p2, &c2p3, &randSt);
+        angle(&c2p3, &c2p1, &randSt);
 
         if (windtriangle(&tri1s[0], c1p1, c1p2, c1p3)) {
             tri2s[0] = maketriangle(c2p1, c2p2, c2p3);
@@ -871,19 +876,19 @@ prepTrig(int const wd,
         c2p4 = newCP[3];
 
         /* check for hor/ver edges */
-        angle (&c1p1, &c1p2);
-        angle (&c1p2, &c1p3);
-        angle (&c1p3, &c1p4);
-        angle (&c1p4, &c1p1);
-        angle (&c1p1, &c1p3);
-        angle (&c1p2, &c1p4);
-
-        angle (&c2p1, &c2p2);
-        angle (&c2p2, &c2p3);
-        angle (&c2p3, &c2p4);
-        angle (&c2p4, &c2p1);
-        angle (&c2p1, &c2p3);
-        angle (&c2p2, &c2p4);
+        angle (&c1p1, &c1p2, &randSt);
+        angle (&c1p2, &c1p3, &randSt);
+        angle (&c1p3, &c1p4, &randSt);
+        angle (&c1p4, &c1p1, &randSt);
+        angle (&c1p1, &c1p3, &randSt);
+        angle (&c1p2, &c1p4, &randSt);
+
+        angle (&c2p1, &c2p2, &randSt);
+        angle (&c2p2, &c2p3, &randSt);
+        angle (&c2p3, &c2p4, &randSt);
+        angle (&c2p4, &c2p1, &randSt);
+        angle (&c2p1, &c2p3, &randSt);
+        angle (&c2p2, &c2p4, &randSt);
 
         /*-------------------------------------------------------------------*/
         /*        -1-      -2-        -3-      -4-        -5-      -6-       */
@@ -979,6 +984,8 @@ prepTrig(int const wd,
                      &tri2s[9], c2p3, c2p1, rtl2, rtr2, rbl2, rbr2);
         nTri = 10;
     }
+
+    pm_randterm(&randSt);
 }
 
 
@@ -1277,7 +1284,7 @@ warpQuad(point   const p2,
 
 
 static void
-setGlobalCP(struct cmdlineInfo const cmdline) {
+setGlobalCP(struct CmdlineInfo const cmdline) {
 
     unsigned int i;
 
@@ -1392,7 +1399,7 @@ pix(tuple **     const tuples,
 int
 main(int argc, const char ** const argv) {
 
-    struct cmdlineInfo cmdline;
+    struct CmdlineInfo cmdline;
     FILE * ifP;
     struct pam inpam, outpam;
     tuple ** inTuples;
@@ -1405,8 +1412,6 @@ main(int argc, const char ** const argv) {
 
     setGlobalCP(cmdline);
 
-    srand(cmdline.randseedSpec ? cmdline.randseed : pm_randseed());
-
     ifP = pm_openr(cmdline.fileName);
 
     inTuples = pnm_readpam(ifP, &inpam, PAM_STRUCT_SIZE(tuple_type));
@@ -1421,7 +1426,8 @@ main(int argc, const char ** const argv) {
     makeAllWhite(&outpam, outTuples);
 
     if (cmdline.tri)
-        prepTrig(inpam.width, inpam.height);
+        prepTrig(inpam.width, inpam.height,
+                 cmdline.randomseedSpec, cmdline.randomseed);
     if (cmdline.quad)
         prepQuad();
 
diff --git a/editor/pbmreduce.c b/editor/pbmreduce.c
index 3a0968fe..70caa581 100644
--- a/editor/pbmreduce.c
+++ b/editor/pbmreduce.c
@@ -11,9 +11,11 @@
 */
 
 #include "pm_c_util.h"
-#include "pbm.h"
 #include "mallocvar.h"
+#include "rand.h"
 #include "shhopt.h"
+#include "pbm.h"
+
 #include <assert.h>
 
 #define SCALE 1024
@@ -105,7 +107,7 @@ parseCommandLine(int argc, const char ** argv,
         unsigned int scale;
 
         scale = strtol(argv[1], &endptr, 10);
-        if (*argv[1] == '\0') 
+        if (*argv[1] == '\0')
             pm_error("Scale argument is a null string.  Must be a number.");
         else if (*endptr != '\0')
             pm_error("Scale argument contains non-numeric character '%c'.",
@@ -115,7 +117,7 @@ parseCommandLine(int argc, const char ** argv,
                      "You specified %d", scale);
         else if (scale > INT_MAX / scale)
             pm_error("Scale argument too large.  You specified %d", scale);
-        else 
+        else
             cmdlineP->scale = scale;
 
         if (argc-1 > 1) {
@@ -145,10 +147,11 @@ struct FS {
 static void
 initializeFloydSteinberg(struct FS  * const fsP,
                          int          const newcols,
-                         unsigned int const seed,
-                         bool         const seedSpec) {
+                         bool         const seedSpec,
+                         unsigned int const seed) {
 
     unsigned int col;
+    struct pm_randSt randSt;
 
     MALLOCARRAY(fsP->thiserr, newcols + 2);
     MALLOCARRAY(fsP->nexterr, newcols + 2);
@@ -156,33 +159,36 @@ initializeFloydSteinberg(struct FS  * const fsP,
     if (fsP->thiserr == NULL || fsP->nexterr == NULL)
         pm_error("out of memory");
 
-    srand(seedSpec ? seed : pm_randseed());
+    pm_randinit(&randSt);
+    pm_srand2(&randSt, seedSpec, seed);
 
     for (col = 0; col < newcols + 2; ++col)
-        fsP->thiserr[col] = (rand() % SCALE - HALFSCALE) / 4;
+        fsP->thiserr[col] = ((int) (pm_rand(&randSt) % SCALE) - HALFSCALE) / 4;
         /* (random errors in [-SCALE/8 .. SCALE/8]) */
+
+    pm_randterm(&randSt);
 }
 
 
 
 /*
     Scanning method
-    
+
     In Floyd-Steinberg dithering mode horizontal direction of scan alternates
     between rows; this is called "serpentine scanning".
-    
+
     Example input (14 x 7), N=3:
-    
+
     111222333444xx    Fractional pixels on the right edge and bottom edge (x)
-    111222333444xx    are ignored; their values do not influence output. 
+    111222333444xx    are ignored; their values do not influence output.
     111222333444xx
     888777666555xx
     888777666555xx
     888777666555xx
     xxxxxxxxxxxxxx
-    
+
     Output (4 x 2):
-    
+
     1234
     8765
 
@@ -196,11 +202,14 @@ enum Direction { RIGHT_TO_LEFT, LEFT_TO_RIGHT };
 static enum Direction
 oppositeDir(enum Direction const arg) {
 
+    enum Direction retval;
+
     switch (arg) {
-    case LEFT_TO_RIGHT: return RIGHT_TO_LEFT;
-    case RIGHT_TO_LEFT: return LEFT_TO_RIGHT;
+    case LEFT_TO_RIGHT: retval = RIGHT_TO_LEFT; break;
+    case RIGHT_TO_LEFT: retval = LEFT_TO_RIGHT; break;
     }
-    assert(false);  /* All cases handled above */
+
+    return retval;
 }
 
 
@@ -240,7 +249,7 @@ main(int argc, const char * argv[]) {
 
     if (cmdline.halftone == QT_FS)
         initializeFloydSteinberg(&fs, newcols,
-                                 cmdline.randomseed, cmdline.randomseedSpec);
+                                 cmdline.randomseedSpec, cmdline.randomseed);
     else {
         /* These variables are meaningless in this case, and the values
            should never be used.
@@ -258,10 +267,10 @@ main(int argc, const char * argv[]) {
         int limitCol;
         int startCol;
         int step;
-   
+
         for (colChar = 0; colChar < colChars; ++colChar)
             newbitrow[colChar] = 0x00;  /* Clear to white */
- 
+
         for (subrow = 0; subrow < cmdline.scale; ++subrow)
             pbm_readpbmrow(ifP, bitslice[subrow], cols, format);
 
@@ -274,7 +283,7 @@ main(int argc, const char * argv[]) {
         case LEFT_TO_RIGHT: {
             startCol = 0;
             limitCol = newcols;
-            step = +1;  
+            step = +1;
         } break;
         case RIGHT_TO_LEFT: {
             startCol = newcols - 1;
diff --git a/editor/pnmflip b/editor/pnmflip
index 07d4ddb9..962198a2 100755
--- a/editor/pnmflip
+++ b/editor/pnmflip
@@ -46,6 +46,13 @@ exec perl -w -x -S -- "$0" "$@"
 use strict;
 use File::Basename;
 use Cwd 'abs_path';
+use IO::Handle;
+
+sub pm_message($) {
+    STDERR->print("pnmflip: $_[0]\n");
+}
+
+
 
 my $xformOpt;
 my @miscOptions;
@@ -78,8 +85,7 @@ foreach (@ARGV) {
     } else {
         # It's a parameter
         if (defined($infile)) {
-            print(STDERR
-                  "You may specify at most one non-option parameter.\n");
+            pm_message("You may specify at most one non-option parameter.");
             exit(10);
         } else {
             $infile = $_;
diff --git a/editor/pnmquant b/editor/pnmquant
index f7af9e7a..0bebce69 100755
--- a/editor/pnmquant
+++ b/editor/pnmquant
@@ -37,9 +37,16 @@ use Getopt::Long;
 use File::Spec;
 #use Fcntl ":seek";  # not available in Perl 5.00503
 use Fcntl;  # gets open flags
+use IO::Handle;
 
 my ($TRUE, $FALSE) = (1,0);
 
+sub pm_message($) {
+    STDERR->print("pnmquant: $_[0]\n");
+}
+
+
+
 my ($SEEK_SET, $SEEK_CUR, $SEEK_END) = (0, 1, 2);
 
 
@@ -102,17 +109,16 @@ sub parseCommandLine(@) {
                                   "plain");
 
     if (!$optsAreValid) {
-        print(STDERR "Invalid option syntax.\n");
+        pm_message("Invalid option syntax");
         exit(1);
     }
     if (@ARGV > 2) {
-        print(STDERR "This program takes at most 2 arguments.  You specified ",
-              scalar(@ARGV), "\n");
+        pm_message("This program takes at most 2 arguments.  You specified " .
+                   scalar(@ARGV));
         exit(1);
     } 
     elsif (@ARGV < 1) {
-        print(STDERR 
-              "You must specify the number of colors as an argument.\n");
+        pm_message("You must specify the number of colors as an argument.");
         exit(1);
     }
     my $infile;
@@ -120,9 +126,8 @@ sub parseCommandLine(@) {
     
     if (!($cmdline{ncolors} =~ m{ ^[[:digit:]]+$ }x ) || 
         $cmdline{ncolors} == 0) {
-        print(STDERR 
-              "Number of colors argument '$cmdline{ncolors}' " .
-              "is not a positive integer.\n");
+        pm_message("Number of colors argument '$cmdline{ncolors}' " .
+                   "is not a positive integer.");
         exit(1);
     }
 
@@ -210,8 +215,8 @@ sub makeColormap($$$$$$$) {
     my ($mapfileFh, $mapfileSpec) = tempFile(".pnm");
 
     if (!defined($mapfileFh)) {
-        print(STDERR "Unable to create temporary file for colormap.  " .
-              "errno = $ERRNO\n");
+        pm_message("Unable to create temporary file for colormap.  " .
+                   "errno = $ERRNO");
         exit(1);
     }
        
@@ -223,8 +228,8 @@ sub makeColormap($$$$$$$) {
         (defined($opt_center)    ? 1 : 0);
 
     if ($colorSummaryOptCt > 1) {
-        print(STDERR "You can specify only one of " .
-              "-meanpixel, -meancolor, and -center\n");
+        pm_message("You can specify only one of " .
+                   "-meanpixel, -meancolor, and -center");
         exit(1);
     }
     if (defined($opt_meanpixel)) {
@@ -240,8 +245,8 @@ sub makeColormap($$$$$$$) {
         (defined($opt_spreadbrightness) ? 1 : 0);
 
     if ($spreadOptCt > 1) {
-        print(STDERR "You can specify only one of " .
-              "-spreadluminosity and -spreadbrightness\n");
+        pm_message("You can specify only one of " .
+                   "-spreadluminosity and -spreadbrightness");
         exit(1);
     }
 
@@ -263,7 +268,7 @@ sub makeColormap($$$$$$$) {
     my $maprc = system("pnmcolormap", $ncolors, @options);
 
     if ($maprc != 0) {
-        print(STDERR "pnmcolormap failed, rc=$maprc\n");
+        pm_message("pnmcolormap failed, rc=$maprc");
         exit(1);
     } 
     return $mapfileSpec;
@@ -287,15 +292,15 @@ sub remap($$$$$$) {
     }
     if ($opt_norandom) {
         if (defined($opt_randomseed)) {
-             print(STDERR "You cannot specify -randomseed with -norandom\n");
+             pm_message("You cannot specify -randomseed with -norandom");
              exit(1);
         }
         push(@options, "-norandom");
     }
     if (defined($opt_randomseed)) {
         if ($opt_randomseed < 0) {
-             print(STDERR "-randomseed value must not be negative.  " .
-                   "You specified $opt_randomseed\n");
+             pm_message("-randomseed value must not be negative.  " .
+                        "You specified $opt_randomseed");
              exit(10);
         }
         push(@options, "-randomseed=$opt_randomseed");
@@ -310,7 +315,7 @@ sub remap($$$$$$) {
     my $remaprc = system("pnmremap", "-mapfile=$mapfileSpec", @options);
     
     if ($remaprc != 0) {
-        print(STDERR "pnmremap failed, rc=$remaprc\n");
+        pm_message("pnmremap failed, rc=$remaprc");
         exit(1);
     }
 }
diff --git a/editor/pnmquantall b/editor/pnmquantall
index aea6cc84..80d00fce 100755
--- a/editor/pnmquantall
+++ b/editor/pnmquantall
@@ -57,11 +57,18 @@ use warnings;
 use English;
 use Fcntl;  # gets open flags
 use File::Copy;
+use IO::Handle;
 
 my $TRUE=1; my $FALSE = 0;
 
 
 
+sub pm_message($) {
+    STDERR->print("pnmquantall: $_[0]\n");
+}
+
+
+
 sub doVersionHack($) {
     my ($argvR) = @_;
 
@@ -84,7 +91,7 @@ sub parseArgs($$$$) {
 
     if (@argv > 0 && $argv[0] eq "-ext") {
         if (@argv < 2) {
-            print STDERR ("-ext requires a value\n");
+            pm_message("-ext requires a value");
         exit(100);
         } else {
             $$extR = $argv[1];
@@ -96,8 +103,8 @@ sub parseArgs($$$$) {
     }
 
     if (@argv < $firstArgPos + 2) {
-        print STDERR ("Not enough arguments.  You need at least the number " .
-                      "of colors and one file name\n");
+        pm_message("Not enough arguments.  You need at least the number " .
+                      "of colors and one file name");
         exit(100);
     }
     
@@ -212,7 +219,7 @@ if (!$progError) {
 my $exitStatus;
 
 if ($progError) {
-    print STDERR ("Failed.  $progError\n");
+    pm_message("Failed.  $progError");
     $exitStatus = 1;
 } else {
     $exitStatus = 0;
diff --git a/editor/pnmremap.c b/editor/pnmremap.c
index 0c0096ba..e5b59d04 100644
--- a/editor/pnmremap.c
+++ b/editor/pnmremap.c
@@ -28,6 +28,7 @@
 #include "pm_c_util.h"
 #include "mallocvar.h"
 #include "nstring.h"
+#include "rand.h"
 #include "shhopt.h"
 #include "pam.h"
 #include "ppm.h"
@@ -399,20 +400,21 @@ randomizeError(long **       const err,
    Set a random error in the range [-1 .. 1] (normalized via FS_SCALE)
    in the error array err[][].
 -----------------------------------------------------------------------------*/
-    unsigned int const seed = (random.init == RANDOM_WITHSEED) ?
-        random.seed : pm_randseed();
-
     unsigned int col;
+    struct pm_randSt randSt;
 
     assert(random.init != RANDOM_NONE);
 
-    srand(seed);
+    pm_randinit(&randSt);
+    pm_srand2(&randSt, random.init == RANDOM_WITHSEED, random.seed);
 
     for (col = 0; col < width; ++col) {
         unsigned int plane;
         for (plane = 0; plane < depth; ++plane)
-            err[plane][col] = rand() % (FS_SCALE * 2) - FS_SCALE;
+            err[plane][col] = pm_rand(&randSt) % (FS_SCALE * 2) - FS_SCALE;
     }
+
+    pm_randterm(&randSt);
 }
 
 
diff --git a/editor/specialty/pampaintspill.c b/editor/specialty/pampaintspill.c
index c7994723..5cd482d5 100644
--- a/editor/specialty/pampaintspill.c
+++ b/editor/specialty/pampaintspill.c
@@ -45,11 +45,11 @@
 
 #include "mallocvar.h"
 #include "nstring.h"
+#include "rand.h"
 #include "shhopt.h"
 #include "pam.h"
 #include "pammap.h"
 
-
 static time_t const timeUpdateDelta = 30;
     /* Seconds between progress updates */
 static int const    minUpdates = 4;
@@ -67,6 +67,8 @@ struct cmdlineInfo {
     unsigned int all;
     float        power;
     unsigned int downsample;
+    unsigned int randomseedSpec;
+    unsigned int randomseed;
 };
 
 struct coords {
@@ -98,16 +100,18 @@ parseCommandLine(int argc, const char ** const argv,
     MALLOCARRAY_NOFAIL(option_def, 100);
     option_def_index = 0;          /* Incremented by OPTENTRY */
 
-    OPTENT3(0, "bgcolor",    OPT_STRING, &cmdlineP->bgcolor,    
+    OPTENT3(0, "bgcolor",    OPT_STRING, &cmdlineP->bgcolor,
             &bgcolorSpec, 0);
     OPTENT3(0, "wrap",       OPT_FLAG,   NULL,
             &cmdlineP->wrap,       0);
     OPTENT3(0, "all",        OPT_FLAG,   NULL,
             &cmdlineP->all,        0);
-    OPTENT3(0, "power",      OPT_FLOAT,  &cmdlineP->power,      
+    OPTENT3(0, "power",      OPT_FLOAT,  &cmdlineP->power,
             &powerSpec, 0);
-    OPTENT3(0, "downsample", OPT_UINT,   &cmdlineP->downsample, 
+    OPTENT3(0, "downsample", OPT_UINT,   &cmdlineP->downsample,
             &downsampleSpec, 0);
+    OPTENT3(0, "randomseed", OPT_UINT,   &cmdlineP->randomseed,
+            &cmdlineP->randomseedSpec, 0);
 
     opt.opt_table = option_def;
     opt.short_allowed = 0;
@@ -223,7 +227,9 @@ locatePaintSources(struct pam *            const pamP,
                    tuple **                const tuples,
                    tuple                   const bgColor,
                    unsigned int            const downsample,
-                   struct paintSourceSet * const paintSourcesP) {
+                   struct paintSourceSet * const paintSourcesP,
+                   bool                    const randomseedSpec,
+                   unsigned int            const randomseed) {
 /*--------------------------------------------------------------------
   Construct a list of all pixel coordinates in the input image that
   represent a non-background color.
@@ -248,21 +254,24 @@ locatePaintSources(struct pam *            const pamP,
     pm_message("Image contains %u background + %u non-background pixels",
                pamP->width * pamP->height - paintSources.size,
                paintSources.size);
-    
+
     /* Reduce the number of paint sources to reduce execution time. */
     if (downsample > 0 && downsample < paintSources.size) {
+        struct pm_randSt randSt;
         unsigned int i;
 
-        srand(pm_randseed());
+        pm_randinit(&randSt);
+        pm_srand2(&randSt, randomseedSpec, randomseed);
 
         for (i = 0; i < downsample; ++i) {
             unsigned int const swapIdx =
-                i + rand() % (paintSources.size - i);
+                i + pm_rand(&randSt) % (paintSources.size - i);
             struct coords const swapVal = paintSources.list[i];
 
             paintSources.list[i] = paintSources.list[swapIdx];
             paintSources.list[swapIdx] = swapVal;
         }
+        pm_randterm(&randSt);
         paintSources.size = downsample;
     }
 
@@ -426,7 +435,7 @@ produceOutputImage(struct pam *          const pamP,
     for (row = 0; row < pamP->height; ++row) {
         struct coords   target;
         double        * newColor;
-        
+
         MALLOCARRAY(newColor, pamP->depth);
 
         target.y = row;
@@ -484,7 +493,8 @@ main(int argc, const char *argv[]) {
                pnm_colorname(&inPam, bgColor, PAM_COLORNAME_HEXOK));
 
     locatePaintSources(&inPam, inTuples, bgColor, cmdline.downsample,
-                       &paintSources);
+                       &paintSources,
+                       cmdline.randomseedSpec, cmdline.randomseed);
 
     produceOutputImage(&inPam, inTuples, bgColor, paintSources, distFunc,
                        cmdline.power, cmdline.all, &outTuples);
@@ -498,3 +508,5 @@ main(int argc, const char *argv[]) {
 
     return 0;
 }
+
+
diff --git a/editor/specialty/ppmshift.c b/editor/specialty/ppmshift.c
index cdb0f173..27cbb78c3 100644
--- a/editor/specialty/ppmshift.c
+++ b/editor/specialty/ppmshift.c
@@ -12,11 +12,11 @@
 #include <stdbool.h>
 
 #include "mallocvar.h"
+#include "rand.h"
 #include "shhopt.h"
 #include "ppm.h"
 
 
-
 struct CmdlineInfo {
     /* All the information the user supplied in the command line,
        in a form easy for the program to use.
@@ -24,6 +24,7 @@ struct CmdlineInfo {
     const char * inputFileName;
 
     unsigned int shift;
+    unsigned int seedSpec;
     unsigned int seed;
 };
 
@@ -42,8 +43,6 @@ parseCommandLine(int argc, const char ** const argv,
 
     unsigned int option_def_index;
 
-    unsigned int seedSpec;
-
     MALLOCARRAY(option_def, 100);
 
     opt.opt_table = option_def;
@@ -52,14 +51,11 @@ parseCommandLine(int argc, const char ** const argv,
 
     option_def_index = 0;   /* incremented by OPTENT3 */
     OPTENT3(0,   "seed",            OPT_UINT,     &cmdlineP->seed,
-            &seedSpec,         0);
+            &cmdlineP->seedSpec,         0);
 
     pm_optParseOptions3(&argc, (char **)argv, opt, sizeof(opt), 0);
         /* Uses and sets argc, argv, and some of *cmdlineP and others. */
 
-    if (!seedSpec)
-        cmdlineP->seed = pm_randseed();
-
     if (argc-1 < 1)
         pm_error("You must specify the shift factor as an argument");
     else {
@@ -84,6 +80,62 @@ parseCommandLine(int argc, const char ** const argv,
 
 
 
+static void
+shiftRow(pixel *            const srcrow,
+         unsigned int       const cols,
+         unsigned int       const shift,
+         pixel *            const destrow,
+         struct pm_randSt * const randStP) {
+
+    /* the range by which a line is shifted lays in the range from */
+    /* -shift/2 .. +shift/2 pixels; however, within this range it is */
+    /* randomly chosen */
+
+    pixel * pP;
+    pixel * pP2;
+    int nowshift;
+
+    if (shift != 0)
+        nowshift = (pm_rand(randStP) % (shift+1)) - ((shift+1) / 2);
+    else
+        nowshift = 0;
+
+    pP  = &srcrow[0];
+    pP2 = &destrow[0];
+
+    /* if the shift value is less than zero, we take the original
+       pixel line and copy it into the destination line translated
+       to the left by x pixels. The empty pixels on the right end
+       of the destination line are filled up with the pixel that
+       is the right-most in the original pixel line.
+    */
+    if (nowshift < 0) {
+        unsigned int col;
+        pP += abs(nowshift);
+        for (col = 0; col < cols; ++col) {
+            PPM_ASSIGN(*pP2, PPM_GETR(*pP), PPM_GETG(*pP), PPM_GETB(*pP));
+            ++pP2;
+            if (col < (cols + nowshift) - 1)
+                ++pP;
+        }
+    } else {
+        unsigned int col;
+        /* The shift value is 0 or positive, so fill the first
+           <nowshift> pixels of the destination line with the
+           first pixel from the source line, and copy the rest of
+           the source line to the dest line
+        */
+        for (col = 0; col < cols; ++col) {
+            PPM_ASSIGN(*pP2, PPM_GETR(*pP), PPM_GETG(*pP), PPM_GETB(*pP));
+            ++pP2;
+            if (col >= nowshift)
+                ++pP;
+        }
+    }
+}
+
+
+
 int
 main(int argc, const char ** argv) {
 
@@ -95,13 +147,15 @@ main(int argc, const char ** argv) {
     pixel * destrow;
     unsigned int row;
     unsigned int shift;
+    struct pm_randSt randSt;
 
     /* parse in 'default' parameters */
     pm_proginit(&argc, argv);
 
     parseCommandLine(argc, argv, &cmdline);
 
-    srand(cmdline.seed);
+    pm_randinit(&randSt);
+    pm_srand2(&randSt, cmdline.seedSpec, cmdline.seed);
 
     ifP = pm_openr(cmdline.inputFileName);
 
@@ -115,67 +169,23 @@ main(int argc, const char ** argv) {
     } else
         shift = cmdline.shift;
 
-    srcrow = ppm_allocrow(cols);
-
+    srcrow  = ppm_allocrow(cols);
     destrow = ppm_allocrow(cols);
 
     ppm_writeppminit(stdout, cols, rows, maxval, 0);
 
-    /** now do the shifting **/
-    /* the range by which a line is shifted lays in the range from */
-    /* -shift/2 .. +shift/2 pixels; however, within this range it is */
-    /* randomly chosen */
     for (row = 0; row < rows; ++row) {
-        pixel * pP;
-        pixel * pP2;
-        unsigned int nowshift;
-
-        if (shift != 0)
-            nowshift = (rand() % (shift+1)) - ((shift+1) / 2);
-        else
-            nowshift = 0;
-
         ppm_readppmrow(ifP, srcrow, cols, maxval, format);
 
-        pP  = &srcrow[0];
-        pP2 = &destrow[0];
-
-        /* if the shift value is less than zero, we take the original
-           pixel line and copy it into the destination line translated
-           to the left by x pixels. The empty pixels on the right end
-           of the destination line are filled up with the pixel that
-           is the right-most in the original pixel line.
-        */
-        if (nowshift < 0) {
-            unsigned int col;
-            pP += abs(nowshift);
-            for (col = 0; col < cols; ++col) {
-                PPM_ASSIGN(*pP2, PPM_GETR(*pP), PPM_GETG(*pP), PPM_GETB(*pP));
-                ++pP2;
-                if (col < (cols + nowshift) - 1)
-                    ++pP;
-            }
-        } else {
-            unsigned int col;
-            /* The shift value is 0 or positive, so fill the first
-               <nowshift> pixels of the destination line with the
-               first pixel from the source line, and copy the rest of
-               the source line to the dest line
-            */
-            for (col = 0; col < cols; ++col) {
-                PPM_ASSIGN(*pP2, PPM_GETR(*pP), PPM_GETG(*pP), PPM_GETB(*pP));
-                ++pP2;
-                if (col >= nowshift)
-                    ++pP;
-            }
-        }
+        shiftRow(srcrow, cols, shift, destrow, &randSt);
 
         ppm_writeppmrow(stdout, destrow, cols, maxval, 0);
     }
 
-    pm_close(ifP);
-    ppm_freerow(srcrow);
     ppm_freerow(destrow);
+    ppm_freerow(srcrow);
+    pm_close(ifP);
+    pm_randterm(&randSt);
 
     return 0;
 }
diff --git a/editor/specialty/ppmspread.c b/editor/specialty/ppmspread.c
index 6753f4fe..7b9558e3 100644
--- a/editor/specialty/ppmspread.c
+++ b/editor/specialty/ppmspread.c
@@ -10,102 +10,158 @@
 
 #include <string.h>
 
+#include "nstring.h"
+#include "rand.h"
+#include "shhopt.h"
 #include "ppm.h"
 
 
+struct CmdlineInfo {
+    /* This structure represents all of the information the user
+       supplied in the command line but in a form that's easy for the
+       program to use.
+    */
+    const char * inputFilename;  /* '-' if stdin */
+    unsigned int spread;
+    unsigned int randomseedSpec;
+    unsigned int randomseed;
+};
+
+
+
+static void
+parseCommandLine(int argc, const char ** const argv,
+                 struct CmdlineInfo * const cmdlineP ) {
+
+    optEntry     * option_def;
+        /* Instructions to OptParseOptions3 on how to parse our options */
+    optStruct3     opt;
+    unsigned int   option_def_index;
+
+    MALLOCARRAY_NOFAIL(option_def, 100);
+    option_def_index = 0;          /* Incremented by OPTENTRY */
+
+    OPTENT3(0, "randomseed", OPT_UINT,   &cmdlineP->randomseed,
+            &cmdlineP->randomseedSpec, 0);
+
+    opt.opt_table = option_def;
+    opt.short_allowed = 0;
+    opt.allowNegNum = 1;
+
+    pm_optParseOptions3( &argc, (char **)argv, opt, sizeof(opt), 0 );
+
+    if (argc-1 < 1)
+        pm_error("You must specify the spread factor as an argument");
+    else {
+        const char * error;
+        pm_string_to_uint(argv[1], &cmdlineP->spread, &error);
+
+        if (error)
+            pm_error("Spread factor '%s' is not an unsigned integer.  %s",
+                     argv[1], error);
+
+        if (argc-1 < 2)
+            cmdlineP->inputFilename = "-";
+        else {
+            cmdlineP->inputFilename = argv[2];
+            if (argc-1 >2)
+                pm_error("Too many arguments: %u.  "
+                         "The only possible arguments are "
+                         "the spread factor and the optional input file name",
+                         argc-1);
+        }
+    }
+}
+
+
+
+static void
+spreadRow(pixel **           const srcarray,
+          unsigned int       const cols,
+          unsigned int       const rows,
+          unsigned int       const spread,
+          unsigned int       const row,
+          pixel **           const destarray,
+          struct pm_randSt * const randStP) {
+
+    unsigned int col;
+
+    for (col = 0; col < cols; ++col) {
+        pixel const p = srcarray[row][col];
+
+        int const xdis = (pm_rand(randStP) % (spread + 1) )
+            - ((spread + 1) / 2);
+        int const ydis = (pm_rand(randStP) % (spread + 1))
+            - ((spread + 1) / 2);
+
+        int const xnew = col + xdis;
+        int const ynew = row + ydis;
+
+        /* only set the displaced pixel if it's within the bounds
+           of the image
+        */
+        if (xnew >= 0 && xnew < cols && ynew >= 0 && ynew < rows) {
+            /* Displacing a pixel is accomplished by swapping it
+               with another pixel in its vicinity.
+            */
+            pixel const p2 = srcarray[ynew][xnew];
+                /* Original value of second pixel */
+
+            /* Set second pixel to new value */
+            PPM_ASSIGN(destarray[ynew][xnew],
+                       PPM_GETR(p), PPM_GETG(p), PPM_GETB(p));
+
+            /* Set first pixel to (old) value of second */
+            PPM_ASSIGN(destarray[row][col],
+                       PPM_GETR(p2), PPM_GETG(p2), PPM_GETB(p2));
+        } else {
+            /* Displaced pixel is out of bounds; leave the old pixel there.
+            */
+            PPM_ASSIGN(destarray[row][col],
+                       PPM_GETR(p), PPM_GETG(p), PPM_GETB(p));
+        }
+    }
+}
+
+
 
 int
-main(int    argc,
-     char * argv[]) {
+main(int          argc,
+     const char * argv[]) {
 
+    struct CmdlineInfo cmdline;
     FILE * ifP;
-    int argn, rows, cols;
+    int rows, cols;
     unsigned int row;
-    pixel ** destarray, ** srcarray;
-    pixel * pP;
-    pixel * pP2;
+    pixel ** destarray;
+    pixel ** srcarray;
     pixval maxval;
-    pixval r1, g1, b1;
-    int amount;
-    const char * const usage = "amount [ppmfile]\n        amount: # of pixels to displace a pixel by at most\n";
-
-    /* parse in 'default' parameters */
-    ppm_init(&argc, argv);
-
-    argn = 1;
-
-    /* parse in amount & seed */
-    if (argn == argc)
-        pm_usage(usage);
-    if (sscanf(argv[argn], "%d", &amount) != 1)
-        pm_usage(usage);
-    if (amount < 0)
-        pm_error("amount should be a positive number");
-    ++argn;
-
-    /* parse in filename (if present, stdin otherwise) */
-    if (argn != argc)
-    {
-        ifP = pm_openr(argv[argn]);
-        ++argn;
-    }
-    else
-        ifP = stdin;
+    struct pm_randSt randSt;
 
-    if (argn != argc)
-        pm_usage(usage);
+    pm_proginit(&argc, argv);
+
+    parseCommandLine(argc, argv, &cmdline);
+
+    ifP = pm_openr(cmdline.inputFilename);
 
     srcarray = ppm_readppm(ifP, &cols, &rows, &maxval);
 
     destarray = ppm_allocarray(cols, rows);
 
+    pm_randinit(&randSt);
+    pm_srand2(&randSt, cmdline.randomseedSpec, cmdline.randomseed);
+
     /* clear out the buffer */
     for (row = 0; row < rows; ++row)
         memset(destarray[row], 0, cols * sizeof(pixel));
 
-    srand(pm_randseed());
-
-    /* start displacing pixels */
+    /* Displace pixels */
     for (row = 0; row < rows; ++row) {
-        unsigned int col;
-        pP = &srcarray[row][0];
-
-        for (col = 0; col < cols; ++col) {
-            int const xdis = (rand() % (amount+1)) - ((amount+1) / 2);
-            int const ydis = (rand() % (amount+1)) - ((amount+1) / 2);
+        spreadRow(srcarray, cols, rows, cmdline.spread, row,
+                  destarray, &randSt);
 
-            int const xnew = col + xdis;
-            int const ynew = row + ydis;
-
-            /* only set the displaced pixel if it's within the bounds
-               of the image
-            */
-            if (xnew >= 0 && xnew < cols && ynew >= 0 && ynew < rows) {
-                /* displacing a pixel is accomplished by swapping it
-                   with another pixel in its vicinity - so, first
-                   store other pixel's RGB
-                */
-                pP2 = &srcarray[ynew][xnew];
-                r1 = PPM_GETR(*pP2);
-                g1 = PPM_GETG(*pP2);
-                b1 = PPM_GETB(*pP2);
-                /* set second pixel to new value */
-                pP2 = &destarray[ynew][xnew];
-                PPM_ASSIGN(*pP2, PPM_GETR(*pP), PPM_GETG(*pP), PPM_GETB(*pP));
-                
-                /* now, set first pixel to (old) value of second */
-                pP2 = &destarray[row][col];
-                PPM_ASSIGN(*pP2, r1, g1, b1);
-            } else {
-                /* displaced pixel is out of bounds; leave the old
-                   pixel there
-                */
-                pP2 = &destarray[row][col];
-                PPM_ASSIGN(*pP2, PPM_GETR(*pP), PPM_GETG(*pP), PPM_GETB(*pP));
-            }
-            ++pP;
-        }
     }
+    pm_randterm(&randSt);
 
     ppm_writeppm(stdout, destarray, cols, rows, maxval, 0);
 
@@ -115,3 +171,5 @@ main(int    argc,
 
     return 0;
 }
+
+