Changeset 28667

Timestamp:

Jul 13, 2010, 2:34:04 PM (16 years ago)

Author:

Paul Price

Message:

It's bothered me for a while that the noise model that we have been using is only kind of spatially variable. The variance map is intended to track the change in the noise properties as a function of position, and while it does keep track of different noise levels due to, e.g., different detectors, vignetting, etc., it has not tracked the different noise that comes from convolution kernels changing over an image. This had been an intentional oversight --- the prescription we were using for calculating the variance and covariance removed the spatial variation (by prescribing that the variance level under a convolution remain unchanged) and we were sweeping it under the carpet so we could calculate the covariance pseudo-matrix and use that for setting the absolute level of the noise.

I've now figured out how to have our cake and eat it too. By changing where the scalings are applied, the variance map can track the changing noise introduced by a varying convolution kernel, while the covariance pseudo-matrix tracks an absolute change. This fix is purely algorithmic --- it does not affect any of the variances/covariances stored on disk, nor is it affected by them. I hope that this means that any Magic recertification required is either minimal or completely unnecessary.

Again, all that has changed is where the scalings are applied.

Location:

trunk

Files:

• archive/noise_model/simulate.c (modified) (18 diffs)
• psLib/src/imageops/psImageCovariance.c (modified) (15 diffs)
• psLib/src/imageops/psImageInterpolate.c (modified) (4 diffs)
• psModules/src/imcombine/pmSubtraction.c (modified) (8 diffs)
• psphot/src/psphotSignificanceImage.c (modified) (1 diff)

trunk/archive/noise_model/simulate.c

@@ -6 +6 @@
 #define DET_SIZE 2000
 #define SKY_SIZE 1000
-#define SIZE 1024
 #define OUTPUT_ROOT "test"
-#define SCALE 0.654321
+#define SCALE 0.7654321
 #define ROT M_PI_2
-#define INTERPOLATION PS_INTERPOLATE_LANCZOS4
+#define INTERPOLATION PS_INTERPOLATE_LANCZOS3
 #define OFFSET 16
 #define SMOOTH_SIGMA 6.54321
@@ -16 +15 @@
 #define DUAL_KERNEL "sub.subkernel"
 #define WARP_NUM 10000
+#define CONV_NUM 1000
 
 static const float variances[] = { 3.0, 10.0, 30.0, 100.0, 300.0, 1000.0, 3000.0, 10000.0 };
@@ -21 +21 @@
 static const char *rootNames[] = { "o5298g0209o.fake.warp", "o5298g0210o.fake.warp", "o5298g0211o.fake.warp",
                                    "o5298g0212o.fake.warp", "o5298g0213o.fake.warp", "o5298g0214o.fake.warp",
+
                                    "o5298g0215o.fake.warp", "o5298g0216o.fake.warp" };
+#if 1
 static const char *kernels[] = { "test.5.kernel", "test.11.kernel", "test.17.kernel", "test.23.kernel",
                                  "test.29.kernel", "test.35.kernel", "test.41.kernel", "test.47.kernel" };
-
+#else
+static const char *kernels[] = { "test.11.kernel", "test.11.kernel", "test.11.kernel", "test.11.kernel",
+                                 "test.11.kernel", "test.11.kernel", "test.11.kernel", "test.11.kernel" };
+#endif
 
 void writeImage(const psImage *image, const char *suffix)
@@ -76 +81 @@
 float meanVar(const psImage *var, const psImage *mask, const psKernel *covar)
 {
-    psRandom *rng = psRandomAlloc(PS_RANDOM_TAUS);
     psStats *stats = psStatsAlloc(PS_STAT_SAMPLE_MEDIAN);
-    psImageBackground(stats, NULL, var, mask, 0xFF, rng);
+    psImageStats(stats, var, mask, 0xFF);
     float variance = stats->sampleMedian * psImageCovarianceFactor(covar);
     psFree(stats);
-    psFree(rng);
     return variance;
 }
@@ -93 +96 @@
         for (int x = 0; x < image->numCols; x++) {
             sn->data.F32[y][x] = image->data.F32[y][x] / sqrtf(var->data.F32[y][x] * varFactor);
+            if (!isfinite(sn->data.F32[y][x])) {
+                mask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] = 0xFF;
+            }
         }
     }
