Changeset 23594 for branches/cnb_branches/cnb_branch_20090301/ppStack/src/ppStackMatch.c

Timestamp:

Mar 29, 2009, 6:15:31 PM (17 years ago)

Author:

beaumont

Message:

merged with head

Location:

branches/cnb_branches/cnb_branch_20090301

Files:

• . (modified) (1 prop)
• ppStack (modified) (1 prop)
• ppStack/src/ppStackMatch.c (modified) (13 diffs)

branches/cnb_branches/cnb_branch_20090301/ppStack/src/ppStackMatch.c

@@ -163 +163 @@
 
 
-bool ppStackMatch(pmReadout *readout, psArray **regions, psArray **kernels, float *chi2, float *weighting,
-                  psArray *sources, const pmPSF *psf, psRandom *rng, const pmConfig *config)
+bool ppStackMatch(pmReadout *readout, ppStackOptions *options, int index, const pmConfig *config)
 {
     assert(readout);
-    assert(regions && !*regions);
-    assert(kernels && !*kernels);
+    assert(options);
     assert(config);
-    *weighting = 0.0;
 
     psMetadata *recipe = psMetadataLookupMetadata(NULL, config->recipes, PPSTACK_RECIPE); // ppStack recipe
@@ -197 +194 @@
     }
 
