Index: trunk/psLib/src/imageops/psImageConvolve.c
===================================================================
--- trunk/psLib/src/imageops/psImageConvolve.c	(revision 10999)
+++ trunk/psLib/src/imageops/psImageConvolve.c	(revision 11153)
@@ -5,6 +5,6 @@
  *  @author Robert DeSonia, MHPCC
  *
- *  @version $Revision: 1.41 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2007-01-09 22:38:52 $
+ *  @version $Revision: 1.42 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2007-01-19 04:30:33 $
  *
  *  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
@@ -671,2 +671,116 @@
     return true;
 }
+
+bool psImageSmoothMaskF32 (psImage *image,
+                           psImage *mask,
+                           psMaskType maskVal,
+                           double  sigma,
+                           double  Nsigma)
+{
+    PS_ASSERT_IMAGE_NON_NULL(image, NULL);
+    PS_ASSERT_IMAGE_NON_NULL(mask, NULL);
+
+    // relevant terms
+    int Nrange = sigma*Nsigma + 0.5;    // Number of pixels either side for convolution kernel
+    int Npixel = 2*Nrange + 1;          // Total number of pixels in convolution kernel
+    int Nx = image->numCols;            // Number of columns
+    int Ny = image->numRows;            // Number of rows
+
+    /* generate normalized gaussian */
+    psVector *gaussnorm = psVectorAlloc(Npixel, PS_TYPE_F32);
+    {
+        double sum = 0.0;
+        double factor = -0.5/(sigma*sigma);
+        for (int i = -Nrange; i < Nrange + 1; i++) {
+            gaussnorm->data.F32[i+Nrange] = exp(factor*i*i);
+            sum += gaussnorm->data.F32[i+Nrange];
+        }
+        for (int i = -Nrange; i < Nrange + 1; i++) {
+            gaussnorm->data.F32[i+Nrange] /= sum;
+        }
+    }
+    psF32 *gauss = &gaussnorm->data.F32[Nrange];
+
+    /** Smooth in X direction **/
+    psVector *calculation = psVectorAlloc(Nx, PS_TYPE_F32);
+    for (int j = 0; j < Ny; j++) {
+        psU8  *vm = mask->data.U8[j];
+        psF32 *vi = image->data.F32[j];
+        psF32 *vo = calculation->data.F32;
+        // loop over all pixels in the row
+        for (int i = 0; i < Nx; i++, vi++, vo++, vm++) {
+            int offset = PS_MIN (i, Nrange);
+            psU8  *sm = vm - offset;
+            psF32 *si = vi - offset;
+            psF32 *sg = gauss - offset;
+            double g = 0.0;
+            double s = 0.0;
+            // loop over all valid pixels in the smoothing kernel
+            int xMin = PS_MAX (i - Nrange, 0);
+            int xMax = PS_MIN (i + Nrange + 1, Nx);
+            for (int n = xMin; n < xMax; n++, sm++, si++, sg++) {
+                if (*sm & maskVal)
+                    continue;
+                s += *sg * *si;
+                g += *sg;
+            }
+            *vo = (g > 0) ? s / g : 0.0;
+        }
+        memcpy(image->data.F32[j], calculation->data.F32, Nx*sizeof(psF32));
+    }
+    psFree(calculation);
+
+    /** Smooth in Y direction  **/
+    // allocate and save Nrange extra row vectors for storage. these will be
+    // written over the output rows only after we are Nrange rows beyond them
+    int Nsave = Nrange + 1;
+    psArray *rows = psArrayAlloc(Nsave);
+    for (int i = 0; i < Nsave; i++) {
+        rows->data[i] = psVectorAlloc(Nx, PS_TYPE_F32);
+    }
+
+    psVector *outsum = psVectorAlloc(Nx, PS_TYPE_F32);
+    for (int j = 0; j < Ny; j++) {
+        psVector *output = rows->data[j % Nsave];
+        memset (output->data.F32, 0, Nx*sizeof(psF32));
+        memset (outsum->data.F32, 0, Nx*sizeof(psF32));
+        int yMin = PS_MAX (j - Nrange, 0);
+        int yMax = PS_MIN (j + Nrange + 1, Ny);
+        for (int n = yMin; n < yMax; n++) {
+            psU8  *vm = mask->data.U8[n];
+            psF32 *vi = image->data.F32[n];
+            psF32 *vo = output->data.F32;
+            psF32 *vs = outsum->data.F32;
+            double g = gauss[n - j];
+            for (int i = 0; i < Nx; i++, vi++, vo++, vm++, vs++) {
+                if (*vm & maskVal)
+                    continue;
+                *vo += *vi * g;
+                *vs += g;
+            }
+        }
+        // renormalize the row
+        psF32 *vo = output->data.F32;
+        psF32 *vs = outsum->data.F32;
+        for (int i = 0; i < Nx; i++, vo++, vs++) {
+            *vo = (*vs > 0) ? *vo / *vs : 0.0;
+        }
+
+        // Write the output row
+        if (j - Nrange >= 0) {
+            int Nout = (j - Nrange) % Nsave;
+            psVector *save = rows->data[Nout];
+            memcpy(image->data.F32[j-Nrange], save->data.F32, Nx*sizeof(psF32));
+        }
+    }
+
+    // Write the remaining output rows
+    for (int j = PS_MAX(0, Ny - Nrange); j < Ny; j++) {
+        psVector *save = rows->data[j % Nsave];
+        memcpy(image->data.F32[j], save->data.F32, Nx*sizeof(psF32));
+    }
+    psFree(rows);
+    psFree(outsum);
+    psFree(gaussnorm);
+    return true;
+}