Changeset 21280

Timestamp:

Feb 3, 2009, 4:55:27 PM (17 years ago)

Author:

Paul Price

Message:

Adding function to return an interpolation kernel. Intend to use this for calculating the covariance (pseudo-)matrix.

Location:

trunk/psLib/src/imageops

Files:

• psImageCovariance.c (modified) (3 diffs)
• psImageCovariance.h (modified) (3 diffs)
• psImageInterpolate.c (modified) (3 diffs)
• psImageInterpolate.h (modified) (3 diffs)

trunk/psLib/src/imageops/psImageCovariance.c

@@ -11 +11 @@
 
 #include "psImageCovariance.h"
+
+psKernel *psImageCovarianceNone(void)
+{
+    psKernel *covar = psKernelAlloc(0, 0, 0, 0); // Covariance pseudo-matrix
+    covar->kernel[0][0] = 1.0;
+    return covar;
+}
 
 
@@ -25 +32 @@
     // where M^x is the covariance matrix for x.
     // Note that the errors in f are correlated (covariance) even if the errors in x are not.
-    //
-    // We don't carry the entire covariance matrix for an image (the size goes as N^2, for N pixels, which
-    // makes storage difficult; and if that's not enough, the time to do the calculation is definitely
-    // impractical).  Since there are (generally) lots of zeros in the covariance matrix, and the same basic
-    // pattern repeats (for background pixels), we can just carry that pattern.  We carry this in a psKernel,
-    // since the values are the covariance between the pixel of consideration (at 0,0 in the kernel) and
-    // neighbouring pixels.  Note that this may not be strictly correct near sources, but this is the best we
-    // can do (and much better than most currently do).
 
     psKernel *covar;                    // Covariance matrix to use
@@ -108 +107 @@
 float psImageCovarianceFactor(const psKernel *covariance)
 {
-    return covariance ? covariance->kernel[0][0] : 1.0;
+    return covariance ? covariance->kernel[0][0] : NAN;
 }
 
+psKernel *psImageCovarianceAverage(const psArray *array)
+{
+    PS_ASSERT_ARRAY_NON_NULL(array, NULL);
+    PS_ASSERT_ARRAY_NON_EMPTY(array, NULL);
+
+    int xMin = INT_MAX, xMax = INT_MIN, yMin = INT_MAX, yMax = INT_MIN; // Range for covariance
+    int num = 0;                        // Number of good matrices to average
+    for (int i = 0; i < array->n; i++) {
+        psKernel *covar = array->data[i]; // Covariance matrix
+        if (!covar) {
+            continue;
+        }
+        xMin = PS_MIN(xMin, covar->xMin);
+        xMax = PS_MAX(xMax, covar->xMax);
+        yMin = PS_MIN(yMin, covar->yMin);
+        yMax = PS_MIN(yMax, covar->yMax);
+        num++;
+    }
+    if (num == 0) {
+        psError(PS_ERR_BAD_PARAMETER_SIZE, true, "No covariance matrices supplied for averaging.");
+        return NULL;
+    }
+
+    psKernel *average = psKernelAlloc(xMin, xMax, yMin, yMax); // Average covariance
+    for (int i = 0; i < array->n; i++) {
+        psKernel *covar = array->data[i]; // Covariance matrix
+        if (!covar) {
+            continue;
+        }
+        for (int y = yMin; y <= yMax; y++) {
+            for (int x = xMin; x <= xMax; x++) {
+                average->kernel[y][x] += covar->kernel[y][x];
+            }
+        }
+    }
+    psBinaryOp(average->image, average->image, "/", psScalarAlloc(num, PS_TYPE_F32));
+
+    return average;
+}

trunk/psLib/src/imageops/psImageCovariance.h

@@ -5 +5 @@
  * @author Paul Price, IfA
  *
- * @version $Revision: 1.1 $ $Name: not supported by cvs2svn $
- * @date $Date: 2009-01-28 22:16:33 $
+ * @version $Revision: 1.2 $ $Name: not supported by cvs2svn $
+ * @date $Date: 2009-02-04 02:55:27 $
  * Copyright 2009 Institute for Astronomy, University of Hawaii
  */
@@ -17 +17 @@
 
 #include <psImageConvolve.h>