@@ -100 +106 @@
     float offset = stats->sampleMean;
     psFree(stats);
+
+    static int number = 0;
+
+    psString snName = NULL;
+    psStringAppend(&snName, "sn.%d.fits", number);
+    writeImage(sn, snName);
+    fprintf(stderr, "Writing S/N image %d\n", number);
 
     psHistogram *hist = psHistogramAlloc(-5, +5, 101);
@@ -121 +134 @@
     psFree(sn);
 
-    static int number = 0;
     psString name = NULL;
     fprintf(stderr, "Writing histogram %d\n", number);
-    psStringAppend(&name, "hist_%d.dat", number++);
+    psStringAppend(&name, "hist_%d.dat", number);
     FILE *file = fopen(name, "w");
     psFree(name);
@@ -136 +148 @@
     psFree(hist);
 
+    number++;
+
     return noise;
 }
@@ -144 +158 @@
 {
     psKernel *kernel2 = psKernelAlloc(kernel->xMin, kernel->xMax, kernel->yMin, kernel->yMax);
-    double sum = 0.0, sum2 = 0.0;
     for (int y = kernel->yMin; y <= kernel->yMax; y++) {
         for (int x = kernel->xMin; x <= kernel->xMax; x++) {
-            sum += kernel->kernel[y][x];
-            sum2 += kernel2->kernel[y][x] = PS_SQR(kernel->kernel[y][x]);
-        }
-    }
-    for (int y = kernel->yMin; y <= kernel->yMax; y++) {
-        for (int x = kernel->xMin; x <= kernel->xMax; x++) {
-            kernel2->kernel[y][x] /= sum2;
+            kernel2->kernel[y][x] = PS_SQR(kernel->kernel[y][x]);
         }
     }
@@ -176 +183 @@
 
     psKernel *kernel = psImageSmoothKernel(SMOOTH_SIGMA, SMOOTH_N_SIGMA); // Kernel used for smoothing
+    double sum2 = 0.0;                                               // Sum of kernel squared
+    for (int y = kernel->yMin; y <= kernel->yMax; y++) {
+        for (int x = kernel->xMin; x <= kernel->xMax; x++) {
+            sum2 += PS_SQR(kernel->kernel[y][x]);
+        }
+    }
+    float factor = 1.0 / sum2;
     psKernel *smoothCovar = psImageCovarianceCalculate(kernel, covar);
     psFree(kernel);
-    psImageCovarianceTransfer(smoothVariance, smoothCovar);
+
+    // Apply square root of significance image scaling factor to image
+    psBinaryOp(smoothImage, smoothImage, "*", psScalarAlloc(sqrtf(factor), PS_TYPE_F32));
+
 #else
     psKernel *kernel = psImageSmoothKernel(SMOOTH_SIGMA, SMOOTH_N_SIGMA); // Kernel used for smoothing
     psKernel *kernel2 = psKernelAlloc(kernel->xMin, kernel->xMax, kernel->yMin, kernel->yMax);
-    double sum = 0.0, sum2 = 0.0;
     for (int y = kernel->yMin; y <= kernel->yMax; y++) {
         for (int x = kernel->xMin; x <= kernel->xMax; x++) {
-            sum += kernel->kernel[y][x];
-            sum2 += kernel2->kernel[y][x] = PS_SQR(kernel->kernel[y][x]);
-        }
-    }
-    fprintf(stderr, "Kernel sum: %f %f\n", sum, sum2);
-    for (int y = kernel->yMin; y <= kernel->yMax; y++) {
-        for (int x = kernel->xMin; x <= kernel->xMax; x++) {
-            kernel2->kernel[y][x] /= sum2;
+            kernel2->kernel[y][x] = PS_SQR(kernel->kernel[y][x]);
         }
     }
@@ -332 +341 @@
 }
 
+float covarianceSum(const psKernel *covar)
+{
+    if (!covar) {
+        return 1.0;
+    }
+    double sum = 0.0;
+    for (int y = covar->yMin; y <= covar->yMax; y++) {
+        for (int x = covar->xMin; x <= covar->xMax; x++) {
+            sum += covar->kernel[y][x];
+        }
+    }
+    return sum;
+}
+
+
 int main(int argc, char *argv[])
 {
@@ -383 +407 @@
         }
 
