Changeset 28006 for trunk

Timestamp:

May 18, 2010, 2:34:09 PM (16 years ago)

Author:

eugene

Message:

merging changes from PAP branch with modifications to the covariance matrix calculation

Location:

trunk

Files:

• . (modified) (1 prop)
• Ohana (modified) (1 prop)
• Ohana/src/relastro (modified) (1 prop)
• archive/noise_model (copied) (copied from branches/pap/archive/noise_model )
• archive/noise_model/simulate.c (copied) (copied from branches/pap/archive/noise_model/simulate.c )
• archive/noise_model/sub.subkernel (copied) (copied from branches/pap/archive/noise_model/sub.subkernel )
• archive/noise_model/test.11.kernel (copied) (copied from branches/pap/archive/noise_model/test.11.kernel )
• archive/noise_model/test.17.kernel (copied) (copied from branches/pap/archive/noise_model/test.17.kernel )
• archive/noise_model/test.23.kernel (copied) (copied from branches/pap/archive/noise_model/test.23.kernel )
• archive/noise_model/test.29.kernel (copied) (copied from branches/pap/archive/noise_model/test.29.kernel )
• archive/noise_model/test.35.kernel (copied) (copied from branches/pap/archive/noise_model/test.35.kernel )
• archive/noise_model/test.41.kernel (copied) (copied from branches/pap/archive/noise_model/test.41.kernel )
• archive/noise_model/test.47.kernel (copied) (copied from branches/pap/archive/noise_model/test.47.kernel )
• archive/noise_model/test.5.kernel (copied) (copied from branches/pap/archive/noise_model/test.5.kernel )
• ippconfig (modified) (1 prop)
• ppSim (modified) (1 prop)
• ppSub/src/ppSubReadoutSubtract.c (modified) (1 diff)
• psLib/src/imageops/psImageCovariance.c (modified) (3 diffs)
• psLib/src/imageops/psImageCovariance.h (modified) (1 diff)
• psModules (modified) (1 prop)
• psphot (modified) (1 prop)
• psphot/src/psphotEfficiency.c (modified) (2 diffs)
• pswarp/src/pswarp.h (modified) (2 diffs)
• pswarp/src/pswarpLoop.c (modified) (2 diffs)
• pswarp/src/pswarpTransformReadout.c (modified) (3 diffs)
• pswarp/src/pswarpTransformTile.c (modified) (7 diffs)

trunk/ppSub/src/ppSubReadoutSubtract.c

@@ -48 +48 @@
     pmReadout *outRO = pmFPAfileThisReadout(config->files, view, "PPSUB.OUTPUT");
     outRO->image = (psImage*)psBinaryOp(outRO->image, minuend->image, "-", subtrahend->image);
-    if (minuend->variance && subtrahend->variance) {
-        outRO->variance = (psImage*)psBinaryOp(outRO->variance, minuend->variance, "+", subtrahend->variance);
-    }
     outRO->mask = (psImage*)psBinaryOp(outRO->mask, minuend->mask, "|", subtrahend->mask);
+    outRO->variance = (psImage*)psBinaryOp(outRO->variance, minuend->variance, "+", subtrahend->variance);
 
+    // Measure the variance scales
+    psStats *varStats = psStatsAlloc(PS_STAT_ROBUST_MEDIAN); // Statistics for variance images
+    psRandom *rng = psRandomAlloc(PS_RANDOM_TAUS);           // Random number generator
+    psImageMaskType maskVal = pmConfigMaskGet("MASK.VALUE", config) |
+        pmConfigMaskGet("BLANK", config); // Bits to mask in inputs
+    psImageBackground(varStats, NULL, minuend->variance, minuend->mask, maskVal, rng);
+    float minuendVar = varStats->robustMedian; // Mean variance for minuend
+    psImageBackground(varStats, NULL, subtrahend->variance, subtrahend->mask, maskVal, rng);
+    float subtrahendVar = varStats->robustMedian; // Mean variance for subtrahend
+    psFree(varStats);
+    psFree(rng);
+
+    // Combine the covariances
+    // These are weighted by the appropriate mean variance.  This is probably not perfectly correct, but it
+    // does seem to reproduce the correct magnitude limit in psphot.
     psArray *covars = psArrayAlloc(2);  // Covariance pseudo-matrices
+    psVector *covarWeights = psVectorAlloc(2, PS_TYPE_F32); // Weights for covariances
     covars->data[0] = psMemIncrRefCounter(minuend->covariance);
     covars->data[1] = psMemIncrRefCounter(subtrahend->covariance);
-    outRO->covariance = psImageCovarianceAverage(covars);
+    covarWeights->data.F32[0] = minuendVar;
+    covarWeights->data.F32[1] = subtrahendVar;
+    outRO->covariance = psImageCovarianceAverageWeighted(covars, covarWeights);
     psFree(covars);
