Index: trunk/psLib/src/imageops/psImageConvolve.c
===================================================================
--- trunk/psLib/src/imageops/psImageConvolve.c	(revision 4544)
+++ trunk/psLib/src/imageops/psImageConvolve.c	(revision 4815)
@@ -5,6 +5,6 @@
  *  @author Robert DeSonia, MHPCC
  *
- *  @version $Revision: 1.22 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2005-07-12 19:33:49 $
+ *  @version $Revision: 1.23 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2005-08-18 21:44:40 $
  *
  *  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
@@ -12,5 +12,5 @@
 
 #include <string.h>
-
+#include <math.h>
 #include "psImageConvolve.h"
 #include "psImageFFT.h"
@@ -279,5 +279,8 @@
 }
 
-psImage* psImageConvolve(psImage* out, const psImage* in, const psKernel* kernel, bool direct)
+psImage* psImageConvolve(psImage* out,
+                         const psImage* in,
+                         const psKernel* kernel,
+                         bool direct)
 {
     if (in == NULL) {
@@ -476,2 +479,74 @@
     return out;
 }
+
+void psImageSmooth (psImage *image,
+                    float sigma,
+                    float Nsigma)
+{
+
+    int Nx, Ny, Npixel, Nrange;
+    float factor, g, s;
+    psVector *temp;
+
+    // relevant terms
+    Nrange = sigma*Nsigma + 0.5;
+    Npixel = 2*Nrange + 1;
+    factor = -0.5/(sigma*sigma);
+
+    Nx = image->numCols;
+    Ny = image->numRows;
+
+    // generate gaussian
+    psVector *gaussnorm = psVectorAlloc (Npixel, PS_TYPE_F32);
+    for (int i = -Nrange; i < Nrange + 1; i++) {
+        gaussnorm->data.F32[i+Nrange] = exp (factor*i*i);
+    }
+    psF32 *gauss = &gaussnorm->data.F32[Nrange];
+
+    // smooth in X direction
+    temp = psVectorAlloc (Nx, PS_TYPE_F32);
+    for (int j = 0; j < Ny; j++) {
+        psF32 *vi = image->data.F32[j];
+        psF32 *vo = temp->data.F32;
+        for (int i = 0; i < Nx; i++) {
+            g = s = 0;
+            for (int n = -Nrange; n < Nrange + 1; n++) {
+                if (i+n < 0)
+                    continue;
+                if (i+n >= Nx)
+                    continue;
+                s += gauss[n]*vi[i+n];
+                g += gauss[n];
+            }
+            vo[i] = s / g;
+        }
+        memcpy (image->data.F32[j], temp->data.F32, Nx*sizeof(psF32));
+    }
+    psFree (temp);
+
+    // smooth in Y direction
+    temp = psVectorAlloc (image->numRows, PS_TYPE_F32);
+    for (int i = 0; i < Nx; i++) {
+        psF32  *vo = temp->data.F32;
+        psF32 **vi = image->data.F32;
+        for (int j = 0; j < Ny; j++) {
+            g = s = 0;
+            for (int n = -Nrange; n < Nrange + 1; n++) {
+                if (j+n < 0)
+                    continue;
+                if (j+n >= Ny)
+                    continue;
+                s += gauss[n]*vi[j+n][i];
+                g += gauss[n];
+            }
+            vo[j] = s / g;
+        }
+        // replace temp in image
+        for (int j = 0; j < Ny; j++) {
+            vi[j][i] = vo[j];
+        }
+    }
+    psFree (temp);
+    psFree (gaussnorm);
+}
+