-        fprintf(stderr, "Input image %d: S/N: %f Covar: %f Var: %f\n",
+        fprintf(stderr, "Input image %d: S/N: %f Covar: %f Var: %f CovarSum: %f\n",
                 i,
                 signoise(inImage, inMask, inVariance, inCovar),
                 psImageCovarianceFactor(inCovar),
-                meanVar(inVariance, inMask, inCovar));
+                meanVar(inVariance, inMask, inCovar),
+                covarianceSum(inCovar));
 
         phot(inImage, inMask, inVariance, inCovar);
@@ -417 +442 @@
         }
 
-        fprintf(stderr, "Warp image %d: S/N: %f Covar: %f Var: %f\n",
+        fprintf(stderr, "Warp image %d: S/N: %f Covar: %f Var: %f CovarSum: %f\n",
                 i,
                 signoise(warpImage, warpMask, warpVariance, warpCovar),
                 psImageCovarianceFactor(warpCovar),
-                meanVar(warpVariance, warpMask, warpCovar));
+                meanVar(warpVariance, warpMask, warpCovar),
+                covarianceSum(warpCovar));
 
         phot(warpImage, warpMask, warpVariance, warpCovar);
@@ -434 +460 @@
         pmReadoutReadSubtractionKernels(ro, fits);
         pmSubtractionKernels *kernels = psMetadataLookupPtr(NULL, ro->analysis, PM_SUBTRACTION_ANALYSIS_KERNEL);
+        int normIndex = PM_SUBTRACTION_INDEX_NORM(kernels);
+        int bgIndex = PM_SUBTRACTION_INDEX_BG(kernels);
+#if 1
+//        kernels->solution1->data.F64[normIndex] += 1.0;
+        kernels->solution1->data.F64[bgIndex] = 0.0;
+#endif
+        fprintf(stderr, "Norm: %f BG: %f\n", kernels->solution1->data.F64[normIndex], kernels->solution1->data.F64[bgIndex]);
 #if 0
         for (int i = 0; i < kernels->num; i++) {
@@ -455 +488 @@
         psFitsClose(fits);
 
+#if 0
+        {
+            psRandom *rng = psRandomAlloc(PS_RANDOM_TAUS);
+            psArray *covariances = psArrayAlloc(CONV_NUM);
+            psVector *factors = psVectorAlloc(CONV_NUM, PS_TYPE_F32);
+            double mean = 0.0;
+            for (int i = 0; i < CONV_NUM; i++) {
+                float x, y;
+                p_pmSubtractionPolynomialNormCoords(&x, &y, psRandomUniform(rng) * SKY_SIZE,
+                                                    psRandomUniform(rng) * SKY_SIZE,
+                                                    kernels->xMin, kernels->xMax,
+                                                    kernels->yMin, kernels->yMax);
+                psKernel *kernel = pmSubtractionKernel(kernels, x, y, false);
+                psKernel *covar = covariances->data[i] = psImageCovarianceCalculate(kernel, inCovar);
+                psFree(kernel);
+                mean += factors->data.F32[i] = psImageCovarianceFactor(covar);
+            }
+            psFree(rng);
+            psFree(covariances);
+
+            mean /= WARP_NUM;
+
+            double stdev = 0.0;
+            for (int i = 0; i < CONV_NUM; i++) {
+                stdev += PS_SQR(factors->data.F32[i] - mean);
+            }
+            stdev = sqrt(stdev/(CONV_NUM-1));
+            fprintf(stderr, "Conv covariance mean: %f stdev: %f\n", mean, stdev);
+            psFree(factors);
+        }
+#endif
+
         pmReadout *conv = pmReadoutAlloc(NULL);
 #if 1
@@ -460 +525 @@
         conv->mask = psImageAlloc(SKY_SIZE, SKY_SIZE, PS_TYPE_IMAGE_MASK);
         conv->variance = psImageAlloc(SKY_SIZE, SKY_SIZE, PS_TYPE_F32);
-        if (!pmSubtractionMatchPrecalc(NULL, conv, NULL, ro, ro->analysis, 32, 0.0, 0.01,
+        if (!pmSubtractionMatchPrecalc(NULL, conv, NULL, ro, ro->analysis, 2 * kernels->size + 1, 0.0, 0.01,
                                        0xFF, 0xF0, 0x0F, 0.1, 1.0)) {
             psErrorStackPrint(stderr, "Error:");
@@ -494 +559 @@
         }
 
-        fprintf(stderr, "Conv Image %d: S/N: %f Covar: %f Var: %f\n",
+        fprintf(stderr, "Conv Image %d: S/N: %f Covar: %f Var: %f CovarSum: %f\n",
                 i,
                 signoise(conv->image, conv->mask, conv->variance, conv->covariance),
                 psImageCovarianceFactor(conv->covariance),
-                meanVar(conv->variance, conv->mask, conv->covariance));
+                meanVar(conv->variance, conv->mask, conv->covariance),
+                covarianceSum(conv->covariance));
 
         phot(conv->image, conv->mask, conv->variance, conv->covariance);
         readouts->data[i] = conv;
 
+        //        exit(1);
     }
 
@@ -624 +691 @@
                 meanVar(diffVariance, diffMask, diffCovar));
 
-        phot(diffImage, diffMask, diffVariance, diffCovar);
-
         writeImage(diffImage, "ssdiff.image.fits");
         writeImage(diffMask, "ssdiff.mask.fits");

trunk/psLib/src/imageops/psImageCovariance.c

@@ -34 +34 @@
 static float imageCovarianceCalculate(const psKernel *covar, // Original covariance matrix
                                       const psKernel *kernel, // Convolution kernel
-                                      int x, int y            // Coordinates in output covariance matrix
+                                      int x, int y,           // Coordinates in output covariance matrix
+                                      float scale             // Scale to apply
                                       )
 {
@@ -83 +84 @@
     }
 
-    return sum;
+    return scale * sum;
 }
 
