Index: branches/pap/psLib/src/imageops/psImageConvolve.c
===================================================================
--- branches/pap/psLib/src/imageops/psImageConvolve.c	(revision 25238)
+++ branches/pap/psLib/src/imageops/psImageConvolve.c	(revision 25337)
@@ -678,4 +678,137 @@
 }
 
+static bool imageSmoothMaskPixels(psVector *out, const psImage *image, const psImage *mask,
+                                  psImageMaskType maskVal, const psVector *x, const psVector *y,
+                                  const psVector *gaussNorm, float minGauss, int size, int start, int stop)
+{
+    const psF32 *gauss = &gaussNorm->data.F32[size]; // Gaussian convolution kernel
+    int numCols = image->numCols, numRows = image->numRows; // Size of image
+    int xLast = numCols - 1, yLast = numRows - 1; // Last index
+    for (int i = start; i < stop; i++) {
+        int xPix = x->data.S32[i], yPix = y->data.S32[i]; // Pixel coordinates for smoothing
+
+        int yMin = PS_MAX(yPix - size, 0);
+        int yMax = PS_MIN(yPix + size, yLast);
+        int xMin = PS_MAX(xPix - size, 0);
+        int xMax = PS_MIN(xPix + size, xLast);
+
+        const float *yGauss = &gauss[yMin - yPix];
+
+        double ySumIG = 0.0, ySumG = 0.0;
+        for (int v = yMin; v <= yMax; v++, yGauss++) {
+            const float *xGauss = &gauss[xMin - xPix];
+            double xSumIG = 0.0, xSumG = 0.0;
+            const psImageMaskType *maskData = &mask->data.PS_TYPE_IMAGE_MASK_DATA[v][xMin];
+            const psF32 *imageData = &image->data.F32[v][xMin];
+            for (int u = xMin; u <= xMax; u++, xGauss++, imageData++, maskData++) {
+                if (*maskData & maskVal) {
+                    continue;
+                }
+                xSumIG += *imageData * *xGauss;
+                xSumG += *xGauss;
+            }
+            if (xSumG > minGauss) {
+                ySumIG += xSumIG * *yGauss;
+                ySumG += xSumG * *yGauss;
+            }
+        }
+
+        out->data.F32[i] = ySumG > minGauss ? ySumIG / ySumG : NAN;
+    }
+
+    return true;
+}
+
+static bool psImageSmoothMaskPixelsThread(psThreadJob *job)
+{
+    PS_ASSERT_THREAD_JOB_NON_NULL(job, false);
+    psAssert(job->args, "programming error: no job arguments");
+    psAssert(job->args->n == 11, "programming error: wrong number of job arguments");
+
+    psVector *out = job->args->data[0]; // Output vector
+    const psImage *image  = job->args->data[1]; // Input image
+    const psImage *mask   = job->args->data[2]; // Input mask
+    psImageMaskType maskVal = PS_SCALAR_VALUE(job->args->data[3], PS_TYPE_IMAGE_MASK_DATA);
+    const psVector *x = job->args->data[4];
+    const psVector *y = job->args->data[5];
+    const psVector *gaussNorm = job->args->data[6];
+    float minGauss = PS_SCALAR_VALUE(job->args->data[7], F32);
+    int size = PS_SCALAR_VALUE(job->args->data[8], S32);
+    int start = PS_SCALAR_VALUE(job->args->data[9], S32);
+    int stop = PS_SCALAR_VALUE(job->args->data[10], S32);
+    return imageSmoothMaskPixels(out, image, mask, maskVal, x, y, gaussNorm,
+                                 minGauss, size, start, stop);
+}
+
+
+psVector *psImageSmoothMaskPixels(const psImage *image, const psImage *mask, psImageMaskType maskVal,
+                                  const psVector *x, const psVector *y,
+                                  float sigma, float numSigma, float minGauss)
+{
+    PS_ASSERT_IMAGE_NON_NULL(image, NULL);
+    PS_ASSERT_IMAGE_NON_NULL(mask, NULL);
+    PS_ASSERT_IMAGE_TYPE(mask, PS_TYPE_IMAGE_MASK, NULL);
+    PS_ASSERT_IMAGES_SIZE_EQUAL(image, mask, NULL);
+    PS_ASSERT_VECTOR_NON_NULL(x, NULL);
+    PS_ASSERT_VECTOR_TYPE(x, PS_TYPE_S32, NULL);
+    PS_ASSERT_VECTOR_NON_NULL(y, NULL);
+    PS_ASSERT_VECTOR_TYPE(y, PS_TYPE_S32, NULL);
+    PS_ASSERT_VECTORS_SIZE_EQUAL(x, y, NULL);
+
+    int size = sigma * numSigma + 0.5;  // Half-size of kernel
+
+    int num = x->n;                     // Number of pixels to smooth
+    psVector *out = psVectorAlloc(num, PS_TYPE_F32); // Output results
+
+    // Generate normalized gaussian
+    IMAGE_SMOOTH_GAUSS(gaussNorm, size, sigma, F32);
+
+    // Columns
+    if (threaded) {
+        int numThreads = psThreadPoolSize(); // Number of threads
+        int delta = (numThreads) ? num / numThreads + 1 : num; // Block of cols to do at once
+        for (int start = 0; start < num; start += delta) {
+            int stop = PS_MIN(start + delta, num);  // End of block
+
+            psThreadJob *job = psThreadJobAlloc("PSLIB_IMAGE_SMOOTHMASK_PIXELS");
+            psArrayAdd(job->args, 1, out);
+            psArrayAdd(job->args, 1, (psImage*)image);
+            psArrayAdd(job->args, 1, (psImage*)mask);
+            PS_ARRAY_ADD_SCALAR(job->args, maskVal, PS_TYPE_IMAGE_MASK);
+            psArrayAdd(job->args, 1, (psVector*)x);
+            psArrayAdd(job->args, 1, (psVector*)y);
+            psArrayAdd(job->args, 1, gaussNorm);
+            PS_ARRAY_ADD_SCALAR(job->args, minGauss, PS_TYPE_F32);
+            PS_ARRAY_ADD_SCALAR(job->args, size, PS_TYPE_S32);
+            PS_ARRAY_ADD_SCALAR(job->args, start, PS_TYPE_S32);
+            PS_ARRAY_ADD_SCALAR(job->args, stop, PS_TYPE_S32);
+            if (!psThreadJobAddPending(job)) {
+                psFree(job);
+                psFree(gaussNorm);
+                psFree(out);
+                return NULL;
+            }
+            psFree(job);
+        }
+    } else if (!imageSmoothMaskPixels(out, image, mask, maskVal, x, y,
+                                      gaussNorm, minGauss, size, 0, num)) {
+        psError(PS_ERR_UNKNOWN, false, "Unable to smooth pixels.");
+        psFree(gaussNorm);
+        psFree(out);
+        return NULL;
+    }
+
+    if (threaded && !psThreadPoolWait(true)) {
+        psError(PS_ERR_UNKNOWN, false, "Error waiting for threads.");
+        psFree(gaussNorm);
+        psFree(out);
+        return NULL;
+    }
+
+    psFree(gaussNorm);
+    return out;
+}
+
+
 // Smooth an image with masked pixels
 // The calculation and calcMask images are *deliberately* backwards (row,col instead of col,row or
