Index: /branches/rel10_ifa/psLib/src/imageops/psImageConvolve.c
===================================================================
--- /branches/rel10_ifa/psLib/src/imageops/psImageConvolve.c	(revision 6481)
+++ /branches/rel10_ifa/psLib/src/imageops/psImageConvolve.c	(revision 6482)
@@ -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.30.4.1 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2006-02-24 19:59:12 $
  *
  *  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
@@ -493,6 +493,6 @@
     PS_ASSERT_IMAGE_NON_NULL(image, NULL);
 
-    int Nx, Ny, Npixel, Nrange;
-    double factor, g, s;
+    int Nx, Ny, Npixel, Nrange, Nr;
+    double factor, g, s, sum;
     psVector *temp;
 
@@ -507,8 +507,13 @@
     #define IMAGESMOOTH_CASE(TYPE) \
 case PS_TYPE_##TYPE: { \
-        /* generate gaussian */ \
+        /* 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]; \
@@ -519,15 +524,37 @@
             ps##TYPE *vi = image->data.TYPE[j]; \
             ps##TYPE *vo = temp->data.TYPE; \
-            for (int i = 0; i < Nx; i++) { \
+            /* 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 = -Nrange; n < Nrange + 1; n++) { \
-                    if (i+n < 0) \
-                        continue; \
-                    if (i+n >= Nx) \
-                        continue; \
-                    s += gauss[n]*vi[i+n]; \
-                    g += gauss[n]; \
+                for (int n = -i; n < Nrange + 1; n++, vr++, vg++) { \
+                    s += *vg * *vr; \
+                    g += *vg; \
                 } \
-                vo[i] = s / g; \
+                *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)); \
@@ -536,26 +563,60 @@
         \
         /* smooth in Y direction */ \
-        temp = psVectorAlloc (image->numRows, PS_TYPE_##TYPE); \
-        for (int i = 0; i < Nx; i++) { \
-            ps##TYPE  *vo = temp->data.TYPE; \
-            ps##TYPE **vi = image->data.TYPE; \
-            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]; \
+        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; \
                 } \
-                vo[j] = s / g; \
-            } \
-            /* replace temp in image */ \
-            for (int j = 0; j < Ny; j++) { \
-                vi[j][i] = vo[j]; \
-            } \
-        } \
-        psFree (temp); \
+            } \
+            /* 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; \