Changeset 27732

Timestamp:

Apr 22, 2010, 6:20:17 PM (16 years ago)

Author:

Paul Price

Message:

Need to rescale the covariance matrix when we rescale the images. The previous covariance matrices were appropriate for the *chip* stage, not the warp stage (the pixel scales are currently different). Simulations show that this produces the expected detection limits (+/- 0.1 mag, which is about the expected error; more precision probably requires inserting and recovering fake sources).

Location:

branches/pap

Files:

• psLib/src/imageops/psImageCovariance.c (modified) (3 diffs)
• psLib/src/imageops/psImageCovariance.h (modified) (1 diff)
• pswarp/src/pswarpLoop.c (modified) (2 diffs)

branches/pap/psLib/src/imageops/psImageCovariance.c

@@ -11 +11 @@
 #include "psMemory.h"
 #include "psConstants.h"
+#include "psImageStructManip.h"
+#include "psImagePixelManip.h"
 #include "psImageConvolve.h"
 #include "psTrace.h"
@@ -16 +18 @@
 #include "psScalar.h"
 #include "psThread.h"
+#include "psImageInterpolate.h"
 
 #include "psImageCovariance.h"
 
 static bool threaded = false;           // Run threaded?
-
 
 psKernel *psImageCovarianceNone(void)
@@ -530 +532 @@
 
 
+psKernel *psImageCovarianceScale(const psKernel *in, float scale)
+{
+    // Trivial cases
+    if (!in) {
+        psKernel *out = psKernelAlloc(0, 0, 0, 0); // Output covariance
+        out->kernel[0][0] = 1.0;
+        return out;
+    }
+    PS_ASSERT_KERNEL_NON_NULL(in, NULL);
+    if (scale == 1.0) {
+        psImage *copy = psImageCopy(NULL, in->image, PS_TYPE_F32); // Copy of input covariance
+        psKernel *out = psKernelAllocFromImage(copy, -in->xMin, -in->yMin); // Output covariance
+        psFree(copy);
+        return out;
+    }
+
+    int xMinIn = in->xMin, xMaxIn = in->xMax, yMinIn = in->yMin, yMaxIn = in->yMax; // Input size
+    int xMinOut = (float)xMinIn / scale - 0.5, xMaxOut = (float)xMaxIn / scale + 0.5;     // Output size in x
+    int yMinOut = (float)yMinIn / scale - 0.5, yMaxOut = (float)yMaxIn / scale + 0.5;     // Output size in y
+
+    // Over-fill the covariance matrix so we're not troubled by edge effects
+    psKernel *overfill = psKernelAlloc(xMinIn - 1, xMaxIn + 1, yMinIn - 1, yMaxIn + 1); // Overfilled covar
+    psImageInit(overfill->image, 0.0);
+    int numOverlay = (xMaxIn - xMinIn + 1) * (yMaxIn - yMinIn + 1); // Number of pixels to overlay
+    if (psImageOverlaySection(overfill->image, in->image, 1, 1, "=") != numOverlay) {
+        psError(psErrorCodeLast(), false, "Unable to overfill covariance matrix.");
+        psFree(overfill);
+        return NULL;
+    }
+
+    psImageInterpolation *interp = psImageInterpolationAlloc(PS_INTERPOLATE_BILINEAR, overfill->image,
+                                                             NULL, NULL, 0, NAN, NAN, 0xFF, 0xFF,
+                                                             0.0, 0); // Interpolation
+    psFree(overfill);
+
+    // In transforming the positions, we get +0.5 to account for the centre of the pixels being at 0.5
+    // and +1 to account for the overfill.
+
+    psKernel *out = psKernelAlloc(xMinOut, xMaxOut, yMinOut, yMaxOut); // Output covariance
+    for (int y = yMinOut; y <= yMaxOut; y++) {
+        float yIn = y * scale + 0.5 - yMinIn + 1; // Position on input image (not the kernel)
+        for (int x = xMinOut; x <= xMaxOut; x++) {
+            float xIn = x * scale + 0.5 - xMinIn + 1; // Position on input (not the kernel)
+            double value;                                     // Value on output
+            if (!psImageInterpolate(&value, NULL, NULL, xIn, yIn, interp)) {
+                psError(psErrorCodeLast(), false, "Unable to interpolate kernel.");
+                return false;
+            }
+            out->kernel[y][x] = value;
+        }
+    }
+
+    psFree(interp);
+
+    return out;
+}
+
+
 bool psImageCovarianceSetThreads(bool set)
 {

branches/pap/psLib/src/imageops/psImageCovariance.h

@@ -90 +90 @@
     );
 