-    pmReadout *output = pmReadoutAlloc(NULL); // Output readout, for holding results temporarily
-
-    static int numInput = -1;            // Index of input file
-    numInput++;
+    // Match the PSF
+    if (options->convolve) {
+        pmReadout *conv = pmReadoutAlloc(NULL); // Conv readout, for holding results temporarily
 #ifdef TESTING
-    // Read previously produced kernel
-    if (psMetadataLookupBool(NULL, config->arguments, "PPSTACK.DEBUG.STACK")) {
-        const char *outName = psMetadataLookupStr(NULL, config->arguments, "OUTPUT"); // Output root
-        assert(outName);
-        // Read convolution kernel
-        psString filename = NULL;   // Output filename
-        psStringAppend(&filename, "%s.%d.kernel", outName, numInput);
-        psString resolved = pmConfigConvertFilename(filename, config, false, false); // Resolved filename
-        psFree(filename);
-        psFits *fits = psFitsOpen(resolved, "r"); // FITS file for subtraction kernel
-        psFree(resolved);
-        if (!fits || !pmReadoutReadSubtractionKernels(output, fits)) {
-            psError(PS_ERR_IO, false, "Unable to read previously produced kernel");
+        // This is a hack to use the temporary convolved images and kernel generated previously.
+        // This makes the 'matching' operation much faster, allowing debugging of the stack process easier.
+        // It implicitly assumes the output root name is the same between invocations.
+
+        // Read previously produced kernel
+        if (psMetadataLookupBool(NULL, config->arguments, "PPSTACK.DEBUG.STACK")) {
+            pmFPAfile *file = pmFPAfileSelectSingle(config->files, "PPSTACK.CONV.KERNEL", index);
+            psAssert(file, "Require file");
+
+            pmFPAview *view = pmFPAviewAlloc(0); // View to readout of interest
+            view->chip = view->cell = view->readout = 0;
+            psString filename = pmFPAfileNameFromRule(filerule->rule, file, view); // Filename of interest
+
+            // Read convolution kernel
+            psString resolved = pmConfigConvertFilename(filename, config, false, false); // Resolved filename
+            psFree(filename);
+            psFits *fits = psFitsOpen(resolved, "r"); // FITS file for subtraction kernel
+            psFree(resolved);
+            if (!fits || !pmReadoutReadSubtractionKernels(conv, fits)) {
+                psError(PS_ERR_IO, false, "Unable to read previously produced kernel");
+                psFitsClose(fits);
+                return false;
+            }
             psFitsClose(fits);
-            return false;
-        }
-        psFitsClose(fits);
-
-        // Add in variance factor
-        pmSubtractionKernels *kernels = psMetadataLookupPtr(NULL, output->analysis,
-                                                            PM_SUBTRACTION_ANALYSIS_KERNEL); // Kernels
-        float vf = pmSubtractionVarianceFactor(kernels, 0.0, 0.0, false); // Variance factor
-        psMetadataItem *vfItem = psMetadataLookup(readout->parent->concepts, "CELL.VARFACTOR");
-        if (!isfinite(vf)) {
-            vf = 1.0;
-        }
-        if (isfinite(vfItem->data.F32)) {
-            vfItem->data.F32 *= vf;
-        } else {
-            vfItem->data.F32 = vf;
-        }
-
-        // Read image, mask, variance
-        const char *tempImage = psMetadataLookupStr(NULL, recipe, "TEMP.IMAGE"); // Suffix for image
-        const char *tempMask = psMetadataLookupStr(NULL, recipe, "TEMP.MASK"); // Suffix for mask
-        const char *tempVariance = psMetadataLookupStr(NULL, recipe, "TEMP.VARIANCE"); // Suffix for variance map
-        psString imageName = NULL, maskName = NULL, varianceName = NULL; // Names for convolved images
-        psStringAppend(&imageName, "%s.%d.%s", outName, numInput, tempImage);
-        psStringAppend(&maskName, "%s.%d.%s", outName, numInput, tempMask);
-        psStringAppend(&varianceName, "%s.%d.%s", outName, numInput, tempVariance);
-
-        if (!readImage(&readout->image, imageName, config) || !readImage(&readout->mask, maskName, config) ||
-            !readImage(&readout->variance, varianceName, config)) {
-            psError(PS_ERR_IO, false, "Unable to read previously produced image.");
+
+            // Add in variance factor
+            pmSubtractionKernels *kernels = psMetadataLookupPtr(NULL, conv->analysis,
+                                                                PM_SUBTRACTION_ANALYSIS_KERNEL); // Kernels
+            float vf = pmSubtractionVarianceFactor(kernels, 0.0, 0.0, false); // Variance factor
+            psMetadataItem *vfItem = psMetadataLookup(readout->parent->concepts, "CELL.VARFACTOR");
+            if (!isfinite(vf)) {
+                vf = 1.0;
+            }
+            if (isfinite(vfItem->data.F32)) {
+                vfItem->data.F32 *= vf;
+            } else {
+                vfItem->data.F32 = vf;
+            }
+
+            if (!readImage(&readout->image, options->imageNames->data[index], config) ||
+                !readImage(&readout->mask, options->maskNames->data[index], config) ||
+                !readImage(&readout->variance, options->varianceNames->data[index], config)) {
+                psError(PS_ERR_IO, false, "Unable to read previously produced image.");
+                psFree(imageName);
+                psFree(maskName);
+                psFree(varianceName);
+                return false;
+            }
             psFree(imageName);
             psFree(maskName);
             psFree(varianceName);
-            return false;
-        }
-        psFree(imageName);
-        psFree(maskName);
-        psFree(varianceName);
-
-        psRegion *region = psMetadataLookupPtr(NULL, output->analysis,
-                                               PM_SUBTRACTION_ANALYSIS_REGION); // Convolution region
-
-        pmSubtractionAnalysis(readout->analysis, kernels, region,
-                              readout->image->numCols, readout->image->numRows);
-
-        psKernel *kernel = pmSubtractionKernel(kernels, 0.0, 0.0, false); // Convolution kernel
-        psKernel *covar = psImageCovarianceCalculate(kernel, readout->covariance); // New covariance matrix
-        psFree(readout->covariance);
-        readout->covariance = covar;
-        psFree(kernel);
-
-    } else {
-#endif
-
-        // Normal operations here
-        if (psMetadataLookupBool(&mdok, config->arguments, "HAVE.PSF")) {
-            assert(psf);
-            assert(sources);
+
+            psRegion *region = psMetadataLookupPtr(NULL, conv->analysis,
+                                                   PM_SUBTRACTION_ANALYSIS_REGION); // Convolution region
+
+            pmSubtractionAnalysis(readout->analysis, kernels, region,
+                                  readout->image->numCols, readout->image->numRows);
+
+            psKernel *kernel = pmSubtractionKernel(kernels, 0.0, 0.0, false); // Convolution kernel
+            psKernel *covar = psImageCovarianceCalculate(kernel, readout->covariance); // Covariance matrix
+            psFree(readout->covariance);
+            readout->covariance = covar;
+            psFree(kernel);
+        } else {
+#endif
+
+            // Normal operations here
+            psAssert(options->psf, "Require target PSF");
+            psAssert(options->sourceLists && options->sourceLists->data[index], "Require source list");
 
             int order = psMetadataLookupS32(NULL, ppsub, "SPATIAL.ORDER"); // Spatial polynomial order
@@ -284 +277 @@
             float rej = psMetadataLookupF32(NULL, ppsub, "REJ"); // Rejection threshold
             float sysError = psMetadataLookupF32(NULL, ppsub, "SYS"); // Relative systematic error in kernel
-            pmSubtractionKernelsType type = pmSubtractionKernelsTypeFromString(
-                psMetadataLookupStr(NULL, ppsub, "KERNEL.TYPE")); // Kernel type
+            const char *typeStr = psMetadataLookupStr(NULL, ppsub, "KERNEL.TYPE"); // Kernel type
+            pmSubtractionKernelsType type = pmSubtractionKernelsTypeFromString(typeStr); // Kernel type
             psVector *widths = psMetadataLookupPtr(NULL, ppsub, "ISIS.WIDTHS"); // ISIS Gaussian widths
             psVector *orders = psMetadataLookupPtr(NULL, ppsub, "ISIS.ORDERS"); // ISIS Polynomial orders
@@ -308 +301 @@
             }
 
-#if 0
-            // Testing the normalisation of the fake image
-            {
-                pmReadout *test = pmReadoutAlloc(NULL); // Test readout
-                psArray *sources = psArrayAlloc(1);     // Array of sources
-                pmSource *source = pmSourceAlloc();     // Source
-                sources->data[0] = source;
-                source->peak = pmPeakAlloc(500, 500, 10000, PM_PEAK_LONE);
-                source->psfMag = -13.0;
-                pmReadoutFakeFromSources(test, 1000, 1000, sources, NULL, NULL, psf, 0.1, 0, false, true);
-                float sum = 0.0;
-                for (int y = 0; y < test->image->numRows; y++) {
-                    for (int x = 0; x < test->image->numCols; x++) {
-                        sum += test->image->data.F32[y][x];
-                    }
-                }
-                fprintf(stderr, "Photometry: %f --> %f = -13.0 ???\n", sum, -2.5*log10(sum));
-
-                psFits *fits = psFitsOpen("testphot.fits", "w");
-                psFitsWriteImage(fits, NULL, test->image, 0, NULL);
-                psFitsClose(fits);
-                exit(0);
-            }
-#endif
-
             pmReadout *fake = pmReadoutAlloc(NULL); // Fake readout with target PSF
 
             // For the sake of stamps, remove nearby sources
-            psArray *stampSources = stackSourcesFilter(sources, footprint); // Filtered list of sources
+            psArray *stampSources = stackSourcesFilter(options->sourceLists->data[index],
+                                                       footprint); // Filtered list of sources
 
             if (!pmReadoutFakeFromSources(fake, readout->image->numCols, readout->image->numRows,
-                                          stampSources, NULL, NULL, psf, NAN, footprint + size,
+                                          stampSources, NULL, NULL, options->psf, NAN, footprint + size,
                                           false, true)) {
                 psError(PS_ERR_UNKNOWN, false, "Unable to generate fake image with target PSF.");
                 psFree(fake);
                 psFree(optWidths);
-                psFree(output);
+                psFree(conv);
                 return false;
             }
@@ -357 +326 @@
                 pmHDU *hdu = pmHDUFromCell(readout->parent);
                 psString name = NULL;
-                psStringAppend(&name, "fake_%03d.fits", numInput);
+                psStringAppend(&name, "fake_%03d.fits", index);
                 pmStackVisualPlotTestImage(fake->image, name);
                 psFits *fits = psFitsOpen(name, "w");
@@ -367 +336 @@
                 pmHDU *hdu = pmHDUFromCell(readout->parent);
                 psString name = NULL;
-                psStringAppend(&name, "real_%03d.fits", numInput);
+                psStringAppend(&name, "real_%03d.fits", index);
                 pmStackVisualPlotTestImage(readout->image, name);
                 psFits *fits = psFitsOpen(name, "w");
@@ -381 +350 @@
 
             // Do the image matching
-            if (!pmSubtractionMatch(output, NULL, readout, fake, footprint, stride, regionSize, spacing,
-                                    threshold, stampSources, stampsName, type, size, order, widths, orders,
-                                    inner, ringsOrder, binning, penalty, optimum, optWidths, optOrder,
-                                    optThresh, iter, rej, sysError, maskVal, maskBad, maskPoor, poorFrac,
-                                    badFrac, PM_SUBTRACTION_MODE_1)) {
-                psError(PS_ERR_UNKNOWN, false, "Unable to match images.");
-                psFree(fake);
-                psFree(optWidths);
-                psFree(stampSources);
-                psFree(output);
-                return false;
+            pmSubtractionKernels *kernel = psMetadataLookupPtr(&mdok, readout->analysis,
+                                                               PM_SUBTRACTION_ANALYSIS_KERNEL); // Conv kernel
+            if (kernel) {
+                if (!pmSubtractionMatchPrecalc(conv, NULL, readout, fake, readout->analysis,
+                                               stride, sysError, maskVal, maskBad, maskPoor,
+                                               poorFrac, badFrac)) {
+                    psError(PS_ERR_UNKNOWN, false, "Unable to convolve images.");
+                    psFree(fake);
+                    psFree(optWidths);
+                    psFree(stampSources);
+                    psFree(conv);
+                    return false;
+                }
+            } else {
+                if (!pmSubtractionMatch(conv, NULL, readout, fake, footprint, stride, regionSize, spacing,
+                                        threshold, stampSources, stampsName, type, size, order, widths,
+                                        orders, inner, ringsOrder, binning, penalty,
+                                        optimum, optWidths, optOrder, optThresh, iter, rej, sysError,
+                                        maskVal, maskBad, maskPoor, poorFrac, badFrac,
+                                        PM_SUBTRACTION_MODE_1)) {
+                    psError(PS_ERR_UNKNOWN, false, "Unable to match images.");
+                    psFree(fake);
+                    psFree(optWidths);
+                    psFree(stampSources);
+                    psFree(conv);
+                    return false;
+                }
             }
 
@@ -398 +383 @@
                 pmHDU *hdu = pmHDUFromCell(readout->parent);
                 psString name = NULL;
-                psStringAppend(&name, "conv_%03d.fits", numInput);
-                pmStackVisualPlotTestImage(output->image, name);
+                psStringAppend(&name, "conv_%03d.fits", index);
+                pmStackVisualPlotTestImage(conv->image, name);
                 psFits *fits = psFitsOpen(name, "w");
                 psFree(name);
-                psFitsWriteImage(fits, hdu->header, output->image, 0, NULL);
+                psFitsWriteImage(fits, hdu->header, conv->image, 0, NULL);
                 psFitsClose(fits);
             }
@@ -408 +393 @@
                 pmHDU *hdu = pmHDUFromCell(readout->parent);
                 psString name = NULL;
-                psStringAppend(&name, "diff_%03d.fits", numInput);
+                psStringAppend(&name, "diff_%03d.fits", index);
                 pmStackVisualPlotTestImage(fake->image, name);
                 psFits *fits = psFitsOpen(name, "w");
                 psFree(name);
-                psBinaryOp(fake->image, output->image, "-", fake->image);
+                psBinaryOp(fake->image, conv->image, "-", fake->image);
                 psFitsWriteImage(fits, hdu->header, fake->image, 0, NULL);
                 psFitsClose(fits);
@@ -428 +413 @@
             // Set the variance factor
             psMetadataItem *vfItem = psMetadataLookup(readout->parent->concepts, "CELL.VARFACTOR");
-            float vf = psMetadataLookupF32(NULL, output->analysis, PM_SUBTRACTION_ANALYSIS_VARFACTOR_1);
+            float vf = psMetadataLookupF32(NULL, conv->analysis, PM_SUBTRACTION_ANALYSIS_VARFACTOR_1);
             if (!isfinite(vf)) {
                 vf = 1.0;
@@ -443 +428 @@
             psFree(readout->variance);
             psFree(readout->covariance);
-            readout->image  = psMemIncrRefCounter(output->image);
-            readout->mask   = psMemIncrRefCounter(output->mask);
-            readout->variance = psMemIncrRefCounter(output->variance);
-            readout->covariance = psImageCovarianceTruncate(output->covariance, COVAR_FRAC);
-        } else {
-            // Fake the convolution
-            psRegion *region = psRegionAlloc(0, readout->image->numCols - 1, 0, readout->image->numRows - 1);
-            psMetadataAddPtr(output->analysis, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_REGION,
-                             PS_DATA_REGION | PS_META_DUPLICATE_OK, "Fake subtraction region", region);
-            psFree(region);
-            pmSubtractionKernels *kernels = pmSubtractionKernelsPOIS(FAKE_SIZE, 0, penalty,
-                                                                     PM_SUBTRACTION_MODE_1);
-            // Set solution to delta function
-            kernels->solution1 = psVectorAlloc(kernels->num + 2, PS_TYPE_F64);
-            psVectorInit(kernels->solution1, 0.0);
-            int normIndex = PM_SUBTRACTION_INDEX_NORM(kernels); // Index for normalisation
-            kernels->solution1->data.F64[normIndex] = 1.0;
-            psMetadataAddPtr(output->analysis, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_KERNEL,
-                             PS_DATA_UNKNOWN | PS_META_DUPLICATE_OK, "Fake subtraction kernel", kernels);
-            psFree(kernels);
-        }
-
+            readout->image  = psMemIncrRefCounter(conv->image);
+            readout->mask   = psMemIncrRefCounter(conv->mask);
+            readout->variance = psMemIncrRefCounter(conv->variance);
+            readout->covariance = psImageCovarianceTruncate(conv->covariance, COVAR_FRAC);
 #ifdef TESTING
-        // Write convolution kernel
+        }
+#endif
+
+        // Extract the regions and solutions used in the image matching
+        // This stops them from being freed when we iterate back up the FPA
+        psArray *regions = options->regions->data[index] = psArrayAllocEmpty(ARRAY_BUFFER); // Match regions
         {
-            const char *outName = psMetadataLookupStr(NULL, config->arguments, "OUTPUT"); // Output root
-            assert(outName);
-
-            psString filename = NULL;   // Output filename
-            psStringAppend(&filename, "%s.%d.kernel", outName, numInput);
-            psString resolved = pmConfigConvertFilename(filename, config, true, false); // Resolved filename
-            psFree(filename);
-            psFits *fits = psFitsOpen(resolved, "w"); // FITS file for subtraction kernel
-            psFree(resolved);
-            pmReadoutWriteSubtractionKernels(output, fits);
-            psFitsClose(fits);
-        }
-    }
-#endif
-
-    readout->analysis = psMetadataCopy(readout->analysis, output->analysis);
-
-// Extract the regions and solutions used in the image matching
-// This stops them from being freed when we iterate back up the FPA
-    *regions = psArrayAllocEmpty(ARRAY_BUFFER); // Array of regions
-    {
-        psString regex = NULL;          // Regular expression
-        psStringAppend(&regex, "^%s$", PM_SUBTRACTION_ANALYSIS_REGION);
-        psMetadataIterator *iter = psMetadataIteratorAlloc(output->analysis, PS_LIST_HEAD, regex); // Iterator
-        psFree(regex);
-        psMetadataItem *item = NULL;// Item from iteration
-        while ((item = psMetadataGetAndIncrement(iter))) {
-            assert(item->type == PS_DATA_REGION);
-            *regions = psArrayAdd(*regions, ARRAY_BUFFER, item->data.V);
-        }
-        psFree(iter);
-    }
-    *kernels = psArrayAllocEmpty(ARRAY_BUFFER); // Array of kernels
-    {
-        psString regex = NULL;          // Regular expression
-        psStringAppend(&regex, "^%s$", PM_SUBTRACTION_ANALYSIS_KERNEL);
-        psMetadataIterator *iter = psMetadataIteratorAlloc(output->analysis, PS_LIST_HEAD, regex); // Iterator
-        psFree(regex);
-        psMetadataItem *item = NULL;// Item from iteration
-        while ((item = psMetadataGetAndIncrement(iter))) {
-            assert(item->type == PS_DATA_UNKNOWN);
-            // Set the normalisation dimensions, since these will be otherwise unavailable when reading the
-            // images by scans.
-            pmSubtractionKernels *kernel = item->data.V; // Kernel used in subtraction
-            kernel->numCols = readout->image->numCols;
-            kernel->numRows = readout->image->numRows;
-
-            *kernels = psArrayAdd(*kernels, ARRAY_BUFFER, kernel);
-        }
-        psFree(iter);
-    }
-    assert((*regions)->n == (*kernels)->n);
-
-    // Record chi^2
-    {
-        *chi2 = 0.0;
-        int num = 0;                    // Number of measurements of chi^2
-        psString regex = NULL;          // Regular expression
-        psStringAppend(&regex, "^%s$", PM_SUBTRACTION_ANALYSIS_KERNEL);
-        psMetadataIterator *iter = psMetadataIteratorAlloc(output->analysis, PS_LIST_HEAD, regex); // Iterator
-        psFree(regex);
-        psMetadataItem *item = NULL;// Item from iteration
-        while ((item = psMetadataGetAndIncrement(iter))) {
-            assert(item->type == PS_DATA_UNKNOWN);
-            pmSubtractionKernels *kernels = item->data.V; // Convolution kernels
-            *chi2 += kernels->mean;
-            num++;
-        }
-        psFree(iter);
-        *chi2 /= psImageCovarianceFactor(readout->covariance) * num;
-    }
-
-    // Reject image completely if the maximum deconvolution fraction exceeds the limit
-    float deconv = psMetadataLookupF32(NULL, output->analysis,
-                                       PM_SUBTRACTION_ANALYSIS_DECONV_MAX); // Maximum deconvolution fraction
-    if (deconv > deconvLimit) {
-        psWarning("Maximum deconvolution fraction (%f) exceeds limit (%f) --- rejecting\n",
-                  deconv, deconvLimit);
-        psFree(output);
-        return NULL;
+            psString regex = NULL;          // Regular expression
+            psStringAppend(&regex, "^%s$", PM_SUBTRACTION_ANALYSIS_REGION);
+            psMetadataIterator *iter = psMetadataIteratorAlloc(conv->analysis, PS_LIST_HEAD, regex);
+            psFree(regex);
+            psMetadataItem *item = NULL;// Item from iteration
+            while ((item = psMetadataGetAndIncrement(iter))) {
+                assert(item->type == PS_DATA_REGION);
+                regions = psArrayAdd(regions, ARRAY_BUFFER, item->data.V);
+            }
+            psFree(iter);
+        }
+        psArray *kernels = options->kernels->data[index] = psArrayAllocEmpty(ARRAY_BUFFER); // Match kernels
+        {
+            psString regex = NULL;          // Regular expression
+            psStringAppend(&regex, "^%s$", PM_SUBTRACTION_ANALYSIS_KERNEL);
+            psMetadataIterator *iter = psMetadataIteratorAlloc(conv->analysis, PS_LIST_HEAD, regex);
+            psFree(regex);
+            psMetadataItem *item = NULL;// Item from iteration
+            while ((item = psMetadataGetAndIncrement(iter))) {
+                assert(item->type == PS_DATA_UNKNOWN);
+                // Set the normalisation dimensions, since these will be otherwise unavailable when reading
+                // the images by scans.
+                pmSubtractionKernels *kernel = item->data.V; // Kernel used in subtraction
+                kernel->numCols = readout->image->numCols;
+                kernel->numRows = readout->image->numRows;
+
+                kernels = psArrayAdd(kernels, ARRAY_BUFFER, kernel);
+            }
+            psFree(iter);
+        }
+        psAssert((regions)->n == (kernels)->n, "Number of match regions and kernels should match");
+
+        // Record chi^2
+        {
+            double sum = 0.0;           // Sum of chi^2
+            int num = 0;                // Number of measurements of chi^2
+            psString regex = NULL;      // Regular expression
+            psStringAppend(&regex, "^%s$", PM_SUBTRACTION_ANALYSIS_KERNEL);
+            psMetadataIterator *iter = psMetadataIteratorAlloc(conv->analysis, PS_LIST_HEAD, regex);
+            psFree(regex);
+            psMetadataItem *item = NULL;// Item from iteration
+            while ((item = psMetadataGetAndIncrement(iter))) {
+                assert(item->type == PS_DATA_UNKNOWN);
+                pmSubtractionKernels *kernels = item->data.V; // Convolution kernels
+                sum += kernels->mean;
+                num++;
+            }
+            psFree(iter);
+            options->matchChi2->data.F32[index] = sum / (psImageCovarianceFactor(readout->covariance) * num);
+        }
+
+        // Reject image completely if the maximum deconvolution fraction exceeds the limit
+        float deconv = psMetadataLookupF32(NULL, conv->analysis,
+                                           PM_SUBTRACTION_ANALYSIS_DECONV_MAX); // Max deconvolution fraction
+        if (deconv > deconvLimit) {
+            psWarning("Maximum deconvolution fraction (%f) exceeds limit (%f) --- rejecting\n",
+                      deconv, deconvLimit);
+            psFree(conv);
+            return NULL;
+        }
+
+        readout->analysis = psMetadataCopy(readout->analysis, conv->analysis);
+
+        psFree(conv);
+    } else {
+        // Match the normalisation
+        float norm = powf(10.0, -0.4 * options->norm->data.F32[index]); // Normalisation
+        psBinaryOp(readout->image, readout->image, "*", psScalarAlloc(norm, PS_TYPE_F32));
+        psBinaryOp(readout->variance, readout->variance, "*", psScalarAlloc(PS_SQR(norm), PS_TYPE_F32));
     }
 
     // Ensure the background value is zero
     psStats *bg = psStatsAlloc(PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV); // Statistics for background
