Index: trunk/psLib/src/fits/psFitsImage.c
===================================================================
--- trunk/psLib/src/fits/psFitsImage.c	(revision 36760)
+++ trunk/psLib/src/fits/psFitsImage.c	(revision 37050)
@@ -345,14 +345,14 @@
         if (newScaleZero) {
             // Choose an appropriate BSCALE and BZERO
- 	  if (!psFitsScaleDetermine(bscale, bzero, boffset, bsoften, blank, image, mask, maskVal, fits)) {
-                // We can't have the write dying for this reason --- try to save it somehow!
-                psWarning("Unable to determine BSCALE and BZERO for image --- refusing to quantise.");
-                psErrorClear();
-                return psMemIncrRefCounter((psImage*)image);
-            }
+	    if (!psFitsScaleDetermine(bscale, bzero, boffset, bsoften, blank, image, mask, maskVal, fits)) {
+		// We can't have the write dying for this reason --- try to save it somehow!
+		psWarning("Unable to determine BSCALE and BZERO for image --- refusing to quantise.");
+		psErrorClear();
+		return psMemIncrRefCounter((psImage*)image);
+	    }
         } else {
-            // Don't want to muck with the current BSCALE and BZERO.  Get the current values and use those.
-            int status = 0;                 // Status of cfitsio
-            if (fits_read_key_dbl(fits->fd, "BSCALE", bscale, NULL, &status) && status != KEY_NO_EXIST) {
+	    // Don't want to muck with the current BSCALE and BZERO.  Get the current values and use those.
+	    int status = 0;                 // Status of cfitsio
+	    if (fits_read_key_dbl(fits->fd, "BSCALE", bscale, NULL, &status) && status != KEY_NO_EXIST) {
                 psFitsError(status, true, "Unable to read header.");
                 return NULL;
@@ -371,5 +371,25 @@
         }
 
-        return psFitsScaleForDisk(image, fits, *bscale, *bzero, *boffset, *bsoften, rng);
+	// XXX this test triggers the pswarp memory corruption problem.  you need to include "psThread.h" above to activate this
+	// for (int nx = 0; nx < 3; nx++) { 
+	//     psList *pending = psThreadGetPendingQueue();
+	//     if (pending && pending->head) {
+	// 	fprintf (stderr, "thread jobs are still pending...\n");
+	//     }
+	//     // XXX can we mess up the memory?
+	//     psImage *temp1 = psImageAlloc(image->numCols, image->numRows, PS_TYPE_F32);
+	//     fprintf (stderr, "corruption check set 1\n");
+	//     psMemCheckCorruption (stderr, true);
+	//     psImageInit (temp1, 2.0);
+	//     psFree (temp1);
+	//     fprintf (stderr, "corruption check set 2\n");
+	//     psMemCheckCorruption (stderr, true);
+	// }
+
+	// fprintf (stderr, "corruption check psFitsInsertImageWithMask\n");
+	// psMemCheckCorruption (stderr, true);
+
+        psImage *out = psFitsScaleForDisk(image, fits, *bscale, *bzero, *boffset, *bsoften, rng);
+        return out;
     }
 
@@ -638,6 +658,8 @@
     long blank = 0;                     // Blank (undefined) value for image
     psFitsFloat floatType;              // Custom floating-point convention type
-    psImage *diskImage = imageToDiskRepresentation(&bscale, &bzero, &boffset, &bsoften, &blank, &floatType, fits, image,
-                                                   mask, maskVal, NULL, true); // Image to write out
+
+    psImage *diskImage = imageToDiskRepresentation(&bscale, &bzero, &boffset, &bsoften, &blank, &floatType, 
+						   fits, image, mask, maskVal, NULL, true); // Image to write out
+
     if (!diskImage) {
         psError(PS_ERR_UNKNOWN, false, "Unable to convert image to desired disk format.");
@@ -666,4 +688,5 @@
         goto INSERT_DONE;
     }
+
     if (cfitsioBzero != 0.0) {
         psAssert(bzero == 0.0 && bscale == 1.0,
@@ -840,4 +863,5 @@
         }
     }
+
     if (psFitsError(status, true, "Could not write image to file.")) {
         success = false;
@@ -856,4 +880,5 @@
  INSERT_DONE:
     psFree(diskImage);
+
     if (compress) {
         // Restore compression state
@@ -906,4 +931,5 @@
     psImage *diskImage = imageToDiskRepresentation(&bscale, &bzero, &boffset, &bsoften, &blank, &floatType, fits, input,
                                                    mask, maskVal, NULL, false); // Image to write out
+    
     if (!diskImage) {
         psError(PS_ERR_UNKNOWN, false, "Unable to convert image to desired disk format.");
@@ -932,4 +958,5 @@
         goto UPDATE_DONE;
     }
+
     if (cfitsioBzero != 0.0) {
         psAssert(bzero == 0.0 && bscale == 1.0,
@@ -979,5 +1006,5 @@
         goto UPDATE_DONE;
     }
-
+    
     // We only want cfitsio to do the scale and zero if the type conversion requires it (e.g., input type is
     // an unsigned integer type).  In all other cases, we have already converted the image to use the
Index: trunk/psLib/src/fits/psFitsScale.c
===================================================================
--- trunk/psLib/src/fits/psFitsScale.c	(revision 36760)
+++ trunk/psLib/src/fits/psFitsScale.c	(revision 37050)
@@ -929,6 +929,5 @@
 
 
-psImage *psFitsScaleForDisk(const psImage *image, const psFits *fits, double bscale, double bzero, double boffset, double bsoften,
-                            psRandom *rng)
+psImage *psFitsScaleForDisk(const psImage *image, const psFits *fits, double bscale, double bzero, double boffset, double bsoften, psRandom *rng)
 {
     PS_ASSERT_IMAGE_NON_NULL(image, NULL);
@@ -971,81 +970,84 @@
     psImage *out = psImageAlloc(numCols, numRows, outType); // Output image
 
-    if (!psMemIncrRefCounter(rng) && options->fuzz) {
-        // Don't blab about which seed we're going to get --- it's not necessary for this purpose
-        rng = psRandomAlloc(PS_RANDOM_TAUS);
-    }
-
-#define SCALE_WRITE_OUT_CASE(IN, INTYPE, OUT, OUTTYPE) \
-    case PS_TYPE_##OUTTYPE: { \
-        double scale = 1.0 / bscale; \
-        double zero = bzero; \
+    // fprintf (stderr, "corruption check psFitsScaleForDisk 1\n");
+    // psMemCheckCorruption (stderr, true);
+
+    if (options->fuzz && !psMemIncrRefCounter (rng)) {
+      // Don't blab about which seed we're going to get --- it's not necessary for this purpose
+      rng = psRandomAlloc(PS_RANDOM_TAUS);
+    }
+
+#define SCALE_WRITE_OUT_CASE(IN, INTYPE, OUT, OUTTYPE)			\
+    case PS_TYPE_##OUTTYPE: {						\
+        double scale = 1.0 / bscale;					\
+        double zero = bzero;						\
         /* Note: BITPIX=8 treated differently, since it uses unsigned values; the rest use signed */ \
-        double min = bitpix == 8 ? 0 : -pow(2.0, options->bitpix - 1); \
+        double min = bitpix == 8 ? 0 : -pow(2.0, options->bitpix - 1);	\
         double max = bitpix == 8 ? 255 : (pow(2.0, options->bitpix - 1) - 1.0); \
-        for (int y = 0; y < numRows; y++) { \
-            for (int x = 0; x < numCols; x++) { \
-	      ps##INTYPE value; \
-	      if ((options->scaling == PS_FITS_SCALE_LOG_RANGE)||	\
-		(options->scaling == PS_FITS_SCALE_LOG_MANUAL)||        \
-		  (options->scaling == PS_FITS_SCALE_LOG_STDEV_POSITIVE)|| \
-		  (options->scaling == PS_FITS_SCALE_LOG_STDEV_NEGATIVE)|| \
-		  (options->scaling == PS_FITS_SCALE_LOG_STDEV_BOTH)) {	\
-		if (isfinite((IN)->data.INTYPE[y][x])) {		\
-		  value = log10( (IN)->data.INTYPE[y][x] - boffset );	\
+        for (int y = 0; y < numRows; y++) {				\
+            for (int x = 0; x < numCols; x++) {				\
+		ps##INTYPE value;					\
+		if ((options->scaling == PS_FITS_SCALE_LOG_RANGE)||	\
+		    (options->scaling == PS_FITS_SCALE_LOG_MANUAL)||	\
+		    (options->scaling == PS_FITS_SCALE_LOG_STDEV_POSITIVE)|| \
+		    (options->scaling == PS_FITS_SCALE_LOG_STDEV_NEGATIVE)|| \
+		    (options->scaling == PS_FITS_SCALE_LOG_STDEV_BOTH)) { \
+		    if (isfinite((IN)->data.INTYPE[y][x])) {		\
+			value = log10( (IN)->data.INTYPE[y][x] - boffset ); \
+		    }							\
+		    else {						\
+			value = (IN)->data.INTYPE[y][x];		\
+		    }							\
+		}							\
+		else if ((options->scaling == PS_FITS_SCALE_ASINH_RANGE)|| \
+			 (options->scaling == PS_FITS_SCALE_ASINH_MANUAL)|| \
+			 (options->scaling == PS_FITS_SCALE_ASINH_STDEV_POSITIVE)|| \
+			 (options->scaling == PS_FITS_SCALE_ASINH_STDEV_NEGATIVE)|| \
+			 (options->scaling == PS_FITS_SCALE_ASINH_STDEV_BOTH)) { \
+		    if (isfinite((IN)->data.INTYPE[y][x])) {		\
+			value = 1.0857362 * (asinh( ((IN)->data.INTYPE[y][x] - boffset) / (2.0 * bsoften))); \
+		    }							\
+		    else {						\
+			value = (IN)->data.INTYPE[y][x];		\
+		    }							\
 		}							\
 		else {							\
-		  value = (IN)->data.INTYPE[y][x];			\
+		    value = (IN)->data.INTYPE[y][x];			\
 		}							\
-	      }								\
-	      else if ((options->scaling == PS_FITS_SCALE_ASINH_RANGE)||	\
-		(options->scaling == PS_FITS_SCALE_ASINH_MANUAL)||        \
-		  (options->scaling == PS_FITS_SCALE_ASINH_STDEV_POSITIVE)|| \
-		  (options->scaling == PS_FITS_SCALE_ASINH_STDEV_NEGATIVE)|| \
-		  (options->scaling == PS_FITS_SCALE_ASINH_STDEV_BOTH)) {	\
-		if (isfinite((IN)->data.INTYPE[y][x])) {		\
-		  value = 1.0857362 * (asinh( ((IN)->data.INTYPE[y][x] - boffset) / (2.0 * bsoften))); \
+		if (!isfinite(value)) {					\
+		    /* This choice of "max" for non-finite pixels is mainly cosmetic --- it has to be */ \
+		    /* something, and "min" would produce holes in the cores of bright stars. */ \
+		    (OUT)->data.OUTTYPE[y][x] = max;			\
+		} else {						\
+		    value = (value - zero) * scale;			\
+		    if (options->fuzz && (value - (int)value != 0.0)) {	\
+			/* Add random factor [0.0,1.0): adds a variance of 1/12, */ \
+			/* but preserves the expectation value given the floor() */ \
+			value += psRandomUniform(rng) ;			\
+		    }							\
+		    /* Check for underflow and overflow; set either to max */ \
+		    (OUT)->data.OUTTYPE[y][x] = (value < min || value > max ? max : floor(value)); \
 		}							\
-		else {							\
-		  value = (IN)->data.INTYPE[y][x];			\
-		}							\
-	      }								\
-		else { \
-		  value = (IN)->data.INTYPE[y][x];	\
-		}					    \
-                if (!isfinite(value)) { \
-                    /* This choice of "max" for non-finite pixels is mainly cosmetic --- it has to be */ \
-                    /* something, and "min" would produce holes in the cores of bright stars. */ \
-                    (OUT)->data.OUTTYPE[y][x] = max; \
-                } else { \
-                    value = (value - zero) * scale; \
-                    if (options->fuzz && (value - (int)value != 0.0)) { \
-                       /* Add random factor [-0.5,0.5): adds a variance of 1/12, */ \
-                       /* but preserves the expectation value */ \
-                        value += psRandomUniform(rng) ; \
-                    } \
-                    /* Check for underflow and overflow; set either to max */ \
-                    (OUT)->data.OUTTYPE[y][x] = (value < min || value > max ? max : floor(value)); \
-                } \
-            } \
-        } \
-        break; \
-    }
-
-#define SCALE_WRITE_IN_CASE(IN, INTYPE, OUT) \
-    case PS_TYPE_##INTYPE: { \
-        switch (outType) { \
-	  SCALE_WRITE_OUT_CASE(IN, INTYPE, OUT, U8);;	\
-	  SCALE_WRITE_OUT_CASE(IN, INTYPE, OUT, S16);;	\
-	  SCALE_WRITE_OUT_CASE(IN, INTYPE, OUT, S32);;	\
-	  SCALE_WRITE_OUT_CASE(IN, INTYPE, OUT, S64);;	\
-          default: \
-            psAbort("Should be unreachable."); \
-        } \
-        break; \
-    }
-
+            }								\
+        }								\
+        break;								\
+    }
+    
+#define SCALE_WRITE_IN_CASE(IN, INTYPE, OUT)			\
+    case PS_TYPE_##INTYPE: {					\
+	switch (outType) {					\
+	    SCALE_WRITE_OUT_CASE(IN, INTYPE, OUT, U8);;		\
+	    SCALE_WRITE_OUT_CASE(IN, INTYPE, OUT, S16);;	\
+	    SCALE_WRITE_OUT_CASE(IN, INTYPE, OUT, S32);;	\
+	    SCALE_WRITE_OUT_CASE(IN, INTYPE, OUT, S64);;	\
+	  default:						\
+	    psAbort("Should be unreachable.");			\
+	}							\
+	break;							\
+    }
+    
     switch (image->type.type) {
-        SCALE_WRITE_IN_CASE(image, F32, out);
-        SCALE_WRITE_IN_CASE(image, F64, out); 
+	SCALE_WRITE_IN_CASE(image, F32, out);
+	SCALE_WRITE_IN_CASE(image, F64, out); 
       default:
         psAbort("Should be unreachable.");
@@ -1053,5 +1055,4 @@
 
     psFree(rng);
-
     return out;
 }
Index: trunk/psLib/src/imageops/psImageBackground.c
===================================================================
--- trunk/psLib/src/imageops/psImageBackground.c	(revision 36760)
+++ trunk/psLib/src/imageops/psImageBackground.c	(revision 37050)
@@ -124,5 +124,5 @@
         }
     } else if (nSubset >= nGoodPixels) {
-	    psTrace ("psLib.imageops", 4, "case 2: nSubset >= nGoodPixels (%d x %d : %d : %d : %d)\n", (int) nx, (int) ny, (int) nSubset, (int) nGoodPixels, (int) nPixels);
+	psTrace ("psLib.imageops", 4, "case 2: nSubset >= nGoodPixels (%d x %d : %d : %d : %d)\n", (int) nx, (int) ny, (int) nSubset, (int) nGoodPixels, (int) nPixels);
 	// in this case, we have to select from all masked pixels just to get the desired
 	// sample size
@@ -157,5 +157,9 @@
 	    double frnd = psRandomUniform(rng);
 	    int pixel = pixelVector->n * frnd;
-		
+	    // psAssert (pixel < pixelVector->nalloc, "oops, we went too far (1)");
+	    // psAssert (pixel >= 0, "oops, we went too far (2)");
+	    // psAssert (pixelVector->n - 1 < pixelVector->nalloc, "oops, we went too far (3)");
+	    // psAssert (pixelVector->n - 1 >= 0, "oops, we went too far (4)");
+
 	    psVectorAppend(values, pixelVector->data.F32[pixel]);
 	    pixelVector->data.F32[pixel] = pixelVector->data.F32[pixelVector->n - 1];
Index: trunk/psLib/src/sys/psMemory.c
===================================================================
--- trunk/psLib/src/sys/psMemory.c	(revision 36760)
+++ trunk/psLib/src/sys/psMemory.c	(revision 37050)
@@ -1315,6 +1315,11 @@
     checkingForCorruption = true;
 
+    // int nPrint = 0;
     psS32 nbad = 0;               // number of bad blocks
     for (psMemBlock *memBlock = (psMemBlock *) lastMemBlockAllocated; memBlock != NULL; memBlock = memBlock->nextBlock) {
+	// if (nPrint < 20) {
+	//     psMemBlockPrint (stderr, memBlock);
+	//     nPrint ++;
+	// }
         if (isBadMemBlock(output, memBlock, __FILE__, __LINE__, __func__)) {
             nbad++;
Index: trunk/psLib/src/sys/psThread.c
===================================================================
--- trunk/psLib/src/sys/psThread.c	(revision 36760)
+++ trunk/psLib/src/sys/psThread.c	(revision 37050)
@@ -152,4 +152,8 @@
     psThreadJob *job = psListGetAndRemove((psList *) done, PS_LIST_HEAD);
     return job;
+}
+
+psList *psThreadGetPendingQueue () {
+  return (psList *) pending;
 }
 
Index: trunk/psLib/src/sys/psThread.h
===================================================================
--- trunk/psLib/src/sys/psThread.h	(revision 36760)
+++ trunk/psLib/src/sys/psThread.h	(revision 37050)
@@ -90,4 +90,7 @@
 psThreadJob *psThreadJobGetDone(void);
 
+// utility function to access the pending queue
+psList *psThreadGetPendingQueue ();
+
 /// Allocate a thread task
 psThreadTask *psThreadTaskAlloc(const char *type, // Type of task