@@ -91 +92 @@
     PS_ASSERT_THREAD_JOB_NON_NULL(job, false);
     psAssert(job->args, "No job arguments");
-    psAssert(job->args->n == 5, "Wrong number of job arguments: %ld", job->args->n);
+    psAssert(job->args->n == 6, "Wrong number of job arguments: %ld", job->args->n);
 
     psKernel *out = job->args->data[0]; // Output covariance matrix
@@ -98 +99 @@
     int x = PS_SCALAR_VALUE(job->args->data[3], S32); // x coordinate in output covariance matrix
     int y = PS_SCALAR_VALUE(job->args->data[4], S32); // y coordinate in output covariance matrix
-
-    out->kernel[y][x] = imageCovarianceCalculate(covar, kernel, x, y);
+    float scale = PS_SCALAR_VALUE(job->args->data[5], F32); // Scaling to apply
+
+    out->kernel[y][x] = imageCovarianceCalculate(covar, kernel, x, y, scale);
 
     return true;
@@ -127 +129 @@
 
     // Check for non-finite elements
+    double sumKernel = 0.0, sumKernel2 = 0.0; // Sum of the kernel
     for (int y = kernel->yMin; y <= kernel->yMax; y++) {
         for (int x = kernel->xMin; x <= kernel->xMax; x++) {
@@ -135 +138 @@
                 return NULL;
             }
+            sumKernel += kernel->kernel[y][x];
+            sumKernel2 += PS_SQR(kernel->kernel[y][x]);
         }
     }
@@ -155 +160 @@
     int yMin = kernel->yMin - kernel->yMax + covar->yMin, yMax = kernel->yMax - kernel->yMin + covar->yMax;
     psKernel *out = psKernelAlloc(xMin, xMax, yMin, yMax); // Covariance matrix for output