+    psFree(covarWeights);
+
+    psImageCovarianceTransfer(outRO->variance, outRO->covariance);
 
     outRO->data_exists = true;

trunk/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)
 {

trunk/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
 ///

trunk/psphot/src/psphotEfficiency.c

@@ -73 +73 @@
     psTrace("psphot.fake", 1, "Limiting flux: %f\n", fluxLim);
     psTrace("psphot.fake", 1, "Limiting mag: %f\n", *magLim);
+    psLogMsg("psphot", PS_LOG_INFO,
+             "Detection efficiency:\n"
+             "  Mean variance: %f * %f\n"
+             "  Threshold: %f\n"
+             "  Normalisation: %f\n"
+             "  Limiting magnitude: %f\n",
+             meanVar, *covarFactor, thresh, *norm, *magLim);
 
     *radius = smoothSigma * smoothNsigma;
@@ -166 +173 @@
     // loop over the available readouts
     for (int i = 0; i < num; i++) {
-        if (i == chisqNum) continue; // skip chisq image
-        if (!psphotEfficiencyReadout (config, view, "PSPHOT.INPUT", i, recipe)) {
+        if (i == chisqNum) continue; // skip chisq image
+        if (!psphotEfficiencyReadout (config, view, "PSPHOT.INPUT", i, recipe)) {
             psError (PSPHOT_ERR_CONFIG, false, "failed to measure detection efficiency for PSPHOT.INPUT entry %d", i);
-            return false;
-        }
+            return false;
+        }
     }
     return true;

trunk/pswarp/src/pswarp.h

@@ -29 +29 @@
 #define PSASTRO_RECIPE "PSASTRO" ///< Name of the recipe to use
 
-#define PSWARP_ANALYSIS_VARFACTOR "PSWARP.VARFACTOR" ///< Name for variance factor in analysis metadata
 #define PSWARP_ANALYSIS_GOODPIX   "PSWARP.GOODPIX" ///< Name for number of good pixels in analysis metadata
 #define PSWARP_ANALYSIS_COVARIANCES "PSWARP.COVARIANCES" ///< Name for covariance matrices on analysis MD
+#define PSWARP_ANALYSIS_JACOBIAN "PSWARP.JACOBIAN" ///< Name for Jacobian on analysis MD
 
 /**
@@ -70 +70 @@
     int xMin, xMax, yMin, yMax;         ///< Bounds of tile
     psKernel *covariance;               ///< Covariance matrix
+    double jacobian;                    ///< (Square root of) local Jacobian
 } pswarpTransformTileArgs;
 

trunk/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
     {
@@ -291 +273 @@
         psAssert(covariances, "Should be there");
         psArray *covars = psListToArray(covariances); // Array of covariance matrices
-        output->covariance = psImageCovarianceAverage(covars);
+        psKernel *covar = psImageCovarianceAverage(covars);
         psFree(covars);
         psMetadataRemoveKey(output->analysis, PSWARP_ANALYSIS_COVARIANCES);
+
+        // Correct covariance matrix scale for the mean (square root of the) Jacobian
+        double jacobian = psMetadataLookupF64(NULL, output->analysis, PSWARP_ANALYSIS_JACOBIAN); // Jacobian
+        int goodPixels = psMetadataLookupS32(NULL, output->analysis, PSWARP_ANALYSIS_GOODPIX);   // Good pixels
+        jacobian /= goodPixels;
+        output->covariance = psImageCovarianceScale(covar, jacobian);
+        psFree(covar);
+
         psImageCovarianceTransfer(output->variance, output->covariance);
     }

trunk/pswarp/src/pswarpTransformReadout.c

@@ -148 +148 @@
     psThreadJob *job = NULL;
     int xMin = output->image->numCols, xMax = 0, yMin = output->image->numRows, yMax = 0; // Bounds
-    int goodPixels = 0;                 // total number of good pixels across all tiles
+    int goodPixels = psMetadataLookupS32(&mdok, output->analysis, PSWARP_ANALYSIS_GOODPIX); // Number of pixels
     psList *covariances = psMetadataLookupPtr(&mdok, output->analysis,
                                               PSWARP_ANALYSIS_COVARIANCES); // Collection of covar. matrices
@@ -157 +157 @@
         psFree(covariances);            // Drop reference; still have the copy on the analysis metadata
     }
+    double jacobian = psMetadataLookupF64(&mdok, output->analysis, PSWARP_ANALYSIS_JACOBIAN); // Jacobian
+    if (!isfinite(jacobian)) {
+        jacobian = 0.0;
+    }
+
     while ((job = psThreadJobGetDone()) != NULL) {
         if (job->args->n < 1) {
@@ -171 +176 @@
                 psListAdd(covariances, PS_LIST_TAIL, args->covariance);
             }
+            if (args->goodPixels > 0 && isfinite(args->jacobian)) {
+                jacobian += args->jacobian * args->goodPixels;
+            }
         }
         psFree(job);
     }
     psFree(grid);
+
+    psMetadataAddS32(output->analysis, PS_LIST_TAIL, PSWARP_ANALYSIS_GOODPIX, PS_META_REPLACE,
+                     "Number of good pixels", goodPixels);
+    psMetadataAddF64(output->analysis, PS_LIST_TAIL, PSWARP_ANALYSIS_JACOBIAN, PS_META_REPLACE,
+                     "Jacobian of transformation", jacobian);
 
     if (xMin < xMax && yMin < yMax) {

trunk/pswarp/src/pswarpTransformTile.c

@@ -44 +44 @@
     args->yMax = PS_MIN_S32;
     args->covariance = NULL;
+    args->jacobian = NAN;
 
     return args;
@@ -50 +51 @@
 bool pswarpTransformTile(pswarpTransformTileArgs *args)
 {
-    psImage *inImage = args->input->image; ///< Input image
-    psImage *outImage = args->output->image; ///< Output image
+    pmReadout *input = args->input;        // Input readout
+    psImage *inImage = input->image, *inMask = input->mask; // Input images
+    pmReadout *output = args->output;      // Output readout
+    psImage *outImage = output->image, *outVariance = output->variance, *outMask = output->mask; // Outputs
 
     int inNumCols = inImage->numCols, inNumRows = inImage->numRows; ///< Size of input image
@@ -61 +64 @@
 
     // Dereference images for convenience
-    psF32 **outImageData     = args->output->image->data.F32;
-    psF32 **outVarData       = (args->output->variance) ? args->output->variance->data.F32 : NULL;
-    psImageMaskType **outMaskData = (args->output->mask)   ? args->output->mask->data.PS_TYPE_IMAGE_MASK_DATA : NULL;
-    psImageMaskType **inMaskData  = (args->input->mask)    ? args->input->mask->data.PS_TYPE_IMAGE_MASK_DATA : NULL;
+    psF32 **outImageData     = outImage->data.F32;
+    psF32 **outVarData       = outVariance ? outVariance->data.F32 : NULL;
+    psImageMaskType **outMaskData = outMask ? outMask->data.PS_TYPE_IMAGE_MASK_DATA : NULL;
+    psImageMaskType **inMaskData  = inMask ? inMask->data.PS_TYPE_IMAGE_MASK_DATA : NULL;
 
     pswarpMap *map = args->grid->maps[args->gridX][args->gridY]; ///< Map for this tile
@@ -74 +77 @@
     int yMin = PS_MAX(minPt.y, 0);
     int yMax = PS_MIN(maxPt.y, outNumRows);
+
+    double jacobian = map->Xx * map->Yy - map->Yx * map->Xy; // Jacobian of transformation
+    double jacobian2 = PS_SQR(jacobian);                     // Square Jacobian
 
     // Iterate over the output image pixels (parent frame)
@@ -87 +93 @@
             // pswarpMapApply converts the output coordinate (x,y) to the input coordinate.
             // both are in the parent frames of the input and output images.
-            double xInRaw, yInRaw;      // Input raw pixel coordinates
-            pswarpMapApply(&xInRaw, &yInRaw, map, x + 0.5, y + 0.5);
-            double xIn = xInRaw - outCol0, yIn = yInRaw - outRow0; // Position on input image
+            double xIn, yIn;            // Input pixel coordinates
+            pswarpMapApply(&xIn, &yIn, map, x + 0.5, y + 0.5);
             if (xIn < 0 || xIn >= inNumCols || yIn < 0 || yIn >= inNumRows) {
                 continue;
@@ -105 +110 @@
 
             if (outImageData) {
-                outImageData[yOut][xOut] = imageValue;
+                outImageData[yOut][xOut] = imageValue * jacobian;
             }
             if (outVarData) {
-                outVarData[yOut][xOut] = varValue;
+                outVarData[yOut][xOut] = varValue * jacobian2;
             }
             if (outMaskData) {
@@ -122 +127 @@
         double xIn, yIn;                // Position of interest on input
         pswarpMapApply(&xIn, &yIn, map, xOut + 0.5, yOut + 0.5);
-        // XXX Why are we subtracting the *output* col0,row0 from the *input* coordinates?
-        // I expect these are zero, so that it makes no difference, but it's distracting.
-        xIn -= outCol0;
-        yIn -= outRow0;
         psKernel *kernel = psImageInterpolationKernel(xIn, yIn, args->interp->mode); // Interpolation kernel
         args->covariance = psImageCovarianceCalculate(kernel, args->input->covariance);
+        args->jacobian = sqrt(jacobian);
         psFree(kernel);
     }