Index: /branches/eam_branches/ipp-20130904/psLib/src/imageops/psImageConvolve.c
===================================================================
--- /branches/eam_branches/ipp-20130904/psLib/src/imageops/psImageConvolve.c	(revision 36274)
+++ /branches/eam_branches/ipp-20130904/psLib/src/imageops/psImageConvolve.c	(revision 36275)
@@ -36,9 +36,6 @@
 #define MIN_GAUSS_FRAC 0.25             // Minimum Gaussian fraction to accept when smoothing
 
-
-
 static bool threaded = false;           // Run image convolution threaded?
 static pthread_mutex_t threadMutex = PTHREAD_MUTEX_INITIALIZER;
-
 
 static void kernelFree(psKernel *kernel)
@@ -2169,3 +2166,2 @@
     return threaded;
 }
-
Index: /branches/eam_branches/ipp-20130904/psLib/src/imageops/psImageConvolve.h
===================================================================
--- /branches/eam_branches/ipp-20130904/psLib/src/imageops/psImageConvolve.h	(revision 36274)
+++ /branches/eam_branches/ipp-20130904/psLib/src/imageops/psImageConvolve.h	(revision 36275)
@@ -35,4 +35,11 @@
     psVector *kernel;
 } psImageSmoothCacheData;
+
+/// a structure to contain data related to image smoothing with a pre-cached 1D gauss kernel
+typedef struct {
+    float Nsigma;
+    int Ns;				// number of pixel radii
+    float *radflux;			// conv kernel in special positions
+} psImageSmooth2dCacheData;
 
 /// A convolution kernel
Index: /branches/eam_branches/ipp-20130904/psLib/src/imageops/psImageConvolve2dCache.c
===================================================================
--- /branches/eam_branches/ipp-20130904/psLib/src/imageops/psImageConvolve2dCache.c	(revision 36275)
+++ /branches/eam_branches/ipp-20130904/psLib/src/imageops/psImageConvolve2dCache.c	(revision 36275)
@@ -0,0 +1,92 @@
+/// @file  psImageConvolve2dCache.c -- specialized 2D convolution for psf matching
+/// @author Eugene Magnier, IfA
+///
+/// @date $Date: 2009-02-05 23:56:14 $
+/// Copyright 2004-2007 Institute for Astronomy, University of Hawaii
+
+#ifdef HAVE_CONFIG_H
+#include "config.h"
+#endif
+
+#include <string.h>
+#include <math.h>
+#include "psAbort.h"
+#include "psMemory.h"
+#include "psLogMsg.h"
+#include "psError.h"
+#include "psAssert.h"
+#include "psScalar.h"
+#include "psBinaryOp.h"
+#include "psImageFFT.h"
+#include "psImageStructManip.h"
+#include "psImagePixelManip.h"
+#include "psTrace.h"
+#include "psThread.h"
+
+#include "psImageConvolve.h"
+
+// fixed set of radii and number of entries in circularly symmetric profile
+# define NRAD_MAX 44
+static float radii2[NRAD_MAX] = {  0.0,  1.0,  2.0,  4.0,  5.0, 
+				   8.0,  9.0, 10.0, 13.0, 16.0, 
+				   17.0, 18.0, 20.0, 25.0, 26.0, 
+				   29.0, 32.0, 34.0, 36.0, 37.0, 
+				   40.0, 41.0, 45.0, 49.0, 50.0, 
+				   52.0, 53.0, 58.0, 61.0, 64.0,
+				   65.0, 68.0, 72.0, 73.0, 74.0, 
+				   80.0, 81.0, 82.0, 85.0, 89.0,
+				   90.0, 97.0, 98.0, 100.0
+};
+
+static int   radiiN[NRAD_MAX] = {    1,  4, 4,  4,  8,
+				     4,  4, 8,  8,  4,
+				     8,  4, 8, 12,  8, 
+				     8,  4, 8,  4,  8,
+				     8,  8, 8,  4, 12,
+				     8,  8, 8,  8,  4, 
+				     16, 8, 4,  8,  8, 
+				      8, 4, 8, 16,  8, 
+				      8, 8, 4, 12
+};
+
+# define ADD_AXIS(RAD,DD) s += radflux[RAD]*(vi[p - DD] + vi[p + DD] + vi[p - DD*Nx] + vi[p + DD*Nx]);
+# define ADD_DIAG(RAD,DD) s += radflux[RAD]*(vi[p - DD - DD*Nx] + vi[p + DD - DD*Nx] + vi[p - DD + DD*Nx] + vi[p + DD + DD*Nx]);
+# define ADD_RAND(RAD,DX,DY) s += radflux[RAD]*(vi[p - DX - DY*Nx] + vi[p + DX - DY*Nx] + vi[p - DX + DY*Nx] + vi[p + DX + DY*Nx] + \
+						vi[p - DY - DX*Nx] + vi[p + DY - DX*Nx] + vi[p - DY + DX*Nx] + vi[p + DY + DX*Nx]);
+
+void psImageSmooth2dCacheDataFree (psImageSmooth2dCacheData *smdata) {
+    psFree (smdata->radflux);
+}
+
+// allocate the psImageSmooth2dCache data structure, but do not define the kernel
+psImageSmooth2dCacheData *psImageSmooth2dCacheAlloc (psImage *image, double sigma, double Nsigma) {
+
+    psImageSmooth2dCacheData *smdata = psAlloc(sizeof(psImageSmooth2dCacheData));
+    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;
+    }
+
+    // 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;
+}