Index: trunk/psLib/src/imageops/psImageConvolve.c
===================================================================
--- trunk/psLib/src/imageops/psImageConvolve.c	(revision 6226)
+++ trunk/psLib/src/imageops/psImageConvolve.c	(revision 6937)
@@ -5,6 +5,6 @@
  *  @author Robert DeSonia, MHPCC
  *
- *  @version $Revision: 1.30 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2006-01-27 20:08:58 $
+ *  @version $Revision: 1.31 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2006-04-21 20:45:56 $
  *
  *  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
@@ -487,4 +487,6 @@
 }
 
+// XXX this is the old version from MHPCC
+# if 0
 bool psImageSmooth (psImage *image,
                     double  sigma,
@@ -576,3 +578,157 @@
     return true;
 }
-
+# else
+
+    // XXX EAM : we MUST use this version : it is a factor of >7 faster than the above
+    bool psImageSmooth (psImage *image,
+                        double  sigma,
+                        double  Nsigma)
+{
+    PS_ASSERT_IMAGE_NON_NULL(image, NULL);
+
+    int Nx, Ny, Npixel, Nrange, Nr;
+    double factor, g, s, sum;
+    psVector *temp;
+
+    // relevant terms
+    Nrange = sigma*Nsigma + 0.5;
+    Npixel = 2*Nrange + 1;
+    factor = -0.5/(sigma*sigma);
+
+    Nx = image->numCols;
+    Ny = image->numRows;
+
+    #define IMAGESMOOTH_CASE(TYPE) \
+case PS_TYPE_##TYPE: { \
+        /* generate normalized gaussian */ \
+        psVector *gaussnorm = psVectorAlloc (Npixel, PS_TYPE_##TYPE); \
+        sum = 0.0; \
+        for (int i = -Nrange; i < Nrange + 1; i++) { \
+            gaussnorm->data.TYPE[i+Nrange] = exp (factor*i*i); \
+            sum += gaussnorm->data.TYPE[i+Nrange]; \
+        } \
+        for (int i = -Nrange; i < Nrange + 1; i++) { \
+            gaussnorm->data.TYPE[i+Nrange] /= sum; \
+        } \
+        ps##TYPE *gauss = &gaussnorm->data.TYPE[Nrange]; \
+        \
+        /* smooth in X direction */ \
+        temp = psVectorAlloc (Nx, PS_TYPE_##TYPE); \
+        for (int j = 0; j < Ny; j++) { \
+            ps##TYPE *vi = image->data.TYPE[j]; \
+            ps##TYPE *vo = temp->data.TYPE; \
+            /* smooth first Nrange pixels, with renorm */ \
+            /* XXX need to check that this does not run over end for narrow images */ \
+            for (int i = 0; i < Nrange; i++, vi++, vo++) { \
+                ps##TYPE *vr = vi - i; \
+                ps##TYPE *vg = gauss - i; \
+                g = s = 0; \
+                for (int n = -i; n < Nrange + 1; n++, vr++, vg++) { \
+                    s += *vg * *vr; \
+                    g += *vg; \
+                } \
+                *vo = s / g; \
+            } \
+            /* smooth middle pixels */ \
+            for (int i = Nrange; i < Nx - Nrange; i++, vi++, vo++) { \
+                ps##TYPE *vr = vi - Nrange; \
+                ps##TYPE *vg = gauss - Nrange; \
+                s = 0; \
+                for (int n = -Nrange; n < Nrange + 1; n++, vr++, vg++) { \
+                    s += *vg * *vr; \
+                } \
+                *vo = s; \
+            } \
+            /* smooth last Nrange pixels, with renorm */ \
+            /* XXX does this miss the last column? */ \
+            for (int i = Nx - Nrange; i < Nx; i++, vi++, vo++) { \
+                ps##TYPE *vr = vi - Nrange; \
+                ps##TYPE *vg = gauss - Nrange; \
+                g = s = 0; \
+                for (int n = -Nrange; n < Nx - i - 1; n++, vr++, vg++) { \
+                    s += *vg * *vr; \
+                    g += *vg; \
+                } \
+                *vo = s / g; \
+            } \
+            memcpy (image->data.TYPE[j], temp->data.TYPE, Nx*sizeof(ps##TYPE)); \
+        } \
+        psFree (temp); \
+        \
+        /* smooth in Y direction */ \
+        psArray *temprows = psArrayAlloc (Nrange); \
+        temprows->n = Nrange; \
+        for (int j = 0; j < Nrange; j++) { \
+            temp = psVectorAlloc (Nx, PS_TYPE_##TYPE); \
+            /* zero the output row */ \
+            memset (temp->data.TYPE, 0, Nx*sizeof(PS_TYPE_##TYPE)); \
+            for (int n = -j, sum = 0; n < Nrange + 1; n++) { sum += gauss[n]; } \
+            for (int n = -j; n < Nrange + 1; n++) { \
+                ps##TYPE *vi = image->data.TYPE[j+n]; \
+                ps##TYPE *vo = temp->data.TYPE; \
+                g = gauss[n]; \
+                for (int i = 0; i < Nx; i++, vi++, vo++) { \
+                    *vo += *vi * g; \
+                } \
+            } \
+            /* save output rows on temp array */ \
+            temprows->data[j] = temp; \
+        } \
+        for (int j = Nrange; j < Ny - Nrange; j++) { \
+            /* save the Nrange-offset output row, then zero */ \
+            Nr = j % Nrange; \
+            temp = temprows->data[Nr]; \
+            memcpy (image->data.TYPE[j-Nrange], temp->data.TYPE, Nx*sizeof(ps##TYPE)); \
+            memset (temp->data.TYPE, 0, Nx*sizeof(PS_TYPE_##TYPE)); \
+            for (int n = -Nrange; n < Nrange + 1; n++) { \
+                ps##TYPE *vi = image->data.TYPE[j+n]; \
+                ps##TYPE *vo = temp->data.TYPE; \
+                g = gauss[n]; \
+                for (int i = 0; i < Nx; i++, vi++, vo++) { \
+                    *vo += *vi * g; \
+                } \
+            } \
+        } \
+        for (int j = Ny - Nrange; j < Ny; j++) { \
+            /* save the Nrange-offset output row, then zero */ \
+            Nr = j % Nrange; \
+            temp = temprows->data[Nr]; \
+            memcpy (image->data.TYPE[j-Nrange], temp->data.TYPE, Nx*sizeof(ps##TYPE)); \
+            memset (temp->data.TYPE, 0, Nx*sizeof(PS_TYPE_##TYPE)); \
+            for (int n = -Nrange, sum = 0; n < Ny - j - 1; n++) { sum += gauss[n]; } \
+            for (int n = -Nrange; n < Ny - j - 1; n++) { \
+                ps##TYPE *vi = image->data.TYPE[j+n]; \
+                ps##TYPE *vo = temp->data.TYPE; \
+                g = gauss[n]; \
+                for (int i = 0; i < Nx; i++, vi++, vo++) { \
+                    *vo += *vi * g; \
+                } \
+            } \
+        } \
+        for (int j = Ny; j < Ny + Nrange; j++) { \
+            /* save the Nrange-offset output row, then zero */ \
+            Nr = j % Nrange; \
+            temp = temprows->data[Nr]; \
+            memcpy (image->data.TYPE[j-Nrange], temp->data.TYPE, Nx*sizeof(ps##TYPE)); \
+        } \
+        psFree (temprows); \
+        psFree (gaussnorm); \
+        break; \
+    }
+
+    switch (image->type.type) {
+        IMAGESMOOTH_CASE(F32);
+        IMAGESMOOTH_CASE(F64);
+    default: {
+            char* typeStr;
+            PS_TYPE_NAME(typeStr,image->type.type);
+            psError(PS_ERR_BAD_PARAMETER_TYPE, true,
+                    PS_ERRORTEXT_psImage_IMAGE_TYPE_UNSUPPORTED,
+                    typeStr);
+            return false;
+        }
+    }
+    return true;
+}
+
+# endif