Changeset 23573

Timestamp:

Mar 27, 2009, 11:49:17 AM (17 years ago)

Author:

Paul Price

Message:

Allow convolution of images to be optional when stacking.

Location:

trunk

Files:

• ippconfig/recipes/ppStack.config (modified) (1 diff)
• ppStack/src/ppStack.h (modified) (3 diffs)
• ppStack/src/ppStackCombineInitial.c (modified) (3 diffs)
• ppStack/src/ppStackConvolve.c (modified) (5 diffs)
• ppStack/src/ppStackMatch.c (modified) (15 diffs)
• ppStack/src/ppStackOptions.c (modified) (3 diffs)
• ppStack/src/ppStackOptions.h (modified) (2 diffs)
• ppStack/src/ppStackPrepare.c (modified) (3 diffs)
• ppStack/src/ppStackReject.c (modified) (3 diffs)
• ppStack/src/ppStackSetup.c (modified) (3 diffs)
• ppStack/src/ppStackSources.c (modified) (13 diffs)

trunk/ippconfig/recipes/ppStack.config

@@ -1 +1 @@
 # Recipe configuration for ppStack (image combination)
 
+CONVOLVE        BOOL    TRUE            # Convolve images when stacking?
 ITER            S32     1               # Number of rejection iterations
 COMBINE.REJ     F32     3.0             # Rejection threshold in combination (sigma)

trunk/ppStack/src/ppStack.h

@@ -7 +7 @@
 #include <pslib.h>
 #include <psmodules.h>
+
+#include "ppStackOptions.h"
 
 // Mask values for inputs
@@ -106 +108 @@
 /// Convolve image to match specified seeing
 bool ppStackMatch(pmReadout *readout,   // Readout to be convolved; replaced with output
-                  psArray **regions,    // Array of regions used in each PSF matching, returned
-                  psArray **kernels,    // Array of kernels used in each PSF matching, returned
-                  float *chi2,          // Chi^2 from the stamps
-                  float *weighting,     // Stack weighting (1/noise^2)
-                  psArray *sources,     // Array of sources
-                  const pmPSF *psf,     // Target PSF
-                  psRandom *rng,        // Random number generator
+                  ppStackOptions *options, // Options for stacking
+                  int index,            // Index of image to match
                   const pmConfig *config // Configuration
     );
@@ -120 +117 @@
 ///
 /// Corrects the source PSF photometry to a common system.  Return the sum of the exposure times.
-float ppStackSourcesTransparency(const psArray *sourceLists, // Sources for each input
-                                 psVector *inputMask, // Indicates bad input
-                                 const pmFPAview *view, // View to readout
-                                 const pmConfig *config // Configuration
+bool ppStackSourcesTransparency(ppStackOptions *options, // Stacking options
+                                const pmFPAview *view, // View to readout
+                                const pmConfig *config // Configuration
     );
 

trunk/ppStack/src/ppStackCombineInitial.c

@@ -16 +16 @@
     psAssert(config, "Require configuration");
 
-    psTimerStart("PPSTACK_FINAL");
+    if (!options->convolve) {
+        // No need to do initial combination when we haven't convolved
+        return true;
+    }
+
+    psTimerStart("PPSTACK_INITIAL");
 
     psMetadata *recipe = psMetadataLookupMetadata(NULL, config->recipes, PPSTACK_RECIPE); // ppStack recipe
@@ -69 +74 @@
         }
 
-        // call: ppStackReadoutInitial(config, outRO, readouts, subRegions, subKernels)
         psThreadJob *job = psThreadJobAlloc("PPSTACK_INITIAL_COMBINE"); // Job to start
         psArrayAdd(job->args, 1, thread);
@@ -128 +132 @@
 
     if (options->stats) {
-        psMetadataAddF32(options->stats, PS_LIST_TAIL, "TIME_FINAL", 0,
-                         "Time to make final stack", psTimerMark("PPSTACK_FINAL"));
+        psMetadataAddF32(options->stats, PS_LIST_TAIL, "TIME_INITIAL", 0,
+                         "Time to make final stack", psTimerMark("PPSTACK_INITIAL"));
     }
 

trunk/ppStack/src/ppStackConvolve.c

@@ -30 +30 @@
     int num = options->num;             // Number of inputs
     options->cells = psArrayAlloc(num); // Cells for convolved images --- a handle for reading again
