Index: trunk/psLib/src/imageops/psImageConvolve.c
===================================================================
--- trunk/psLib/src/imageops/psImageConvolve.c	(revision 20460)
+++ trunk/psLib/src/imageops/psImageConvolve.c	(revision 20462)
@@ -7,6 +7,6 @@
 /// @author Eugene Magnier, IfA
 ///
-/// @version $Revision: 1.79 $ $Name: not supported by cvs2svn $
-/// @date $Date: 2008-10-29 19:46:31 $
+/// @version $Revision: 1.80 $ $Name: not supported by cvs2svn $
+/// @date $Date: 2008-10-29 20:43:18 $
 ///
 /// Copyright 2004-2007 Institute for Astronomy, University of Hawaii
@@ -711,6 +711,4 @@
 }
 
-
-
 psImage *psImageSmoothMask(psImage *output,
                            const psImage *image,
@@ -751,6 +749,329 @@
 
     switch (image->type.type) {
-        IMAGESMOOTH_MASK_CASE(F32);
-        IMAGESMOOTH_MASK_CASE(F64);
+      case PS_TYPE_F32: {
+	  psImage *calculation = psImageAlloc(numRows, numCols, PS_TYPE_F32); /* Calculation image; BW */
+	  psImage *calcMask = psImageAlloc(numRows, numCols, PS_TYPE_MASK); /* Mask for calculation image; BW */
+   
+	  /** Smooth in X direction **/
+	  for (int j = 0; j < numRows; j++) {
+	      for (int i = 0; i < numCols; i++) {
+		  int xMin = PS_MAX(i - size, 0);
+		  int xMax = PS_MIN(i + size, xLast);
+		  const psMaskType *maskData = &mask->data.PS_TYPE_MASK_DATA[j][xMin];
+		  const psF32 *imageData = &image->data.F32[j][xMin];
+		  int uMin = - PS_MIN(i, size); /* Minimum kernel index */
+		  const psF32 *gaussData = &gauss[uMin];
+		  double sumIG = 0.0, sumG = 0.0; /* Sums for convolution */
+		  for (int x = xMin; x <= xMax; x++, maskData++, imageData++, gaussData++) {
+		      if (*maskData & maskVal) {
+			  continue;
+		      }
+		      // if ((i == 1000) && (j > 10) && (j < 15)) {
+		      // 	  fprintf (stderr, "%d %d  %d  %f  %f   %f %f\n", i, j, x, *imageData, *gaussData, sumIG, sumG);
+		      // }
+		      sumIG += *imageData * *gaussData;
+		      sumG += *gaussData;
+		  }
+		  if (sumG > minGauss) {
+		      /* BW */
+		      calculation->data.F32[i][j] = sumIG / sumG;
+		      calcMask->data.PS_TYPE_MASK_DATA[i][j] = 0;
+		  } else {
+		      /* BW */
+		      calcMask->data.PS_TYPE_MASK_DATA[i][j] = 0xFF;
+		  }
+	      }
+	  }
+   
+	  output = psImageRecycle(output, numCols, numRows, PS_TYPE_F32);
+   
+	  /** Smooth in Y direction  **/
+	  for (int i = 0; i < numCols; i++) {
+	      for (int j = 0; j < numRows; j++) {
+		  int yMin = PS_MAX(j - size, 0);
+		  int yMax = PS_MIN(j + size, yLast);
+		  const psMaskType *maskData = &calcMask->data.PS_TYPE_MASK_DATA[i][yMin]; /* BW */
+		  const psF32 *imageData = &calculation->data.F32[i][yMin]; /* BW */
+		  int vMin = - PS_MIN(j, size); /* Minimum kernel index */
+		  const psF32 *gaussData = &gauss[vMin];
+		  double sumIG = 0.0, sumG = 0.0; /* Sums for convolution */
+		  for (int y = yMin; y <= yMax; y++, maskData++, imageData++, gaussData++) {
+		      if (*maskData) {
+			  continue;
+		      }
+		      sumIG += *imageData * *gaussData;
+		      sumG += *gaussData;
+		  }
+		  output->data.F32[j][i] = (sumG > minGauss) ? sumIG / sumG : NAN;
+	      }
+	  }
+   
+	  psFree(calculation);
+	  psFree(calcMask);
+	  psFree(gaussNorm);
+   
+	  return output;
+      }
+	IMAGESMOOTH_MASK_CASE(F64);
+      default:
+	psFree(gaussNorm);
+	char *typeStr;
+	PS_TYPE_NAME(typeStr,image->type.type);
+	psError(PS_ERR_BAD_PARAMETER_TYPE, true,
+		_("Specified psImage type, %s, is not supported."),
+		typeStr);
+	return false;
+    }
+
+    psAbort("Should never reach here.");
+    return false;
+}
+
+bool psImageSmoothMask_ScanRows_F32 (psImage *calculation, 
+				     psImage *calcMask, 
+				     const psImage *image, 
+				     const psImage *mask, 
+				     psMaskType maskVal,
+				     psVector *gaussNorm, 
+				     float minGauss, 
+				     int size, 
+				     int rowStart, 
+				     int rowStop) 
+{
+    const psF32 *gauss = &gaussNorm->data.F32[size]; // Gaussian convolution kernel
+
+    int numCols = image->numCols;
+    int xLast = numCols - 1; // Last index
+
+    for (int j = rowStart; j < rowStop; j++) {
+	for (int i = 0; i < numCols; i++) {
+	    int xMin = PS_MAX(i - size, 0);
+	    int xMax = PS_MIN(i + size, xLast);
+	    const psMaskType *maskData = &mask->data.PS_TYPE_MASK_DATA[j][xMin];
+	    const psF32 *imageData = &image->data.F32[j][xMin];
+	    int uMin = - PS_MIN(i, size); /* Minimum kernel index */
+	    const psF32 *gaussData = &gauss[uMin];
+	    double sumIG = 0.0, sumG = 0.0; /* Sums for convolution */
+	    for (int x = xMin; x <= xMax; x++, maskData++, imageData++, gaussData++) {
+		if (*maskData & maskVal) {
+		    continue;
+		}
+		sumIG += *imageData * *gaussData;
+		sumG += *gaussData;
+	    }
+	    if (sumG > minGauss) {
+		/* BW */
+		calculation->data.F32[i][j] = sumIG / sumG;
+		calcMask->data.PS_TYPE_MASK_DATA[i][j] = 0;
+	    } else {
+		/* BW */
+		calcMask->data.PS_TYPE_MASK_DATA[i][j] = 0xFF;
+	    }
+	}
+    }
+    return true;
+}
+
+bool psImageSmoothMask_ScanCols_F32 (psImage *output,
+				     psImage *calculation, 
+				     psImage *calcMask, 
+				     psMaskType maskVal,
+				     psVector *gaussNorm, 
+				     float minGauss, 
+				     int size, 
+				     int colStart, 
+				     int colStop) 
+{
+    const psF32 *gauss = &gaussNorm->data.F32[size]; // Gaussian convolution kernel
+
+    int numRows = output->numRows;
+    int yLast = numRows - 1; // Last index
+
+    for (int i = colStart; i < colStop; i++) {
+	for (int j = 0; j < numRows; j++) {
+	    int yMin = PS_MAX(j - size, 0);
+	    int yMax = PS_MIN(j + size, yLast);
+	    const psMaskType *maskData = &calcMask->data.PS_TYPE_MASK_DATA[i][yMin]; /* BW */
+	    const psF32 *imageData = &calculation->data.F32[i][yMin]; /* BW */
+	    int vMin = - PS_MIN(j, size); /* Minimum kernel index */
+	    const psF32 *gaussData = &gauss[vMin];
+	    double sumIG = 0.0, sumG = 0.0; /* Sums for convolution */
+	    for (int y = yMin; y <= yMax; y++, maskData++, imageData++, gaussData++) {
+		if (*maskData) {
+		    continue;
+		}
+		sumIG += *imageData * *gaussData;
+		sumG += *gaussData;
+	    }
+	    output->data.F32[j][i] = (sumG > minGauss) ? sumIG / sumG : NAN;
+	}
+    }
+    return true;
+}
+   
+bool psImageSmoothMask_ScanRows_F32_Threaded(psThreadJob *job)
+{
+    PS_ASSERT_THREAD_JOB_NON_NULL(job, false);
+    psAssert (job->args, "programming error: job args not alloced?");
+    psAssert (job->args->n == 10, "programming error: job Nargs mismatch");
+
+    psImage *calculation  = job->args->data[0]; // calculation image
+    psImage *calcMask     = job->args->data[1]; // calculation mask
+    const psImage *image  = job->args->data[2]; // input image
+    const psImage *mask   = job->args->data[3]; // input mask
+
+    psMaskType maskVal 	  = PS_SCALAR_VALUE(job->args->data[4],U8);
+    psVector *gaussNorm	  = job->args->data[5]; // gauss kernel
+    float minGauss    	  = PS_SCALAR_VALUE(job->args->data[6],F32);
+    int size           	  = PS_SCALAR_VALUE(job->args->data[7],S32);
+    int rowStart       	  = PS_SCALAR_VALUE(job->args->data[8],S32);
+    int rowStop        	  = PS_SCALAR_VALUE(job->args->data[9],S32);
+    return psImageSmoothMask_ScanRows_F32(calculation, calcMask, image, mask, maskVal, 
+					  gaussNorm, minGauss, size, 
+					  rowStart, rowStop);
+}
+
+bool psImageSmoothMask_ScanCols_F32_Threaded(psThreadJob *job)
+{
+    PS_ASSERT_THREAD_JOB_NON_NULL(job, false);
+    psAssert (job->args, "programming error: job args not alloced?");
+    psAssert (job->args->n == 9, "programming error: job Nargs mismatch");
+
+    psImage *output       = job->args->data[0]; // input image
+    psImage *calculation  = job->args->data[1]; // calculation image
+    psImage *calcMask     = job->args->data[2]; // calculation mask
+    psMaskType maskVal 	  = PS_SCALAR_VALUE(job->args->data[3],U8);
+
+    psVector *gaussNorm	  = job->args->data[4]; // gauss kernel
+    float minGauss    	  = PS_SCALAR_VALUE(job->args->data[5],F32);
+    int size           	  = PS_SCALAR_VALUE(job->args->data[6],S32);
+    int colStart       	  = PS_SCALAR_VALUE(job->args->data[7],S32);
+    int colStop        	  = PS_SCALAR_VALUE(job->args->data[8],S32);
+    return psImageSmoothMask_ScanCols_F32(output, calculation, calcMask, maskVal, 
+					  gaussNorm, minGauss, size,
+					  colStart, colStop);
+}
+
+psImage *psImageSmoothMask_Threaded(psImage *output,
+				    const psImage *image,
+				    const psImage *mask,
+				    psMaskType maskVal,
+				    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_MASK, NULL);
+    PS_ASSERT_IMAGES_SIZE_EQUAL(image, mask, NULL);
+
+    int numCols = image->numCols, numRows = image->numRows; // Size of image
+
+    int nThreads = psThreadPoolSize();
+    int scanCols = nThreads ? (numCols / 4 / nThreads) : numCols;
+    int scanRows = nThreads ? (numRows / 4 / nThreads) : numRows;
+
+    // relevant terms
+    int size = sigma * numSigma + 0.5;  // Half-size of kernel
+    int fullSize = 2*size + 1;          // Total number of pixels in convolution kernel
+
+    // Generate normalized gaussian
+    psVector *gaussNorm = psVectorAlloc(fullSize + 1, PS_TYPE_F32); // Normalised Gaussian
+    {
+        double sum = 0.0;
+        double norm = -0.5/(sigma*sigma);
+        for (int i = 0, u = -size; i <= fullSize; i++, u++) {
+            float value = expf(norm*PS_SQR(u));
+            gaussNorm->data.F32[i] = value;
+            sum += value;
+        }
+        sum = 1.0 / sum;
+        for (int i = 0; i <= fullSize; i++) {
+            gaussNorm->data.F32[i] *= sum;
+        }
+    }
+
+    switch (image->type.type) {
+	// Smooth an image with masked pixels
+	// The calculation and calcMask images are *deliberately* backwards (row,col instead of col,row or
+	// [col][row] instead of [row][col]) to optimise the iteration over rows; such instances are marked "BW"
+
+      case PS_TYPE_F32: {
+	  psImage *calculation = psImageAlloc(numRows, numCols, PS_TYPE_F32); /* Calculation image; BW */
+	  psImage *calcMask = psImageAlloc(numRows, numCols, PS_TYPE_MASK); /* Mask for calculation image; BW */
+   
+	  /** Smooth in X direction **/
+	  for (int rowStart = 0; rowStart < numRows; rowStart+=scanRows) {
+	      int rowStop = PS_MIN (rowStart + scanRows, numRows);
+
+	      // allocate a job, construct the arguments for this job
+	      psThreadJob *job = psThreadJobAlloc("PSLIB_IMAGE_SMOOTHMASK_SCANROWS");
+	      psArrayAdd(job->args, 1, calculation);
+	      psArrayAdd(job->args, 1, calcMask);
+	      psArrayAdd(job->args, 1, (psImage *) image); // cast away const
+	      psArrayAdd(job->args, 1, (psImage *) mask); // cast away const
+	      PS_ARRAY_ADD_SCALAR(job->args, maskVal,  PS_TYPE_U8);
+	      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, rowStart, PS_TYPE_S32);
+	      PS_ARRAY_ADD_SCALAR(job->args, rowStop,  PS_TYPE_S32);
+	      // -> psImageSmoothMask_ScanRows_F32 (calculation, calcMask, image, mask, maskVal, gauss, minGauss, size, rowStart, rowStop);
+
+	      // if threading is not active, we simply run the job and return
+	      if (!psThreadJobAddPending(job)) {
+		  psError(PS_ERR_UNKNOWN, false, "Unable to smooth image");
+		  psFree(job);
+		  return false;
+	      }
+	      psFree(job);
+
+	  }
+
+	  // wait here for the threaded jobs to finish (NOP if threading is not active)
+	  if (!psThreadPoolWait(true)) {
+	      psError(PS_ERR_UNKNOWN, false, "Unable to smooth image");
+	      return false;
+	  }
+
+	  output = psImageRecycle(output, numCols, numRows, PS_TYPE_F32);
+   
+	  /** Smooth in Y direction  **/
+	  for (int colStart = 0; colStart < numCols; colStart+=scanCols) {
+	      int colStop = PS_MIN (colStart + scanCols, numCols);
+
+	      // allocate a job, construct the arguments for this job
+	      psThreadJob *job = psThreadJobAlloc("PSLIB_IMAGE_SMOOTHMASK_SCANCOLS");
+	      psArrayAdd(job->args, 1, output);
+	      psArrayAdd(job->args, 1, calculation);
+	      psArrayAdd(job->args, 1, calcMask);
+	      PS_ARRAY_ADD_SCALAR(job->args, maskVal,  PS_TYPE_U8);
+	      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, colStart, PS_TYPE_S32);
+	      PS_ARRAY_ADD_SCALAR(job->args, colStop,  PS_TYPE_S32);
+	      // -> psImageSmoothMask_ScanCols_F32 (output, calculation, calcMask, maskVal, gauss, minGauss, size, colStart, colStop);
+
+	      // if threading is not active, we simply run the job and return
+	      if (!psThreadJobAddPending(job)) {
+		  psError(PS_ERR_UNKNOWN, false, "Unable to smooth image");
+		  psFree(job);
+		  return false;
+	      }
+	      psFree(job);
+	  }
+
+	  // wait here for the threaded jobs to finish (NOP if threading is not active)
+	  if (!psThreadPoolWait(true)) {
+	      psError(PS_ERR_UNKNOWN, false, "Unable to smooth image");
+	      return false;
+	  }
+	  psFree(calculation);
+	  psFree(calcMask);
+	  psFree(gaussNorm);
+   
+	  return output;
+      }
       default:
         psFree(gaussNorm);