+    float scale = 1.0 / sumKernel2;          // Scaling to apply
 
     for (int y = yMin; y <= yMax; y++) {
@@ -165 +171 @@
                 PS_ARRAY_ADD_SCALAR(job->args, x, PS_TYPE_S32);
                 PS_ARRAY_ADD_SCALAR(job->args, y, PS_TYPE_S32);
+                PS_ARRAY_ADD_SCALAR(job->args, scale, PS_TYPE_F32);
                 if (!psThreadJobAddPending(job)) {
                     psFree(covar);
@@ -170 +177 @@
                 }
             } else {
-                out->kernel[y][x] = imageCovarianceCalculate(covar, kernel, x, y);
-            }
-        }
-    }
-    psFree(covar);
+                out->kernel[y][x] = imageCovarianceCalculate(covar, kernel, x, y, scale);
+            }
+        }
+    }
 
     if (threaded && !psThreadPoolWait(true)) {
@@ -180 +186 @@
         return false;
     }
+
+    psFree(covar);
 
     return out;
@@ -194 +202 @@
 
     // Check for non-finite elements
+    double sumKernel2 = 0.0; // Sum of the squared kernel
     for (int y = kernel->yMin; y <= kernel->yMax; y++) {
         for (int x = kernel->xMin; x <= kernel->xMax; x++) {
@@ -202 +211 @@
                 return NAN;
             }
+            sumKernel2 += PS_SQR(kernel->kernel[y][x]);
         }
     }
@@ -215 +225 @@
     }
 
-    float factor = imageCovarianceCalculate(covar, kernel, 0, 0); // Covariance factor
+    float scale = 1.0 / sumKernel2;     // Scale to apply
+    float factor = imageCovarianceCalculate(covar, kernel, 0, 0, scale); // Covariance factor
     psFree(covar);
     return factor;
@@ -338 +349 @@
         }
     }
-    psFree(covar);
-
     if (threaded && !psThreadPoolWait(true)) {
         psError(PS_ERR_UNKNOWN, false, "Error waiting for threads.");
         return false;
     }
+    psFree(covar);
 
     return out;
@@ -616 +626 @@
     if (set && !threaded) {
         {
-            psThreadTask *task = psThreadTaskAlloc("PSLIB_IMAGE_COVARIANCE_CALCULATE", 5);
+            psThreadTask *task = psThreadTaskAlloc("PSLIB_IMAGE_COVARIANCE_CALCULATE", 6);
             task->function = &imageCovarianceCalculateThread;
             psThreadTaskAdd(task);

trunk/psLib/src/imageops/psImageInterpolate.c

@@ -427 +427 @@
 
 // Determine the result of the interpolation after all the math has been done
-#define INTERPOLATE_RESULT() \
-    psImageInterpolateStatus status = PS_INTERPOLATE_STATUS_ERROR; /* Status of interpolation */ \
-    *imageValue = sumKernel > 0 ? sumImage / sumKernel : interp->badImage; \
-    if (wantVariance) { \
-        *varianceValue = sumVariance / (sumKernel2 - sumBad); \
-    } \
-    if (sumKernel == 0.0) { \
-        /* No kernel contributions */ \
-        if (haveMask && maskValue) { \
-            *maskValue |= interp->badMask; \
-        } \
-        status = PS_INTERPOLATE_STATUS_BAD; \
-    } else if (sumBad == 0) { \
-        /* Completely good pixel */ \
-        status = PS_INTERPOLATE_STATUS_GOOD; \
-    } else if (sumBad < PS_SQR(interp->poorFrac) * sumKernel2) { \
-        /* Some pixels masked: poor pixel */ \
-        if (haveMask && maskValue) { \
-            *maskValue |= interp->poorMask; \
-        } \
-        status = PS_INTERPOLATE_STATUS_POOR; \
-    } else { \
-        /* Many pixels (or a few important pixels) masked: bad pixel */ \
-        if (haveMask && maskValue) { \
-            *maskValue |= interp->badMask; \
-        } \
-        status = PS_INTERPOLATE_STATUS_BAD; \
-    }
+static psImageInterpolateStatus interpolateResult(const psImageInterpolation *interp,
+                                                  double *imageValue, double *varianceValue,
+                                                  psImageMaskType *maskValue,
+                                                  double sumImage, double sumVariance, double sumBad,
+                                                  double sumKernel, double sumKernel2,
+                                                  bool wantVariance, bool haveMask)
+{
+    *imageValue = sumKernel > 0 ? sumImage / sumKernel : interp->badImage;
+    if (wantVariance) {
+        if (sumBad > 0) {
+            sumVariance *= sumKernel2 / (sumKernel2 - sumBad);
+        }
+        *varianceValue = sumVariance / PS_SQR(sumKernel);
+    }
+    if (sumKernel == 0.0) {
+        // No kernel contributions at all
+        if (haveMask && maskValue) {
+            *maskValue |= interp->badMask;
+        }
+        return PS_INTERPOLATE_STATUS_BAD;
+    }
+    if (sumBad == 0) {
+        // Completely good pixel
+        return PS_INTERPOLATE_STATUS_GOOD;
+    }
+    if (sumBad < PS_SQR(interp->poorFrac) * sumKernel2) {
+        // Some pixels masked: poor pixel
+        if (haveMask && maskValue) {
+            *maskValue |= interp->poorMask;
+        }
+        return PS_INTERPOLATE_STATUS_POOR;
+    }
+    // Many pixels (or a few important pixels) masked: bad pixel
+    if (haveMask && maskValue) {
+        *maskValue |= interp->badMask;
+    }
+    return PS_INTERPOLATE_STATUS_BAD;
+}
 
 // Interpolation engine for separable interpolation kernels
@@ -703 +713 @@
     }
 
-    INTERPOLATE_RESULT();
-
     psFree(xKernelNew);
     psFree(yKernelNew);
@@ -710 +718 @@
     psFree(yKernel2New);
 
-    return status;
+    return interpolateResult(interp, imageValue, varianceValue, maskValue, sumImage, sumVariance, sumBad,
+                             sumKernel, sumKernel2, wantVariance, haveMask);
 }
 