-    options->subKernels = psArrayAlloc(num); // Subtraction kernels --- required in the stacking
-    options->subRegions = psArrayAlloc(num); // Subtraction regions --- required in the stacking
+    options->kernels = psArrayAlloc(num); // PSF-matching kernels --- required in the stacking
+    options->regions = psArrayAlloc(num); // PSF-matching regions --- required in the stacking
     int numGood = 0;                    // Number of good frames
     options->numCols = 0;
@@ -38 +38 @@
     psVectorInit(options->matchChi2, NAN);
     options->weightings = psVectorAlloc(num, PS_TYPE_F32); // Combination weightings for images (1/noise^2)
-    psVectorInit(options->weightings, NAN);
+    psVectorInit(options->weightings, 0.0);
     options->covariances = psArrayAlloc(num); // Covariance matrices
 
@@ -78 +78 @@
 
         // Background subtraction, scaling and normalisation is performed automatically by the image matching
-        psArray *regions = NULL, *kernels = NULL; // Regions and kernels used in subtraction
         psTimerStart("PPSTACK_MATCH");
-
-        if (!ppStackMatch(readout, &regions, &kernels, &options->matchChi2->data.F32[i],
-                          &options->weightings->data.F32[i], options->sourceLists->data[i],
-                          options->psf, rng, config)) {
+        if (!ppStackMatch(readout, options, i, config)) {
             psErrorStackPrint(stderr, "Unable to match image %d --- ignoring.", i);
-            options->inputMask->data.PS_TYPE_VECTOR_MASK_DATA[i] = PPSTACK_MASK_MATCH;
+            options->inputMask->data.PS_TYPE_VECTOR_MASK_DATA[i] |= PPSTACK_MASK_MATCH;
             psErrorClear();
             continue;
@@ -92 +88 @@
 
         if (options->stats) {
-            pmSubtractionKernels *kernels = psMetadataLookupPtr(NULL, readout->analysis,
-                                                                PM_SUBTRACTION_ANALYSIS_KERNEL);// Conv kernel
             psMetadataAddF32(options->stats, PS_LIST_TAIL, "TIME_MATCH", PS_META_DUPLICATE_OK,
                              "Time to match PSF", psTimerMark("PPSTACK_MATCH"));
-            psMetadataAddF32(options->stats, PS_LIST_TAIL, "STAMP.MEAN", PS_META_DUPLICATE_OK,
-                             "Mean deviation for stamps", kernels->mean);
-            psMetadataAddF32(options->stats, PS_LIST_TAIL, "STAMP.RMS", PS_META_DUPLICATE_OK,
-                             "RMS deviation for stamps", kernels->rms);
-            psMetadataAddF32(options->stats, PS_LIST_TAIL, "STAMP.NUM", PS_META_DUPLICATE_OK,
-                             "Number of stamps", kernels->numStamps);
-            float deconv = psMetadataLookupF32(NULL, readout->analysis,
-                                               PM_SUBTRACTION_ANALYSIS_DECONV_MAX); // Deconvolution fraction
-            psMetadataAddF32(options->stats, PS_LIST_TAIL, "KERNEL.DECONV", PS_META_DUPLICATE_OK,
-                             "Deconvolution fraction for kernel", deconv);
             psMetadataAddF32(options->stats, PS_LIST_TAIL, "PPSTACK.WEIGHTING", PS_META_DUPLICATE_OK,
                              "Weighting for image", options->weightings->data.F32[i]);
+
+            if (options->convolve) {
+                // Pull parameters out of convolution kernel
+                pmSubtractionKernels *kernels = psMetadataLookupPtr(NULL, readout->analysis,
+                                                                    PM_SUBTRACTION_ANALYSIS_KERNEL);
+                psMetadataAddF32(options->stats, PS_LIST_TAIL, "STAMP.MEAN", PS_META_DUPLICATE_OK,
+                                 "Mean deviation for stamps", kernels->mean);
+                psMetadataAddF32(options->stats, PS_LIST_TAIL, "STAMP.RMS", PS_META_DUPLICATE_OK,
+                                 "RMS deviation for stamps", kernels->rms);
+                psMetadataAddF32(options->stats, PS_LIST_TAIL, "STAMP.NUM", PS_META_DUPLICATE_OK,
+                                 "Number of stamps", kernels->numStamps);
+                float deconv = psMetadataLookupF32(NULL, readout->analysis,
+                                                   PM_SUBTRACTION_ANALYSIS_DECONV_MAX);
+                psMetadataAddF32(options->stats, PS_LIST_TAIL, "KERNEL.DECONV", PS_META_DUPLICATE_OK,
+                                 "Deconvolution fraction for kernel", deconv);
+            }
         }
         psLogMsg("ppStack", PS_LOG_INFO, "Time to match image %d: %f sec", i, psTimerClear("PPSTACK_MATCH"));
 
-        options->subRegions->data[i] = regions;
-        options->subKernels->data[i] = kernels;
-
         // Write the temporary convolved files
-        pmHDU *hdu = readout->parent->parent->parent->hdu; // HDU for convolved image
-        assert(hdu);
-        ppStackWriteImage(options->imageNames->data[i], hdu->header, readout->image, config);
-        psMetadata *maskHeader = psMetadataCopy(NULL, hdu->header); // Copy of header, for mask
-        pmConfigMaskWriteHeader(config, maskHeader);
-        ppStackWriteImage(options->maskNames->data[i], maskHeader, readout->mask, config);
-        psFree(maskHeader);
-        psImageCovarianceTransfer(readout->variance, readout->covariance);
-        ppStackWriteImage(options->varianceNames->data[i], hdu->header, readout->variance, config);
+        if (options->convolve) {
+            pmHDU *hdu = readout->parent->parent->parent->hdu; // HDU for convolved image
+            assert(hdu);
+            ppStackWriteImage(options->imageNames->data[i], hdu->header, readout->image, config);
+            psMetadata *maskHeader = psMetadataCopy(NULL, hdu->header); // Copy of header, for mask
+            pmConfigMaskWriteHeader(config, maskHeader);
+            ppStackWriteImage(options->maskNames->data[i], maskHeader, readout->mask, config);
+            psFree(maskHeader);
+            psImageCovarianceTransfer(readout->variance, readout->covariance);
+            ppStackWriteImage(options->varianceNames->data[i], hdu->header, readout->variance, config);
 #ifdef TESTING
-        {
-            psString name = NULL;
-            psStringAppend(&name, "covariance_%d.fits", i);
-            ppStackWriteImage(name, hdu->header, readout->covariance->image, config);
-            pmStackVisualPlotTestImage(readout->covariance->image, name);
-            psFree(name);
-        }
+            {
+                psString name = NULL;
+                psStringAppend(&name, "covariance_%d.fits", i);
+                ppStackWriteImage(name, hdu->header, readout->covariance->image, config);
+                pmStackVisualPlotTestImage(readout->covariance->image, name);
+                psFree(name);
+            }
 #endif
+        }
 
         pmCell *inCell = readout->parent; // Input cell
@@ -152 +151 @@
     psFree(rng);
 
+    psFree(options->norm); options->norm = NULL;
     psFree(options->sourceLists); options->sourceLists = NULL;
     psFree(options->psf); options->psf = NULL;

trunk/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");
@@ -384 +353 @@
                                                                PM_SUBTRACTION_ANALYSIS_KERNEL); // Conv kernel
             if (kernel) {
-                if (!pmSubtractionMatchPrecalc(output, NULL, readout, fake, readout->analysis,
+                if (!pmSubtractionMatchPrecalc(conv, NULL, readout, fake, readout->analysis,
                                                stride, sysError, maskVal, maskBad, maskPoor,
                                                poorFrac, badFrac)) {
@@ -391 +360 @@
                     psFree(optWidths);
                     psFree(stampSources);
-                    psFree(output);
+                    psFree(conv);
                     return false;
                 }
             } else {
-                if (!pmSubtractionMatch(output, NULL, readout, fake, footprint, stride, regionSize, spacing,
+                if (!pmSubtractionMatch(conv, NULL, readout, fake, footprint, stride, regionSize, spacing,
                                         threshold, stampSources, stampsName, type, size, order, widths,
                                         orders, inner, ringsOrder, binning, penalty,
@@ -405 +374 @@
                     psFree(optWidths);
                     psFree(stampSources);
-                    psFree(output);
+                    psFree(conv);
                     return false;
                 }
@@ -414 +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);
             }
@@ -424 +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);
@@ -444 +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;
@@ -459 +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.");
@@ -581 +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
@@ -652 +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;

trunk/ppStack/src/ppStackOptions.c

@@ -18 +18 @@
     psFree(options->inputMask);
     psFree(options->sourceLists);
+    psFree(options->norm);
     psFree(options->cells);
-    psFree(options->subKernels);
-    psFree(options->subRegions);
+    psFree(options->kernels);
+    psFree(options->regions);
     psFree(options->matchChi2);
     psFree(options->weightings);
@@ -35 +36 @@
     psMemSetDeallocator(options, (psFreeFunc)stackOptionsFree);
 
+    options->convolve = true;
     options->stats = NULL;
     options->statsFile = NULL;
@@ -44 +46 @@
     options->inputMask = NULL;
     options->sourceLists = NULL;
+    options->norm = NULL;
     options->cells = NULL;
-    options->subKernels = NULL;
-    options->subRegions = NULL;
+    options->kernels = NULL;
+    options->regions = NULL;
     options->numCols = 0;
     options->numRows = 0;

trunk/ppStack/src/ppStackOptions.h

@@ -8 +8 @@
 typedef struct {
     // Setup
+    bool convolve;                      // Convolve images?
     psMetadata *stats;                  // Statistics for output
     FILE *statsFile;                    // File to which to write statistics
@@ -17 +18 @@
     psVector *inputMask;                // Mask for inputs
     psArray *sourceLists;               // Individual lists of sources for matching
+    psVector *norm;                     // Normalisation for each image
     // Convolve
     psArray *cells;                     // Cells for convolved images --- a handle for reading again
-    psArray *subKernels;                // Subtraction kernels --- required in the stacking
-    psArray *subRegions;                // Subtraction regions --- required in the stacking
+    psArray *kernels;                   // PSF-matching kernels --- required in the stacking
+    psArray *regions;                   // PSF-matching regions --- required in the stacking
     int numCols, numRows;               // Size of image
     psVector *matchChi2;                // chi^2 for stamps from matching

trunk/ppStack/src/ppStackPrepare.c

@@ -127 +127 @@
     psVectorInit(options->inputMask, 0);
 
-    if (psMetadataLookupBool(NULL, config->arguments, "HAVE.PSF")) {
-        pmFPAfileActivate(config->files, false, NULL);
-        ppStackFileActivation(config, PPSTACK_FILES_PREPARE, true);
-        pmFPAview *view = ppStackFilesIterateDown(config);
-        if (!view) {
-            return false;
-        }
-
-        // Generate list of PSFs
-        psMetadataIterator *fileIter = psMetadataIteratorAlloc(config->files, PS_LIST_HEAD,
-                                                               "^PPSTACK.INPUT$"); // Iterator
-        psMetadataItem *fileItem; // Item from iteration
-        psArray *psfs = psArrayAlloc(num); // PSFs for PSF envelope
-        int numCols = 0, numRows = 0;   // Size of image
-        int index = 0;                  // Index for file
-        while ((fileItem = psMetadataGetAndIncrement(fileIter))) {
-            assert(fileItem->type == PS_DATA_UNKNOWN);
-            pmFPAfile *inputFile = fileItem->data.V; // An input file
+    pmFPAfileActivate(config->files, false, NULL);
+    ppStackFileActivation(config, PPSTACK_FILES_PREPARE, true);
+    pmFPAview *view = ppStackFilesIterateDown(config);
+    if (!view) {
+        return false;
+    }
+
+    psMetadataIterator *fileIter = psMetadataIteratorAlloc(config->files, PS_LIST_HEAD, "^PPSTACK.INPUT$");
+    psMetadataItem *fileItem; // Item from iteration
+    psArray *psfs = psArrayAlloc(num); // PSFs for PSF envelope
+    int numCols = 0, numRows = 0;   // Size of image
+    int index = 0;                  // Index for file
+    while ((fileItem = psMetadataGetAndIncrement(fileIter))) {
+        assert(fileItem->type == PS_DATA_UNKNOWN);
+        pmFPAfile *inputFile = fileItem->data.V; // An input file
+
+        // Get list of PSFs, to determine target PSF
+        if (options->convolve) {
             pmChip *chip = pmFPAviewThisChip(view, inputFile->fpa); // The chip: holds the PSF
             pmPSF *psf = psMetadataLookupPtr(NULL, chip->analysis, "PSPHOT.PSF"); // PSF
@@ -173 +173 @@
                 numRows = naxis2;
             }
-
-            pmReadout *ro = pmFPAviewThisReadout(view, inputFile->fpa); // Readout with sources
-            psArray *sources = psMetadataLookupPtr(NULL, ro->analysis, "PSPHOT.SOURCES"); // Sources
-            if (!sources) {
-                psError(PS_ERR_UNEXPECTED_NULL, false, "Unable to find sources in readout.");
-                return NULL;
-            }
-            options->sourceLists->data[index] = psMemIncrRefCounter(sources);
-
-            // Re-do photometry if we don't trust the source lists
-            if (psMetadataLookupBool(NULL, recipe, "PHOT")) {
-                psTrace("ppStack", 2, "Photometering input %d of %d....\n", index, num);
-                pmFPAfileActivate(config->files, false, NULL);
-                ppStackFileActivationSingle(config, PPSTACK_FILES_CONVOLVE, true, index);
-                pmFPAfile *file = pmFPAfileSelectSingle(config->files, "PPSTACK.INPUT", index); // File
-                pmFPAview *photView = ppStackFilesIterateDown(config);
-                if (!photView) {
-                    psFree(view);
-                    return false;
-                }
-
-                pmReadout *ro = pmFPAviewThisReadout(view, file->fpa); // Readout of interest
-
-                if (!ppStackInputPhotometer(ro, sources, config)) {
-                    psError(PS_ERR_UNKNOWN, false, "Unable to do photometry on input sources");
-                    psFree(view);
-                    psFree(photView);
-                    return false;
-                }
-
+        }
+
+
+        pmReadout *ro = pmFPAviewThisReadout(view, inputFile->fpa); // Readout with sources
+        psArray *sources = psMetadataLookupPtr(NULL, ro->analysis, "PSPHOT.SOURCES"); // Sources
+        if (!sources) {
+            psError(PS_ERR_UNEXPECTED_NULL, false, "Unable to find sources in readout.");
+            return NULL;
+        }
+        options->sourceLists->data[index] = psMemIncrRefCounter(sources);
+
+        // Re-do photometry if we don't trust the source lists
+        if (psMetadataLookupBool(NULL, recipe, "PHOT")) {
+            psTrace("ppStack", 2, "Photometering input %d of %d....\n", index, num);
+            pmFPAfileActivate(config->files, false, NULL);
+            ppStackFileActivationSingle(config, PPSTACK_FILES_CONVOLVE, true, index);
+            pmFPAfile *file = pmFPAfileSelectSingle(config->files, "PPSTACK.INPUT", index); // File
+            pmFPAview *photView = ppStackFilesIterateDown(config);
+            if (!photView) {
+                psFree(view);
+                return false;
+            }
+
+            pmReadout *ro = pmFPAviewThisReadout(view, file->fpa); // Readout of interest
+
+            if (!ppStackInputPhotometer(ro, sources, config)) {
+                psError(PS_ERR_UNKNOWN, false, "Unable to do photometry on input sources");
+                psFree(view);
                 psFree(photView);
-                if (!ppStackFilesIterateUp(config)) {
-                    psFree(view);
-                    return false;
-                }
-                pmFPAfileActivate(config->files, false, NULL);
-                ppStackFileActivation(config, PPSTACK_FILES_PREPARE, true);
-            }
-
-            index++;
-        }
-        psFree(fileIter);
-
-        // Zero point calibration
-        options->sumExposure = ppStackSourcesTransparency(options->sourceLists, options->inputMask,
-                                                          view, config);
-        if (!isfinite(options->sumExposure) || options->sumExposure <= 0) {
-            psError(PS_ERR_UNKNOWN, false, "Unable to calculate transparency differences");
-            psFree(view);
-            return false;
-        }
-
-        // Generate target PSF
+                return false;
+            }
+
+            psFree(photView);
+            if (!ppStackFilesIterateUp(config)) {
+                psFree(view);
+                return false;
+            }
+            pmFPAfileActivate(config->files, false, NULL);
+            ppStackFileActivation(config, PPSTACK_FILES_PREPARE, true);
+        }
+
+        index++;
+    }
+    psFree(fileIter);
+
+    // Generate target PSF
+    if (options->convolve) {
         options->psf = ppStackPSF(config, numCols, numRows, psfs, options->inputMask);
         psFree(psfs);
@@ -240 +234 @@
                          "Target PSF", options->psf);
         outChip->data_exists = true;
-
+    }
+
+    // Zero point calibration
+    if (!ppStackSourcesTransparency(options, view, config)) {
+        psError(PS_ERR_UNKNOWN, false, "Unable to calculate transparency differences");
         psFree(view);
-
-        if (!ppStackFilesIterateUp(config)) {
-            return false;
-        }
+        return false;
+    }
+    psFree(view);
+
+    if (!ppStackFilesIterateUp(config)) {
+        return false;
     }
 

trunk/ppStack/src/ppStackReject.c

@@ -14 +14 @@
     psAssert(options, "Require options");
     psAssert(config, "Require configuration");
+
+    if (!options->convolve) {
+        // No need to do complicated rejection when we haven't convolved
+        return true;
+    }
 
     int num = options->num;             // Number of inputs
@@ -86 +91 @@
 
         psPixels *reject = pmStackReject(options->inspect->data[i], options->numCols, options->numRows,
-                                         threshold, poorFrac, stride, options->subRegions->data[i],
-                                         options->subKernels->data[i]); // Rejected pixels
+                                         threshold, poorFrac, stride, options->regions->data[i],
+                                         options->kernels->data[i]); // Rejected pixels
 
 #ifdef TESTING
@@ -141 +146 @@
 
     psFree(options->inspect); options->inspect = NULL;
-    psFree(options->subKernels); options->subKernels = NULL;
-    psFree(options->subRegions); options->subRegions = NULL;
+    psFree(options->kernels); options->kernels = NULL;
+    psFree(options->regions); options->regions = NULL;
 
     if (numRejected >= num - 1) {

trunk/ppStack/src/ppStackSetup.c

@@ -11 +11 @@
 #include "ppStackLoop.h"
 
+#define BUFFER 16                       // Buffer for name array
+
+// Generate an array of input filenames
+static psArray *stackNameArray(pmConfig *config, // Configuration
+                               const char *name // Name of file
+    )
+{
+    psAssert(config, "Require configuration");
+    psAssert(config, "Require file name");
+
+    psString regex = NULL;             // Regular expression for selecting file
+    psStringAppend(&regex, "^%s$", name);
+    psMetadataIterator *iter = psMetadataIteratorAlloc(config->files, PS_LIST_HEAD, regex); // Iterator
+    psFree(regex);
+    psMetadataItem *item;               // Item from iteration
+    psArray *array = psArrayAlloc(BUFFER); // Array of file names
+    while ((item = psMetadataGetAndIncrement(iter))) {
+        psAssert(item->type == PS_DATA_UNKNOWN, "Should be this type");
+        pmFPAfile *file = item->data.V; // An input file
+
+        psString filename = psStringCopy(file->filename); // Filename of interest
+        psArrayAdd(array, array->n, filename);
+        psFree(filename);               // Drop reference
+    }
+    psFree(iter);
+
+    return array;
+}
+
+
 bool ppStackSetup(ppStackOptions *options, pmConfig *config)
 {
@@ -18 +48 @@
     psMetadata *recipe = psMetadataLookupMetadata(NULL, config->recipes, PPSTACK_RECIPE); // ppStack recipe
     psAssert(recipe, "We've thrown an error on this before.");
+
+    options->convolve = psMetadataLookupBool(NULL, config->arguments, "CONVOLVE"); // Convolve images?
+    if (!psMetadataLookupBool(NULL, config->arguments, "HAVE.PSF")) {
+        psWarning("No PSFs provided --- unable to convolve to common PSF.");
+        options->convolve = false;
+    }
 
     int num = psMetadataLookupS32(NULL, config->arguments, "INPUTS.NUM"); // Number of inputs
@@ -36 +72 @@
     }
 
-    const char *tempDir = psMetadataLookupStr(NULL, recipe, "TEMP.DIR"); // Directory for temporary images
-    if (!tempDir) {
-        psError(PS_ERR_BAD_PARAMETER_VALUE, false, "Unable to find TEMP.DIR in recipe");
-        return false;
+    // Generate temporary names for convolved images
+    if (options->convolve) {
+        const char *tempDir = psMetadataLookupStr(NULL, recipe, "TEMP.DIR"); // Directory for temporary images
+        if (!tempDir) {
+            psError(PS_ERR_BAD_PARAMETER_VALUE, false, "Unable to find TEMP.DIR in recipe");
+            return false;
+        }
+
+        psString outputName = psStringCopy(psMetadataLookupStr(NULL, config->arguments,
+                                                               "OUTPUT")); // Name for temporary files
+        const char *tempName = basename(outputName);
+        if (!tempName) {
+            psError(PS_ERR_BAD_PARAMETER_VALUE, false, "Unable to construct basename for temporary files.");
+            psFree(outputName);
+            return false;
+        }
+
+        const char *tempImage = psMetadataLookupStr(NULL, recipe, "TEMP.IMAGE"); // Suffix for images
+        const char *tempMask = psMetadataLookupStr(NULL, recipe, "TEMP.MASK"); // Suffix for masks
+        const char *tempVariance = psMetadataLookupStr(NULL, recipe, "TEMP.VARIANCE"); // Suffix for var maps
+        if (!tempImage || !tempMask || !tempVariance) {
+            psError(PS_ERR_BAD_PARAMETER_VALUE, false,
+                    "Unable to find TEMP.IMAGE, TEMP.MASK and TEMP.VARIANCE in recipe");
+            psFree(outputName);
+            return false;
+        }
+
+        options->imageNames = psArrayAlloc(num);
+        options->maskNames = psArrayAlloc(num);
+        options->varianceNames = psArrayAlloc(num);
+        for (int i = 0; i < num; i++) {
+            psString imageName = NULL, maskName = NULL, varianceName = NULL; // Names for convolved images
+            psStringAppend(&imageName, "%s/%s.%d.%s", tempDir, tempName, i, tempImage);
+            psStringAppend(&maskName, "%s/%s.%d.%s", tempDir, tempName, i, tempMask);
+            psStringAppend(&varianceName, "%s/%s.%d.%s", tempDir, tempName, i, tempVariance);
+            psTrace("ppStack", 5, "Temporary files: %s %s %s\n", imageName, maskName, varianceName);
+            options->imageNames->data[i] = imageName;
+            options->maskNames->data[i] = maskName;
+            options->varianceNames->data[i] = varianceName;
+        }
+        psFree(outputName);
+    } else {
+        options->imageNames = stackNameArray(config, "PPSTACK.INPUT");
+        options->maskNames = stackNameArray(config, "PPSTACK.INPUT.MASK");
+        options->varianceNames = stackNameArray(config, "PPSTACK.INPUT.VARIANCE");
     }
-
-    psString outputName = psStringCopy(psMetadataLookupStr(NULL, config->arguments,
-                                                           "OUTPUT")); // Name for temporary files
-    const char *tempName = basename(outputName);
-    if (!tempName) {
-        psError(PS_ERR_BAD_PARAMETER_VALUE, false, "Unable to construct basename for temporary files.");
-        psFree(outputName);
-        return false;
-    }
-
-    const char *tempImage = psMetadataLookupStr(NULL, recipe, "TEMP.IMAGE"); // Suffix for temporary images
-    const char *tempMask = psMetadataLookupStr(NULL, recipe, "TEMP.MASK"); // Suffix for temporary masks
-    const char *tempVariance = psMetadataLookupStr(NULL, recipe, "TEMP.VARIANCE"); // Suffix for temp var maps
-    if (!tempImage || !tempMask || !tempVariance) {
-        psError(PS_ERR_BAD_PARAMETER_VALUE, false,
-                "Unable to find TEMP.IMAGE, TEMP.MASK and TEMP.VARIANCE in recipe");
-        psFree(outputName);
-        return false;
-    }
-
-    options->imageNames = psArrayAlloc(num);
-    options->maskNames = psArrayAlloc(num);
-    options->varianceNames = psArrayAlloc(num);
-    for (int i = 0; i < num; i++) {
-        psString imageName = NULL, maskName = NULL, varianceName = NULL; // Names for convolved images
-        psStringAppend(&imageName, "%s/%s.%d.%s", tempDir, tempName, i, tempImage);
-        psStringAppend(&maskName, "%s/%s.%d.%s", tempDir, tempName, i, tempMask);
-        psStringAppend(&varianceName, "%s/%s.%d.%s", tempDir, tempName, i, tempVariance);
-        psTrace("ppStack", 5, "Temporary files: %s %s %s\n", imageName, maskName, varianceName);
-        options->imageNames->data[i] = imageName;
-        options->maskNames->data[i] = maskName;
-        options->varianceNames->data[i] = varianceName;
-    }
-    psFree(outputName);
 
     if (!pmConfigMaskSetBits(NULL, NULL, config)) {

trunk/ppStack/src/ppStackSources.c

@@ -9 +9 @@
 //#define TESTING                         // Enable debugging output
 
+#ifdef TESTING
 // Size of fake image; set by hand because it's trouble to get it from other places
 #define FAKE_COLS 4861
 #define FAKE_ROWS 4913
+#endif
 
 #ifdef TESTING
@@ -53 +55 @@
 
 
-float ppStackSourcesTransparency(const psArray *sourceLists, psVector *inputMask,
-                                 const pmFPAview *view, const pmConfig *config)
+bool ppStackSourcesTransparency(ppStackOptions *options, const pmFPAview *view, const pmConfig *config)
 {
-    PS_ASSERT_ARRAY_NON_NULL(sourceLists, NAN);
-    PS_ASSERT_VECTOR_NON_NULL(inputMask, NAN);
-    PS_ASSERT_VECTOR_TYPE(inputMask, PS_TYPE_U8, NAN);
-    PS_ASSERT_VECTOR_SIZE(inputMask, sourceLists->n, NAN);
-    PS_ASSERT_PTR_NON_NULL(view, NAN);
-    PS_ASSERT_PTR_NON_NULL(config, NAN);
+    PS_ASSERT_PTR_NON_NULL(options, false);
+    PS_ASSERT_PTR_NON_NULL(view, false);
+    PS_ASSERT_PTR_NON_NULL(config, false);
+
+    psArray *sourceLists = options->sourceLists; // Source lists for each input
+    psVector *inputMask = options->inputMask; // Mask for inputs
+
+    PS_ASSERT_ARRAY_NON_NULL(sourceLists, false);
+    PS_ASSERT_VECTOR_NON_NULL(inputMask, false);
+    PS_ASSERT_VECTOR_TYPE(inputMask, PS_TYPE_U8, false);
+    PS_ASSERT_VECTOR_SIZE(inputMask, sourceLists->n, false);
 
 #if defined(TESTING) && 0
@@ -100 +106 @@
     if (!airmassZP) {
         psError(PS_ERR_UNKNOWN, false, "Unable to find ZP.AIRMASS in recipe.");
-        return NAN;
+        return false;
     }
 
@@ -139 +145 @@
                     exptime, airmass, expFilter);
             psFree(zp);
-            return NAN;
+            return false;
         }
 
@@ -149 +155 @@
                         "Unable to find airmass term (ZP.AIRMASS) for filter %s", filter);
                 psFree(zp);
-                return NAN;
+                return false;
             }
         } else if (strcmp(filter, expFilter) != 0) {
             psError(PS_ERR_BAD_PARAMETER_VALUE, false, "Filters don't match: %s vs %s", filter, expFilter);
             psFree(zp);
-            return NAN;
+            return false;
         }
 
@@ -160 +166 @@
         sumExpTime += exptime;
     }
+
+    options->sumExposure = sumExpTime;
 
     psArray *matches = pmSourceMatchSources(sourceLists, radius, true); // List of matches
@@ -165 +173 @@
         psError(PS_ERR_UNKNOWN, false, "Unable to match sources");
         psFree(zp);
-        return NAN;
+        return false;
     }
 