@@ -1047,8 +1368,7 @@
     if (threaded) {
         int numThreads = psThreadPoolSize(); // Number of threads
-        float cols = (float)numCols / (float)numThreads; // Number of cols to do at once
-        for (int i = 0; i < numThreads; i++) {
-            int start = i * cols;       // Starting colunms
-            int stop = (i + 1) * cols;  // Stopping columns
+        int deltaRows = (numThreads) ? numRows / numThreads : 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
 
             psThreadJob *job = psThreadJobAlloc("PSLIB_IMAGE_CONVOLVE_MASK");
@@ -1060,5 +1380,5 @@
             PS_ARRAY_ADD_SCALAR(job->args, xMin, PS_TYPE_S32);
             PS_ARRAY_ADD_SCALAR(job->args, xMax, PS_TYPE_S32);
-            PS_ARRAY_ADD_SCALAR(job->args, 0x00, PS_TYPE_U8);
+            PS_ARRAY_ADD_SCALAR(job->args, 0x00, PS_TYPE_U8); // specify rows
             if (!psThreadJobAddPending(job)) {
                 psFree(job);
@@ -1086,12 +1406,11 @@
     if (threaded) {
         int numThreads = psThreadPoolSize(); // Number of threads
-        float cols = (float)numCols / (float)numThreads; // Number of columns to do at once
-        for (int i = 0; i < numThreads; i++) {
-            int start = i * cols;       // Starting column
-            int stop = (i + 1) * cols;  // Stopping column
+        int deltaCols = (numThreads) ? numCols / numThreads : 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
 
             psThreadJob *job = psThreadJobAlloc("PSLIB_IMAGE_CONVOLVE_MASK");
+            psArrayAdd(job->args, 1, out);
             psArrayAdd(job->args, 1, conv);
-            psArrayAdd(job->args, 1, out);
             PS_ARRAY_ADD_SCALAR(job->args, start, PS_TYPE_S32);
             PS_ARRAY_ADD_SCALAR(job->args, stop, PS_TYPE_S32);
@@ -1099,5 +1418,5 @@
             PS_ARRAY_ADD_SCALAR(job->args, yMin, PS_TYPE_S32);
             PS_ARRAY_ADD_SCALAR(job->args, yMax, PS_TYPE_S32);
-            PS_ARRAY_ADD_SCALAR(job->args, 0xff, PS_TYPE_U8);
+            PS_ARRAY_ADD_SCALAR(job->args, 0xff, PS_TYPE_U8); // specify cols
             if (!psThreadJobAddPending(job)) {
                 psFree(job);
@@ -1186,5 +1505,6 @@
 }
 
-
+// XXX for now, either thread all or none
+// have to call this before calling psImageSmoothMask_Threaded
 void psImageConvolveSetThreads(bool set)
 {
@@ -1196,10 +1516,22 @@
             psFree(task);
         }
+        {
+            psThreadTask *task = psThreadTaskAlloc("PSLIB_IMAGE_SMOOTHMASK_SCANROWS", 10);
+            task->function = &psImageSmoothMask_ScanRows_F32_Threaded;
+            psThreadTaskAdd(task);
+            psFree(task);
+        }
+        {
+            psThreadTask *task = psThreadTaskAlloc("PSLIB_IMAGE_SMOOTHMASK_SCANCOLS", 9);
+            task->function = &psImageSmoothMask_ScanCols_F32_Threaded;
+            psThreadTaskAdd(task);
+            psFree(task);
+        }
     } else if (!set && threaded) {
         psThreadTaskRemove("PSLIB_IMAGE_CONVOLVE_MASK");
-    }
-
+        psThreadTaskRemove("PSLIB_IMAGE_SMOOTHMASK_SCANROWS");
+        psThreadTaskRemove("PSLIB_IMAGE_SMOOTHMASK_SCANCOLS");
+    }
     threaded = set;
-
 }
 
Index: trunk/psLib/src/imageops/psImageConvolve.h
===================================================================
--- trunk/psLib/src/imageops/psImageConvolve.h	(revision 20460)
+++ trunk/psLib/src/imageops/psImageConvolve.h	(revision 20462)
@@ -5,6 +5,6 @@
  * @author Robert DeSonia, MHPCC
  *
- * @version $Revision: 1.36 $ $Name: not supported by cvs2svn $
- * @date $Date: 2008-10-17 01:32:49 $
+ * @version $Revision: 1.37 $ $Name: not supported by cvs2svn $
+ * @date $Date: 2008-10-29 20:43:18 $
  * Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
  */
@@ -208,4 +208,12 @@
 
 
+psImage *psImageSmoothMask_Threaded(psImage *output,
+				    const psImage *image,
+				    const psImage *mask,
+				    psMaskType maskVal,
+				    float sigma,
+				    float numSigma,
+				    float minGauss);
+
 bool psImageSmoothMaskF32(
     psImage *image,                    ///< the image to be smoothed