+    psRandom *rng = psRandomAlloc(PS_RANDOM_TAUS); // Random number generator
     if (!psImageBackground(bg, NULL, readout->image, readout->mask, maskVal | maskBad, rng)) {
         psWarning("Can't measure background for image.");
@@ -565 +527 @@
     if (!psImageBackground(bg, NULL, readout->variance, readout->mask, maskVal | maskBad, rng)) {
         psError(PS_ERR_UNKNOWN, false, "Can't measure mean variance for image.");
-        psFree(output);
+        psFree(rng);
+        psFree(bg);
         return false;
     }
-    *weighting = 1.0 / (psStatsGetValue(bg, PS_STAT_ROBUST_MEDIAN) *
-                        psImageCovarianceFactor(readout->covariance));
+    options->weightings->data.F32[index] = 1.0 / (psStatsGetValue(bg, PS_STAT_ROBUST_MEDIAN) *
+                                                  psImageCovarianceFactor(readout->covariance));
     psMetadataAddF32(readout->analysis, PS_LIST_TAIL, "PPSTACK.WEIGHTING", 0,
-                     "Weighting by 1/noise^2 for stack", *weighting);
-
+                     "Weighting by 1/noise^2 for stack", options->weightings->data.F32[index]);
+
+    psFree(rng);
     psFree(bg);