@@ -861 +870 @@
     }
 
-    INTERPOLATE_RESULT();
-
-    return status;
+    return interpolateResult(interp, imageValue, varianceValue, maskValue, sumImage, sumVariance, sumBad,
+                             sumKernel, sumKernel2, wantVariance, haveMask);
 }
 

trunk/psModules/src/imcombine/pmSubtraction.c

@@ -52 +52 @@
 
     // Take the square of the normal kernel
-    double sumVariance = 0.0; // Sum of the variance kernels
     for (int v = yMin; v <= yMax; v++) {
         for (int u = xMin; u <= xMax; u++) {
-            sumVariance += out->kernel[v][u] = PS_SQR(normalKernel->kernel[v][u]);
-        }
-    }
-
-    // Normalise so that the sum of the variance kernel is the square of the sum of the normal kernel
-    // This is required to keep the relative scaling between the image and the variance map
-    psBinaryOp(out->image, out->image, "*", psScalarAlloc(1.0 / sumVariance, PS_TYPE_F32));
-
+            out->kernel[v][u] = PS_SQR(normalKernel->kernel[v][u]);
+        }
+    }
     return out;
 }
@@ -271 +265 @@
 // Convolve an image using FFT
 static void convolveVarianceFFT(psImage *target,// Place the result in here
-                              psImage *variance, // Variance map to convolve
-                              psImage *kernelErr, // Kernel error image
-                              psImage *mask, // Mask image
-                              psImageMaskType maskVal, // Value to mask
-                              const psKernel *kernel, // Kernel by which to convolve
-                              psRegion region,// Region of interest
-                              int size        // Size of (square) kernel
-                              )
+                                psImage *variance, // Variance map to convolve
+                                psImage *kernelErr, // Kernel error image
+                                psImage *mask, // Mask image
+                                psImageMaskType maskVal, // Value to mask
+                                const psKernel *kernel, // Kernel by which to convolve
+                                psRegion region,// Region of interest
+                                int size        // Size of (square) kernel
+                                )
 {
     psRegion border = psRegionSet(region.x0 - size, region.x1 + size,
@@ -348 +342 @@
                                   psImage *image, // Image to convolve
                                   psImage *variance, // Variance map to convolve, or NULL
+                                  const psKernel *covar,               // Covariance, or NULL
                                   psImage *kernelErr, // Kernel error image, or NULL
                                   psImage *subMask, // Subtraction mask
@@ -393 +388 @@
         if (variance) {
             convolveDirect(convVariance, variance, *kernelVariance, region, 0.0, kernels->size);
+        }
+    }
+
+    if (variance && covar) {
+        // Apply covariance factor to variance map, to allow for spatial variation
+        float factor = psImageCovarianceCalculateFactor(*kernelImage, covar); // Factor to apply
+        for (int y = region.y0; y < region.y1; y++) {
+            for (int x = region.x0; x < region.x1; x++) {
+                convVariance->data.F32[y][x] *= factor;
+            }
         }
     }
@@ -1085 +1090 @@
     if (kernels->mode == PM_SUBTRACTION_MODE_1 || kernels->mode == PM_SUBTRACTION_MODE_DUAL) {
         convolveRegion(out1->image, out1->variance, out1->mask, &kernelImage, &kernelVariance,
-                       ro1->image, ro1->variance, kernelErr1, subMask, kernels, polyValues, background,
-                       *region, maskBad, maskPoor, poorFrac, useFFT, false);
+                       ro1->image, ro1->variance, ro1->covariance, kernelErr1, subMask, kernels,
+                       polyValues, background, *region, maskBad, maskPoor, poorFrac, useFFT, false);
     }
     if (kernels->mode == PM_SUBTRACTION_MODE_2 || kernels->mode == PM_SUBTRACTION_MODE_DUAL) {
         convolveRegion(out2->image, out2->variance, out2->mask, &kernelImage, &kernelVariance,
-                       ro2->image, ro2->variance, kernelErr2, subMask, kernels, polyValues, background,
-                       *region, maskBad, maskPoor, poorFrac, useFFT,
+                       ro2->image, ro2->variance, ro2->covariance, kernelErr2, subMask, kernels,
+                       polyValues, background, *region, maskBad, maskPoor, poorFrac, useFFT,
                        kernels->mode == PM_SUBTRACTION_MODE_DUAL);
     }