@@ -1325,7 +1458,7 @@
     if (threaded) {
         int numThreads = psThreadPoolSize(); // Number of threads
-        int deltaRows = (numThreads) ? numRows / numThreads : numRows; // Block of rows to do at once
+        int deltaRows = (numThreads) ? numRows / numThreads + 1 : numRows; // Block of rows to do at once
         for (int start = 0; start < numRows; start += deltaRows) {
-            int stop = PS_MIN (start + deltaRows, numRows);  // end of row block
+            int stop = PS_MIN(start + deltaRows, numRows);  // end of row block
 
             psThreadJob *job = psThreadJobAlloc("PSLIB_IMAGE_CONVOLVE_MASK");
@@ -1363,7 +1496,7 @@
     if (threaded) {
         int numThreads = psThreadPoolSize(); // Number of threads
-        int deltaCols = (numThreads) ? numCols / numThreads : numCols; // Block of cols to do at once
+        int deltaCols = (numThreads) ? numCols / numThreads + 1 : numCols; // Block of cols to do at once
         for (int start = 0; start < numCols; start += deltaCols) {
-            int stop = PS_MIN (start + deltaCols, numCols);  // end of col block
+            int stop = PS_MIN(start + deltaCols, numCols);  // end of col block
 
             psThreadJob *job = psThreadJobAlloc("PSLIB_IMAGE_CONVOLVE_MASK");
@@ -1486,4 +1619,10 @@
             psFree(task);
         }
+        {
+            psThreadTask *task = psThreadTaskAlloc("PSLIB_IMAGE_SMOOTHMASK_PIXELS", 9);
+            task->function = &psImageSmoothMaskPixelsThread;
+            psThreadTaskAdd(task);
+            psFree(task);
+        }
     } else if (!set && threaded) {
         psThreadTaskRemove("PSLIB_IMAGE_CONVOLVE_MASK");