+
+// We don't carry the entire covariance matrix for an image (the size goes as N^2, for N pixels, which makes
+// storage difficult; and if that's not enough, the time to do the calculation is definitely impractical).
+// Since there are (generally) lots of zeros in the covariance matrix, and the same basic pattern repeats (for
+// background pixels), we can just carry that pattern.  We carry this in a psKernel, since the values are the
+// covariance between the pixel of consideration (at 0,0 in the kernel) and neighbouring pixels.  Note that
+// this may not be strictly correct near sources, but this is the best we can do (and much better than most
+// currently do).
+
+/// Allocate a covariance pseudo-matrix with no covariance
+psKernel *psImageCovarianceNone(void);
 
 /// Calculate the covariance pseudo-matrix for a convolution kernel
@@ -29 +40 @@
     );
 
+/// Average many covariance pseudo-matrices
+psKernel *psImageCovarianceAverage(
+    const psArray *array                ///< Array of covariance pseudo-matrices
+    );
+
 
 /// @}

trunk/psLib/src/imageops/psImageInterpolate.c

@@ -7 +7 @@
  *  @author Paul Price, IfA
  *
- *  @version $Revision: 1.31 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2009-01-27 06:39:37 $
+ *  @version $Revision: 1.32 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2009-02-04 02:55:27 $
  *
  *  Copyright 2004-2007 Institute for Astronomy, University of Hawaii
@@ -866 +866 @@
 
 
-psImageInterpolateStatus psImageInterpolate(double *imageValue, double *varianceValue, psImageMaskType *maskValue,
-                                            float x, float y, const psImageInterpolation *interp)
+psImageInterpolateStatus psImageInterpolate(double *imageValue, double *varianceValue,
+                                            psImageMaskType *maskValue, float x, float y,
+                                            const psImageInterpolation *interp)
 {
     PS_ASSERT_PTR_NON_NULL(imageValue, PS_INTERPOLATE_STATUS_ERROR);
@@ -963 +964 @@
 
 
+psKernel *psImageInterpolationKernel(float x, float y, psImageInterpolateMode mode)
+{
+    int size = kernelSizes[mode];       // Size of kernel
+
+    // Kernel basics
+    INTERPOLATE_SETUP(x, y);
+    xExact = yExact = false;
+
+    psF32 xKernel[size], yKernel[size]; // Interpolation kernels
+    interpolationKernel(xKernel, xFrac, mode);
+    interpolationKernel(yKernel, yFrac, mode);
+
+    int min = -size/2, max = (size - 1) / 2; // Range for kernel
+    psKernel *kernel = psKernelAlloc(min, max, min, max); // Kernel to return
+
+    for (int y = 0; y < size; y++) {
+        for (int x = 0; x < size; x++) {
+            kernel->kernel[y][x] = yKernel[y] * xKernel[x];
+        }
+    }
+
+    return kernel;
+}
+
+
 psImageInterpolateMode psImageInterpolateModeFromString(const char *name)
 {

trunk/psLib/src/imageops/psImageInterpolate.h

@@ -7 +7 @@
  * @author Paul Price, Institute for Astronomy
  *
- * @version $Revision: 1.8 $ $Name: not supported by cvs2svn $
- * @date $Date: 2009-01-27 06:39:37 $
+ * @version $Revision: 1.9 $ $Name: not supported by cvs2svn $
+ * @date $Date: 2009-02-04 02:55:27 $
  * Copyright 2004-2007 Institute for Astronomy, University of Hawaii
  */
@@ -50 +50 @@
     const psImage *variance;            ///< Variance image for interpolation
     const psImage *mask;                ///< Mask image for interpolation
-    psImageMaskType maskVal;            ///< Value to mask
+    psImageMaskType maskVal;            ///< Value to mask
     double badImage;                    ///< Image value if x,y location is not good
     double badVariance;                 ///< Variance value if x,y location is not good
-    psImageMaskType badMask;            ///< Mask value to give bad pixels
-    psImageMaskType poorMask;           ///< Mask value to give poor pixels
+    psImageMaskType badMask;            ///< Mask value to give bad pixels
+    psImageMaskType poorMask;           ///< Mask value to give poor pixels
     float poorFrac;                     ///< Fraction of flux in bad pixels before output is marked bad
     bool shifting;                      ///< Shifting images? Don't interpolate if the shift is exact.
@@ -101 +101 @@
     );
 
+/// Generate the appropriate interpolation kernel
+psKernel *psImageInterpolationKernel(float x, float y, ///< Position of interest
+                                     psImageInterpolateMode mode ///< Interpolation mode
+    );
+
 #endif