Index: /branches/eam_branches/ipp-20130904/psLib/src/imageops/psImageConvolve2dCache.c
===================================================================
--- /branches/eam_branches/ipp-20130904/psLib/src/imageops/psImageConvolve2dCache.c	(revision 36276)
+++ /branches/eam_branches/ipp-20130904/psLib/src/imageops/psImageConvolve2dCache.c	(revision 36277)
@@ -65,28 +65,157 @@
     psMemSetDeallocator(smdata, (psFreeFunc) psImageSmooth2dCacheDataFree);
 
-    smdata->kernel = NULL;
-
-    if (!image) {
-	// relevant terms
-	smdata->Nrange = sigma*Nsigma + 0.5;    // Number of pixels either side for convolution kernel
-	smdata->Nx = 0;
-	smdata->Ny = 0;
-	smdata->resultX = NULL;
-	smdata->resultY = NULL;
-	return smdata;
-    }
+    smdata->radflux = NULL;
+
+    return smdata;
+}
+
+// generate a 2D smoothing kernel for supplied sigma & kappa (PS1_V1 profile).  sigma here
+// does not need to match that used to allocate the structure, but it is recommended
+bool psImageSmooth2dCacheKernel_PS1_V1 (psImageSmooth2dCacheData *smdata, float sigma, float kappa) {
+    // check for NULL structure elements?
+
+    int size = smdata->Nrange;
+
+    smdata->sigma = sigma;
+
+    smdata->Ns = (int)(smdata->Nsigma * sigma);
+    smdata->Ns = MAX (3, MIN (smdata->Ns, 10));
+
+    int Ns2 = Ns * Ns;
+
+    // we are going to use a hard-wired radial profile
+    smdata->radflux = psAlloc(sizeof(float)*NRAD_MAX);
+
+    float sum = 0.0;
+    for (i = 0; i < NRAD_MAX; i++) {
+      float z = radii2[i] / SQ(sigma);
+      smdata->radflux[i] = 1.0 / (1.0 + kappa*z + pow(z,1.666));
+      if (radii2[i] > Ns2) continue;
+      sum += radiiN[i] * smdata->radflux[i];
+    }
+    for (i = 0; i < NRAD_MAX; i++) {
+      smdata->radflux[i] = smdata->radflux[i] / sum;
+    }
+
+    return true;
+}
+
+// we can use the same DATA structure on multiple images of the same size
+bool psImageSmooth2dCache_F32(psImage *image, psImageSmooth2dCacheData *smdata)
+{
+    PS_ASSERT_IMAGE_NON_NULL(image, false);
+    PS_ASSERT_NON_NULL(smdata->radflux, false);
+    // assert on data type
 
     // relevant terms
-    smdata->Nrange = sigma*Nsigma + 0.5;    // Number of pixels either side for convolution kernel
-    smdata->Nx = image->numCols;            // Number of columns
-    smdata->Ny = image->numRows;            // Number of rows
-
-    // use a temp running buffer for X and Y directions.
-    smdata->resultX = psAlloc(smdata->Nx * sizeof(psF32));
-    memset (smdata->resultX, 0, smdata->Nx*sizeof(psF32));
-
-    smdata->resultY = psAlloc(smdata->Ny * sizeof(psF32));
-    memset (smdata->resultY, 0, smdata->Ny*sizeof(psF32));
-
-    return smdata;
-}
+    int Nrange = smdata->Nrange;    // Number of pixels either side for convolution kernel
+    int Nx = smdata->Nx;            // Number of columns
+    int Ny = smdata->Ny;            // Number of rows
+
+    psF32 *gauss = &smdata->kernel->data.F32[Nrange];
+    psF32 *resultX = smdata->resultX;
+    psF32 *resultY = smdata->resultY;
+       
+    /* Smooth in X direction */
+    {
+	for (int j = 0; j < Ny; j++) {
+	    psF32 *vi = image->data.F32[j];
+	    int xMax = PS_MIN(Nrange, Nx);
+	    /* Smooth first Nrange pixels, with renorm */
+	    for (int i = 0; i < xMax; i++, vi++) {
+		int convRange = PS_MIN(Nrange + 1, Nx - i);
+		psF32 *vr = vi - i;
+		psF32 *vg = gauss - i;
+		double g = 0.0;
+		double s = 0.0;
+		for (int n = -i; n < convRange; n++, vr++, vg++) {
+		    s += *vg * *vr;
+		    g += *vg;
+		}
+		resultX[i] = s / g;
+	    }
+	    /* If that's all the pixels we have, then we're done already */
+	    if (Nx > Nrange) {
+		/* Smooth middle pixels; if Nx < 2*Nrange, this pass is skipped */
+		for (int i = Nrange; i < Nx - Nrange; i++, vi++) {
+		    psF32 *vr = vi - Nrange;
+		    psF32 *vg = gauss - Nrange;
+		    double s = 0;
+		    for (int n = -Nrange; n < Nrange + 1; n++, vr++, vg++) {
+			s += *vg * *vr;
+		    }
+		    resultX[i] = s;
+		}
+		/* Smooth last Nrange pixels, with renorm */
+		// if Nx < 2*Nrange, this pass starts at i == Nrange
+		int xMin = PS_MAX(Nx - Nrange, Nrange);
+		for (int i = xMin; i < Nx; i++, vi++) {
+		    psF32 *vr = vi - Nrange;
+		    psF32 *vg = gauss - Nrange;
+		    double g = 0.0;
+		    double s = 0.0;
+		    for (int n = -Nrange; n < Nx - i; n++, vr++, vg++) {
+			s += *vg * *vr;
+			g += *vg;
+		    }
+		    resultX[i] = s / g;
+		}
+	    }
+	    memcpy(image->data.F32[j], resultX, Nx*sizeof(psF32));
+	}
+    }
+       
+    // this section probably hits the cache poorly for large images, but is probably OK for small ones
+    /* Smooth in Y direction */
+    {
+	for (int i = 0; i < Nx; i++) {
+	    int yMax = PS_MIN(Nrange, Ny);
+	    /* Smooth first Nrange pixels, with renorm */
+	    for (int j = 0; j < yMax; j++) {
+		int convRange = PS_MIN(Nrange + 1, Ny - j);
+		psF32 *vg = gauss - j;
+		double g = 0.0;
+		double s = 0.0;
+		for (int n = -j; n < convRange; n++, vg++) {
+		    psF32 vr = image->data.F32[j+n][i];
+		    s += *vg * vr;
+		    g += *vg;
+		}
+		resultY[j] = s / g;
+	    }
+	    /* If that's all the pixels we have, then we're done already */
+	    if (Ny > Nrange) {
+		/* Smooth middle pixels */
+		for (int j = Nrange; j < Ny - Nrange; j++) {
+		    psF32 *vg = gauss - Nrange;
+		    double s = 0;
+		    for (int n = -Nrange; n < Nrange + 1; n++, vg++) {
+			psF32 vr = image->data.F32[j+n][i]; 
+			s += *vg * vr;
+		    }
+		    resultY[j] = s;
+		}
+		/* Smooth last Nrange pixels, with renorm */
+		// if Ny < 2*Nrange, this pass starts at j == Nrange
+		int yMin = PS_MAX(Ny - Nrange, Nrange);
+		for (int j = yMin; j < Ny; j++) {
+		    psF32 *vg = gauss - Nrange;
+		    double g = 0.0;
+		    double s = 0.0;
+		    for (int n = -Nrange; n < Ny - j; n++, vg++) {
+			psF32 vr = image->data.F32[j+n][i];
+			s += *vg * vr;
+			g += *vg;
+		    }
+		    resultY[j] = s / g;
+		}
+	    }
+	    // loop here 
+	    for (int j = 0; j < Ny; j++) {
+		image->data.F32[j][i] = resultY[j];
+	    }
+	}
+    }
+    return true;
+}
+
Index: /branches/eam_branches/ipp-20130904/psModules/src/objects/pmPCMdata.c
===================================================================
--- /branches/eam_branches/ipp-20130904/psModules/src/objects/pmPCMdata.c	(revision 36276)
+++ /branches/eam_branches/ipp-20130904/psModules/src/objects/pmPCMdata.c	(revision 36277)
@@ -59,6 +59,6 @@
     psFree (pcm->constraint);
 
-    psFree (pcm->smdata1); // pre-allocated data for psImageSmooth_PreAlloc
-    psFree (pcm->smdata2); // pre-allocated data for psImageSmooth_PreAlloc
+    psFree (pcm->smdata); // pre-allocated data for psImageSmooth_PreAlloc
+    psFree (pcm->smdata2d); // pre-allocated data for psImageSmooth_PreAlloc
     return;
 }
@@ -90,6 +90,6 @@
     }
 
-    pcm->smdata1 = NULL;
-    pcm->smdata2 = NULL;
+    pcm->smdata = NULL;
+    pcm->smdata2d = NULL;
 
     pcm->modelConv = NULL;