+
+/// Rescale a covariance matrix following a change in plate scale
+///
+/// The covariance matrix is stretched or shrunk to match the new plate scale.
+psKernel *psImageCovarianceScale(
+    const psKernel *in,                 ///< Input covariance pseudo-matrix
+    float scale                         ///< Scale factor (output plate scale relative to input plate scale)
+    );
+
 /// Control threading for image covariance functions
 ///

branches/pap/pswarp/src/pswarpLoop.c

@@ -267 +267 @@
     }
 
-    // Correct image for change in the plate scale
-    {
-        psAssert(input && input->fpa && input->fpa->toSky, "Require astrometry for input");
-        psAssert(outFPA && outFPA && outFPA->toSky, "Require astrometry for output");
-
-        double inScale = input->fpa->toSky->Xs + input->fpa->toSky->Ys; // Plate scale for input
-        double outScale = outFPA->toSky->Xs + outFPA->toSky->Ys; // Plate scale for output
-        float correction = PS_SQR(outScale / inScale); // Correction factor to apply to image
-        psLogMsg("pswarp", PS_LOG_INFO, "Correcting flux by %f to account for pixel scales", correction);
-
-        psBinaryOp(output->image, output->image, "*", psScalarAlloc(correction, PS_TYPE_F32));
-        if (output->variance) {
-            psBinaryOp(output->variance, output->variance, "*",
-                       psScalarAlloc(PS_SQR(correction), PS_TYPE_F32));
-        }
-    }
-
-
     // Set covariance matrix for output
     {
@@ -292 +274 @@
         psArray *covars = psListToArray(covariances); // Array of covariance matrices
         output->covariance = psImageCovarianceAverage(covars);
+        fprintf(stderr, "Covariance: %f\n", psImageCovarianceFactor(output->covariance));
         psFree(covars);
         psMetadataRemoveKey(output->analysis, PSWARP_ANALYSIS_COVARIANCES);
         psImageCovarianceTransfer(output->variance, output->covariance);
+    }
+
+    // Correct image for change in the plate scale
+    {
+        psAssert(input && input->fpa && input->fpa->toSky, "Require astrometry for input");
+        psAssert(outFPA && outFPA && outFPA->toSky, "Require astrometry for output");
+
+        double inScale = input->fpa->toSky->Xs + input->fpa->toSky->Ys; // Plate scale for input
+        double outScale = outFPA->toSky->Xs + outFPA->toSky->Ys; // Plate scale for output
+        float correction = PS_SQR(outScale / inScale); // Correction factor to apply to image
+        psLogMsg("pswarp", PS_LOG_INFO, "Correcting flux by %f to account for pixel scales", correction);
+        psBinaryOp(output->image, output->image, "*", psScalarAlloc(correction, PS_TYPE_F32));
+        if (output->variance) {
+            psBinaryOp(output->variance, output->variance, "*",
+                       psScalarAlloc(PS_SQR(correction), PS_TYPE_F32));
+        }
+        psKernel *covar = psImageCovarianceScale(output->covariance, outScale / inScale); // Scaled covariance
+        psFree(output->covariance);
+        output->covariance = covar;
     }