@@ -177 +185 @@
     if (!trans) {
         psError(PS_ERR_UNKNOWN, false, "Unable to measure transparencies");
-        return NAN;
+        return false;
     }
 
@@ -186 +194 @@
     for (int i = 0; i < trans->n; i++) {
         if (!isfinite(trans->data.F32[i])) {
-            inputMask->data.U8[i] = PPSTACK_MASK_CAL;
+            inputMask->data.U8[i] |= PPSTACK_MASK_CAL;
         }
     }
@@ -223 +231 @@
     // m_0 = m_1 + zp_1 + trans_1 - c1 * airmass_0 - 2.5log(t_0)
     // We don't need to know the magnitude zero point for the filter, since it cancels out
+    options->norm = psVectorAlloc(num, PS_TYPE_F32);
     for (int i = 0; i < num; i++) {
         if (!isfinite(trans->data.F32[i])) {
@@ -229 +238 @@
         psArray *sources = sourceLists->data[i]; // Sources of interest
         float magCorr = airmassTerm - 2.5*log10(sumExpTime) - zp->data.F32[i] - trans->data.F32[i];
+        options->norm->data.F32[i] = magCorr;
         psLogMsg("ppStack", PS_LOG_INFO, "Applying magnitude correction to image %d: %f\n", i, magCorr);
 
@@ -248 +258 @@
             psError(PS_ERR_UNKNOWN, false, "Unable to match sources");
             psFree(zp);
-            return NAN;
+            return false;
         }
         psVector *trans = pmSourceMatchRelphot(matches, zp, iter, tol, starLimit, transIter, transRej,
@@ -259 +269 @@
 #endif
 
-    return sumExpTime;
+    return true;
 }