@@ -1325 +1330 @@
 
     // Calculate covariances
-    // This can be fairly involved, so we only do it for a single instance
-    // Enable threads for covariance calculation, since we're not threading on top of it.
+    // This can be fairly involved, so we only do it for a small number of instances
     float position[NUM_COVAR_POS] = { -1.0, -0.5, 0.0, +0.5, +1.0 }; // Positions for covariance calculations
+    // Enable threads for covariance calculation, since we're not threading on top of it
     oldThreads = psImageCovarianceSetThreads(true);
     if (kernels->mode == PM_SUBTRACTION_MODE_1 || kernels->mode == PM_SUBTRACTION_MODE_DUAL) {
@@ -1342 +1347 @@
         out1->covariance = psImageCovarianceAverage(covars);
         psFree(covars);
+        // Remove covariance factor from covariance, since we've put it in the variance map already
+        float factor = psImageCovarianceFactor(out1->covariance);
+        psBinaryOp(out1->covariance->image, out1->covariance->image, "/", psScalarAlloc(factor, PS_TYPE_F32));
     }
     if (kernels->mode == PM_SUBTRACTION_MODE_2 || kernels->mode == PM_SUBTRACTION_MODE_DUAL) {
@@ -1356 +1364 @@
         out2->covariance = psImageCovarianceAverage(covars);
         psFree(covars);
+        // Remove covariance factor from covariance, since we've put it in the variance map already
+        float factor = psImageCovarianceFactor(out2->covariance);
+        psBinaryOp(out2->covariance->image, out2->covariance->image, "/", psScalarAlloc(factor, PS_TYPE_F32));
     }
     psImageCovarianceSetThreads(oldThreads);

trunk/psphot/src/psphotSignificanceImage.c

@@ -76 +76 @@
     // Calculate correction factor for the covariance produced by the (potentially multiple) smoothing
     psKernel *kernel = psImageSmoothKernel(SIGMA_SMTH, NSIGMA_SMTH); // Kernel used for smoothing
-    float factor = 1.0 / psImageCovarianceCalculateFactor(kernel, readout->covariance);
+    double sum2 = 0.0;                                               // Sum of kernel squared
+    for (int y = kernel->yMin; y <= kernel->yMax; y++) {
+        for (int x = kernel->xMin; x <= kernel->xMax; x++) {
+            sum2 += PS_SQR(kernel->kernel[y][x]);
+        }
+    }
+    float factor = 1.0 / (sum2 * psImageCovarianceCalculateFactor(kernel, readout->covariance));
     psFree(kernel);