Index: trunk/psLib/src/imageops/psImageInterpolate.c
===================================================================
--- trunk/psLib/src/imageops/psImageInterpolate.c	(revision 18156)
+++ trunk/psLib/src/imageops/psImageInterpolate.c	(revision 18162)
@@ -7,6 +7,6 @@
  *  @author Paul Price, IfA
  *
- *  @version $Revision: 1.11 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2008-06-10 02:42:41 $
+ *  @version $Revision: 1.12 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2008-06-17 21:24:58 $
  *
  *  Copyright 2004-2007 Institute for Astronomy, University of Hawaii
@@ -132,26 +132,24 @@
 }
 
-// Interpolation engine using interpolation kernel
-static psImageInterpolateStatus interpolateKernel(double *imageValue, double *varianceValue,
-                                                  psMaskType *maskValue, float x, float y,
-                                                  const psImageInterpolateOptions *options)
-{
-    // Parameters have been checked by psImageInterpolate()
-
-    int xNum, yNum;                     // Number of interpolation kernel pixels
-    int xCentral, yCentral;             // Central pixel of the convolution
-    switch (options->mode) {
+// Setup for interpolation by kernel
+static inline void interpolateKernelSetup(int *xNum, int *yNum, // Size of interpolation kernel, returned
+                                          int *xCentral, int *yCentral, // Central pixel of convolution
+                                          float x, float y, // Coordinates of interest
+                                          psImageInterpolateMode mode // Mode for interpolation
+                                          )
+{
+    switch (mode) {
       case PS_INTERPOLATE_BILINEAR:
-        xNum = yNum = 2;
+        *xNum = *yNum = 2;
         // Central pixel is the pixel below the point of interest
-        xCentral = floor(x - 0.5 + FLT_EPSILON);
-        yCentral = floor(y - 0.5 + FLT_EPSILON);
+        *xCentral = floor(x - 0.5 + FLT_EPSILON);
+        *yCentral = floor(y - 0.5 + FLT_EPSILON);
         break;
       case PS_INTERPOLATE_BICUBE:
       case PS_INTERPOLATE_GAUSS:
-        xNum = yNum = 3;
+        *xNum = *yNum = 3;
         // Central pixel is the closest pixel to the point of interest
-        xCentral = x;
-        yCentral = y;
+        *xCentral = x;
+        *yCentral = y;
         break;
       case PS_INTERPOLATE_FLAT:
@@ -162,25 +160,15 @@
         psAbort("Invalid interpolation mode.");
     }
-
-    const psImage *image = options->image; // Image of interest
-    int xLast = image->numCols - 1;     // Last pixel in x
-    int yLast = image->numRows - 1;     // Last pixel in y
-
-    if (xCentral - (xNum - 1) / 2 < 0 || xCentral + xNum / 2 > xLast ||
-        yCentral - (yNum - 1) / 2 < 0 || yCentral + yNum / 2 > yLast) {
-        // At least one pixel of the interpolation kernel is off the image
-        if (imageValue) {
-            *imageValue = options->badImage;
-        }
-        if (varianceValue) {
-            *varianceValue = options->badVariance;
-        }
-        if (maskValue) {
-            *maskValue = options->badMask;
-        }
-        return PS_INTERPOLATE_STATUS_OFF;
-    }
-    double kernel[yNum][xNum];          // Interpolation kernel for straight interpolation
-    switch (options->mode) {
+}
+
+// Generate the interpolation kernel; it should be normalised to unity
+static inline void interpolateKernelGenerate(int xNum, int yNum, // Size of interpolation kernel
+                                             double kernel[xNum][yNum], // Kernel, to be set
+                                             int xCentral, int yCentral, // Central pixel of convolution
+                                             float x, float y, // Coordinates of interest
+                                             psImageInterpolateMode mode // Mode for interpolation
+                                             )
+{
+    switch (mode) {
       case PS_INTERPOLATE_BILINEAR: {
           double xFrac = x - 0.5 - xCentral; // Fraction of pixel in x
@@ -238,4 +226,38 @@
         psAbort("Invalid interpolation mode.");
     }
+}
+
+// Interpolation engine using interpolation kernel
+static psImageInterpolateStatus interpolateKernel(double *imageValue, double *varianceValue,
+                                                  psMaskType *maskValue, float x, float y,
+                                                  const psImageInterpolateOptions *options)
+{
+    // Parameters have been checked by psImageInterpolate()
+
+    int xNum, yNum;                     // Number of interpolation kernel pixels
+    int xCentral, yCentral;             // Central pixel of the convolution
+    interpolateKernelSetup(&xNum, &yNum, &xCentral, &yCentral, x, y, options->mode);
+
+    const psImage *image = options->image; // Image of interest
+    int xLast = image->numCols - 1;     // Last pixel in x
+    int yLast = image->numRows - 1;     // Last pixel in y
+
+    if (xCentral - (xNum - 1) / 2 < 0 || xCentral + xNum / 2 > xLast ||
+        yCentral - (yNum - 1) / 2 < 0 || yCentral + yNum / 2 > yLast) {
+        // At least one pixel of the interpolation kernel is off the image
+        if (imageValue) {
+            *imageValue = options->badImage;
+        }
+        if (varianceValue) {
+            *varianceValue = options->badVariance;
+        }
+        if (maskValue) {
+            *maskValue = options->badMask;
+        }
+        return PS_INTERPOLATE_STATUS_OFF;
+    }
+
+    double kernel[yNum][xNum];          // Interpolation kernel for straight interpolation
+    interpolateKernelGenerate(xNum, yNum, kernel, xCentral, yCentral, x, y, options->mode);
 
     // Image interpolation, according to image type
@@ -380,21 +402,23 @@
 }
 
-// Interpolation engine for separable interpolation kernels (either for good reasons or for practical reasons)
-static psImageInterpolateStatus interpolateSeparate(double *imageValue, double *varianceValue,
-                                                    psMaskType *maskValue, float x, float y,
-                                                    const psImageInterpolateOptions *options)
-{
-    // Parameters have been checked by psImageInterpolate()
-
-    int xNum, yNum;                     // Number of interpolation kernel pixels
-    switch (options->mode) {
+// Setup for interpolation by separable kernels
+static inline void interpolateSeparateSetup(int *xNum, int *yNum, // Size of interpolation kernel, returned
+                                            int *xCentral, int *yCentral, // Central pixel of convolution
+                                            float x, float y, // Coordinates of interest
+                                            psImageInterpolateMode mode // Mode for interpolation
+                                            )
+{
+    // Central pixel is the pixel below the point of interest
+    *xCentral = floor(x - 0.5);
+    *yCentral = floor(y - 0.5);
+    switch (mode) {
       case PS_INTERPOLATE_LANCZOS2:
-        xNum = yNum = 4;
+        *xNum = *yNum = 4;
         break;
       case PS_INTERPOLATE_LANCZOS3:
-        xNum = yNum = 6;
+        *xNum = *yNum = 6;
         break;
       case PS_INTERPOLATE_LANCZOS4:
-        xNum = yNum = 8;
+        *xNum = *yNum = 8;
         break;
       case PS_INTERPOLATE_FLAT:
@@ -405,54 +429,23 @@
         psAbort("Invalid interpolation mode.");
     }
-
-    // Central pixel is the pixel below the point of interest
-    int xCentral = floor(x - 0.5), yCentral = floor(y - 0.5); // Central pixel of the convolution
-    const psImage *image = options->image; // Image of interest
-    int xLast = image->numCols - 1;     // Last pixel in x
-    int yLast = image->numRows - 1;     // Last pixel in y
-
-    if (xCentral - (xNum - 1) / 2 < 0 || xCentral + xNum / 2 > xLast ||
-        yCentral - (yNum - 1) / 2 < 0 || yCentral + yNum / 2 > yLast) {
-        // At least one pixel of the interpolation kernel is off the image
-        if (imageValue) {
-            *imageValue = options->badImage;
-        }
-        if (varianceValue) {
-            *varianceValue = options->badVariance;
-        }
-        if (maskValue) {
-            *maskValue = options->badMask;
-        }
-        return PS_INTERPOLATE_STATUS_OFF;
-    }
-
-//    bool xExact, yExact;                // Is the shift exactly on?
-    double xKernel[xNum], yKernel[yNum];// Interpolation kernels in x and y
-    switch (options->mode) {
+}
+
+// Generate the interpolation kernels for separable case; they should be normalised to unity
+static inline void interpolateSeparateGenerate(int xNum, int yNum, // Size of interpolation kernel
+                                               double xKernel[xNum], double yKernel[yNum], // Kernels
+                                               int xCentral, int yCentral, // Central pixel of convolution
+                                               float x, float y, // Coordinates of interest
+                                               psImageInterpolateMode mode // Mode for interpolation
+                                               )
+{
+    // XXX Could put in an "exact shift" (i.e., xFrac = 0.0) version
+    switch (mode) {
       case PS_INTERPOLATE_LANCZOS2:
       case PS_INTERPOLATE_LANCZOS3:
       case PS_INTERPOLATE_LANCZOS4: {
           double xFrac = x - xCentral - 0.5; // Fraction of pixel in x
-#if 0
-          if (fabs(xFrac) < DBL_EPSILON) {
-              xExact = true;
-          } else {
-#endif
-              lanczos(xKernel, xNum, xFrac);
-#if 0
-              xExact = false;
-          }
-#endif
+          lanczos(xKernel, xNum, xFrac);
           double yFrac = y - yCentral - 0.5; // Fraction of pixel in y
-#if 0
-          if (fabs(yFrac) < DBL_EPSILON) {
-              yExact = true;
-          } else {
-#endif
-              lanczos(yKernel, yNum, yFrac);
-#if 0
-              yExact = false;
-          }
-#endif
+          lanczos(yKernel, yNum, yFrac);
           break;
       }
@@ -464,4 +457,38 @@
         psAbort("Invalid interpolation mode.");
     }
+}
+
+// Interpolation engine for separable interpolation kernels (either for good reasons or for practical reasons)
+static psImageInterpolateStatus interpolateSeparate(double *imageValue, double *varianceValue,
+                                                    psMaskType *maskValue, float x, float y,
+                                                    const psImageInterpolateOptions *options)
+{
+    // Parameters have been checked by psImageInterpolate()
+
+    int xNum, yNum;                     // Number of interpolation kernel pixels
+    int xCentral, yCentral; // Central pixel of the convolution
+    interpolateSeparateSetup(&xNum, &yNum, &xCentral, &yCentral, x, y, options->mode);
+
+    const psImage *image = options->image; // Image of interest
+    int xLast = image->numCols - 1;     // Last pixel in x
+    int yLast = image->numRows - 1;     // Last pixel in y
+
+    if (xCentral - (xNum - 1) / 2 < 0 || xCentral + xNum / 2 > xLast ||
+        yCentral - (yNum - 1) / 2 < 0 || yCentral + yNum / 2 > yLast) {
+        // At least one pixel of the interpolation kernel is off the image
+        if (imageValue) {
+            *imageValue = options->badImage;
+        }
+        if (varianceValue) {
+            *varianceValue = options->badVariance;
+        }
+        if (maskValue) {
+            *maskValue = options->badMask;
+        }
+        return PS_INTERPOLATE_STATUS_OFF;
+    }
+
+    double xKernel[xNum], yKernel[yNum]; // Interpolation kernels
+    interpolateSeparateGenerate(xNum, yNum, xKernel, yKernel, xCentral, yCentral, x, y, options->mode);
 
     // Image interpolation, according to image type
@@ -631,4 +658,98 @@
 }
 
+
+static float varianceFactorFlat(float x, float y, const psImageInterpolateOptions *options)
+{
+    // There's no smearing
+    return 1.0;
+}
+
+static float varianceFactorKernel(float x, float y, const psImageInterpolateOptions *options)
+{
+    int xNum, yNum;                     // Number of interpolation kernel pixels
+    int xCentral, yCentral;             // Central pixel of the convolution
+    interpolateKernelSetup(&xNum, &yNum, &xCentral, &yCentral, x, y, options->mode);
+
+    const psImage *image = options->image; // Image of interest
+    int xLast = image->numCols - 1;     // Last pixel in x
+    int yLast = image->numRows - 1;     // Last pixel in y
+
+    if (xCentral - (xNum - 1) / 2 < 0 || xCentral + xNum / 2 > xLast ||
+        yCentral - (yNum - 1) / 2 < 0 || yCentral + yNum / 2 > yLast) {
+        // At least one pixel of the interpolation kernel is off the image
+        return NAN;
+    }
+
+    double kernel[yNum][xNum];          // Interpolation kernel for straight interpolation
+    interpolateKernelGenerate(xNum, yNum, kernel, xCentral, yCentral, x, y, options->mode);
+
+    double sumKernel2 = 0.0;            // Sum of kernel squares
+    for (int j = 0; j < yNum; j++) {
+        for (int i = 0; i < xNum; i++) {
+            sumKernel2 += PS_SQR(kernel[j][i]);
+        }
+    }
+
+    return sumKernel2;
+}
+
+static float varianceFactorSeparate(float x, float y, const psImageInterpolateOptions *options)
+{
+    int xNum, yNum;                     // Number of interpolation kernel pixels
+    int xCentral, yCentral;             // Central pixel of the convolution
+    interpolateSeparateSetup(&xNum, &yNum, &xCentral, &yCentral, x, y, options->mode);
+
+    const psImage *image = options->image; // Image of interest
+    int xLast = image->numCols - 1;     // Last pixel in x
+    int yLast = image->numRows - 1;     // Last pixel in y
+
+    if (xCentral - (xNum - 1) / 2 < 0 || xCentral + xNum / 2 > xLast ||
+        yCentral - (yNum - 1) / 2 < 0 || yCentral + yNum / 2 > yLast) {
+        // At least one pixel of the interpolation kernel is off the image
+        return NAN;
+    }
+
+    double xKernel[xNum], yKernel[yNum]; // Interpolation kernels for separable interpolation
+    interpolateSeparateGenerate(xNum, yNum, xKernel, yKernel, xCentral, yCentral, x, y, options->mode);
+
+    double ySumKernel2 = 0.0;           // Sum of kernel squared in y
+    for (int j = 0; j < yNum; j++) {
+        double xSumKernel2 = 0.0;       // Sum of kernel squared in x
+        for (int i = 0; i < xNum; i++) {
+            xSumKernel2 += PS_SQR(xKernel[i]);
+        }
+        ySumKernel2 += xSumKernel2 * PS_SQR(yKernel[j]);
+    }
+
+    return ySumKernel2;
+}
+
+float psImageInterpolateVarianceFactor(float x, float y, const psImageInterpolateOptions *options)
+{
+    PS_ASSERT_PTR_NON_NULL(options, PS_INTERPOLATE_STATUS_ERROR);
+
+    switch (options->mode) {
+      case PS_INTERPOLATE_FLAT:
+        return varianceFactorFlat(x, y, options);
+      case PS_INTERPOLATE_BILINEAR:
+      case PS_INTERPOLATE_BICUBE:
+      case PS_INTERPOLATE_GAUSS:
+        return varianceFactorKernel(x, y, options);
+      case PS_INTERPOLATE_LANCZOS2:
+      case PS_INTERPOLATE_LANCZOS3:
+      case PS_INTERPOLATE_LANCZOS4:
+        return varianceFactorSeparate(x, y, options);
+      default:
+        psError(PS_ERR_BAD_PARAMETER_VALUE, true,
+                _("Specified interpolation method (%d) is not supported."),
+                options->mode);
+        return NAN;
+    }
+
+    psAbort("Should never reach here.");
+    return NAN;
+}
+
+
 psImageInterpolateMode psImageInterpolateModeFromString(const char *name)
 {
Index: trunk/psLib/src/imageops/psImageInterpolate.h
===================================================================
--- trunk/psLib/src/imageops/psImageInterpolate.h	(revision 18156)
+++ trunk/psLib/src/imageops/psImageInterpolate.h	(revision 18162)
@@ -7,6 +7,6 @@
  * @author Paul Price, Institute for Astronomy
  *
- * @version $Revision: 1.3 $ $Name: not supported by cvs2svn $
- * @date $Date: 2007-08-09 01:40:07 $
+ * @version $Revision: 1.4 $ $Name: not supported by cvs2svn $
+ * @date $Date: 2008-06-17 21:24:58 $
  * Copyright 2004-2007 Institute for Astronomy, University of Hawaii
  */
@@ -84,4 +84,12 @@
     );
 
+/// Return the variance factor for the appropriate position
+///
+/// psImageInterpolate sets the variance appropriate for extended regions (on the scale of the interpolation
+/// kernel), but this is not appropriate for pixel-to-pixel statistics.  This function returns the conversion
+/// factor.
+float psImageInterpolateVarianceFactor(float x, float y, ///< Position of interest
+                                       const psImageInterpolateOptions *options ///< Interpolation options
+    );
 
 #endif