-
-#if 0
-#define RADIUS 10                       // Radius of photometry
-#define MIN_ERR 0.05                    // Minimum photometric error, mag
-#define MAX_MAG -13                     // Maximum magnitude for source
-
-    // Ensure the normalisation is correct
-    // XXX Ideally, would like to do proper PSF-fit photometry, but will settle for aperture photometry
-    int numSources = sources->n;        // Number of sources
-    psVector *ratio = psVectorAlloc(numSources, PS_TYPE_F32); // Flux ratios for sources
-    psVector *ratioMask = psVectorAlloc(numSources, PS_TYPE_MASK); // Mask for flux ratios
-    psVectorInit(ratioMask, 0xFF);
-    psImage *image = readout->image;    // Image of interest
-    psImage *mask = readout->mask;      // Mask of interest
-    int numCols = image->numCols, numRows = image->numRows; // Size of image
-    for (int i = 0; i < numSources; i++) {
-        pmSource *source = sources->data[i]; // Source of interest
-        if (!source || source->mode & SOURCE_MASK || !isfinite(source->psfMag) || !isfinite(source->errMag) ||
-            source->errMag > MIN_ERR || source->psfMag > MAX_MAG) {
-            continue;
-        }
-
-        float xSrc, ySrc;              // Source coordinates
-        coordsFromSource(&xSrc, &ySrc, source);
-        int xMin = PS_MAX(0, xSrc - RADIUS), xMax = PS_MIN(numCols - 1, xSrc + RADIUS); // Bounds in x
-        int yMin = PS_MAX(0, ySrc - RADIUS), yMax = PS_MIN(numRows - 1, ySrc + RADIUS); // Bounds in y
-        int numPix = 0;                 // Number of pixels
-        float sum = 0.0;                // Sum of pixels
-        for (int y = yMin; y <= yMax; y++) {
-            for (int x = xMin; x <= xMax; x++) {
-                if (mask->data.PS_TYPE_MASK_DATA[y][x] & maskBad) {
-                    continue;
-                }
-                sum += image->data.F32[y][x];
-                numPix++;
-            }
-        }
-        if (sum >= 0 && numPix > 0) {
-            float mag = -2.5 * log10(sum * M_PI * PS_SQR(RADIUS) / numPix); // Instrumental magnitude
-            ratio->data.F32[i] = mag - source->psfMag;
-            ratioMask->data.PS_TYPE_MASK_DATA[i] = 0;
-        }
-    }
-
-    psStats *stats = psStatsAlloc(PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV); // Statistics
-    if (!psVectorStats(stats, ratio, NULL, ratioMask, 0xFF)) {
-        psWarning("Unable to measure normalisation --- assuming correct.");
-    } else {
-        psLogMsg("ppStack", PS_LOG_INFO, "Renormalising image by %f (+/- %f) mag\n",
-                 stats->robustMedian, stats->robustStdev);
-        float norm = powf(10.0, -0.4 * stats->robustMedian); // Normalisation to apply
-        psBinaryOp(image, image, "*", psScalarAlloc(norm, PS_TYPE_F32));
-    }
-    psFree(stats);
-    psFree(ratio);
-    psFree(ratioMask);
-#endif
 
 #ifdef TESTING
@@ -636 +543 @@
         pmHDU *hdu = pmHDUFromCell(readout->parent);
         psString name = NULL;
-        psStringAppend(&name, "convolved_%03d.fits", numInput);
-        pmStackVisualPlotTestImage(output->image, name);
+        psStringAppend(&name, "convolved_%03d.fits", index);
+        pmStackVisualPlotTestImage(readout->image, name);
         psFits *fits = psFitsOpen(name, "w");
         psFree(name);
-        psFitsWriteImage(fits, hdu->header, output->image, 0, NULL);
+        psFitsWriteImage(fits, hdu->header, readout->image, 0, NULL);
         psFitsClose(fits);
     }
 #endif
-
-    psFree(output);
 
     return true;