Changeset 26076 for trunk/ppStack

Timestamp:

Nov 9, 2009, 2:53:12 PM (17 years ago)

Author:

Paul Price

Message:

Merging branches/pap (unconvolved stacks, reworked combinePixels function into clearer roles, only throw out most variant pixel on each iteration, throw suspect pixels out as first rejection step) stack development. I think I've got everything, but not entirely sure, since I've already merged this branch once before (for dual convolution).

Location:

trunk/ppStack

Files:

• . (modified) (1 prop)
• src/ppStack.h (modified) (2 diffs)
• src/ppStackArguments.c (modified) (3 diffs)
• src/ppStackCamera.c (modified) (3 diffs)
• src/ppStackCleanup.c (modified) (1 diff)
• src/ppStackCombineFinal.c (modified) (3 diffs)
• src/ppStackCombineInitial.c (modified) (5 diffs)
• src/ppStackCombinePrepare.c (modified) (3 diffs)
• src/ppStackConvolve.c (modified) (9 diffs)
• src/ppStackFiles.c (modified) (1 diff)
• src/ppStackLoop.c (modified) (7 diffs)
• src/ppStackLoop.h (modified) (1 diff)
• src/ppStackMatch.c (modified) (8 diffs)
• src/ppStackOptions.c (modified) (4 diffs)
• src/ppStackOptions.h (modified) (2 diffs)
• src/ppStackReadout.c (modified) (13 diffs)
• src/ppStackReject.c (modified) (7 diffs)
• src/ppStackSetup.c (modified) (2 diffs)
• src/ppStackSources.c (modified) (1 diff)
• src/ppStackThread.c (modified) (9 diffs)
• src/ppStackThread.h (modified) (1 diff)

trunk/ppStack/src/ppStack.h

@@ -59 +59 @@
 // Perform stacking on a readout
 //
-// Returns an array of pixels to inspect for each input image
+// Returns two arrays: pixels to inspect for each input image, and pixels to reject for each input image.
 psArray *ppStackReadoutInitial(const pmConfig *config,   // Configuration
                                pmReadout *outRO,   // Output readout
@@ -83 +83 @@
                          const psArray *rejected, // Array with pixels rejected in each image
                          const psVector *weightings, // Weighting factors for each image
-                         const psVector *addVariance // Additional variance for rejection
+                         const psVector *addVariance, // Additional variance for rejection
+                         bool safety,                 // Enable safety switch?
+                         const psVector *norm         // Normalisations to apply
     );
 

trunk/ppStack/src/ppStackArguments.c

@@ -148 +148 @@
     psMetadataAddStr(arguments, PS_LIST_TAIL, "-stamps", 0, "Stamps file with x,y,flux per line", NULL);
     psMetadataAddStr(arguments, PS_LIST_TAIL, "-stats", 0, "Statistics file", NULL);
-    psMetadataAddS32(arguments, PS_LIST_TAIL, "-iter", 0, "Number of rejection iterations", 0);
+    psMetadataAddF32(arguments, PS_LIST_TAIL, "-combine-iter", 0, "Number of rejection iterations per input", NAN);
     psMetadataAddF32(arguments, PS_LIST_TAIL, "-combine-rej", 0,
                      "Combination rejection thresold (sigma)", NAN);
@@ -250 +250 @@
     }
 
-    VALUE_ARG_RECIPE_INT("-iter",              "ITER",            S32, 0);
+    VALUE_ARG_RECIPE_FLOAT("-combine-iter",    "COMBINE.ITER",    F32);
     VALUE_ARG_RECIPE_FLOAT("-combine-rej",     "COMBINE.REJ",     F32);
     VALUE_ARG_RECIPE_FLOAT("-combine-sys",     "COMBINE.SYS",     F32);
@@ -260 +260 @@
     VALUE_ARG_RECIPE_FLOAT("-poor-frac",       "POOR.FRACTION",   F32);
 
-    valueArgRecipeStr(arguments, recipe, "-mask-val",  "MASK.VAL",  recipe);
+    valueArgRecipeStr(arguments, recipe, "-mask-val",  "MASK.IN",   recipe);
     valueArgRecipeStr(arguments, recipe, "-mask-bad",  "MASK.BAD",  recipe);
     valueArgRecipeStr(arguments, recipe, "-mask-poor", "MASK.POOR", recipe);

trunk/ppStack/src/ppStackCamera.c

@@ -233 +233 @@
 
     // Output image
-    pmFPA *fpa = pmFPAConstruct(config->camera, config->cameraName); // FPA to contain the output
-    if (!fpa) {
+    pmFPA *outFPA = pmFPAConstruct(config->camera, config->cameraName); // FPA to contain the output
+    if (!outFPA) {
         psError(PS_ERR_UNEXPECTED_NULL, false, "Unable to construct an FPA from camera configuration.");
         return false;
     }
-    pmFPAfile *output = pmFPAfileDefineOutput(config, fpa, "PPSTACK.OUTPUT");
-    psFree(fpa);                        // Drop reference
+    pmFPAfile *output = pmFPAfileDefineOutput(config, outFPA, "PPSTACK.OUTPUT");
+    psFree(outFPA);                        // Drop reference
     if (!output) {
         psError(PS_ERR_IO, false, _("Unable to generate output file from PPSTACK.OUTPUT"));
@@ -250 +250 @@
     output->save = true;
 
-    if (!pmFPAAddSourceFromFormat(fpa, "Stack", output->format)) {
+    if (!pmFPAAddSourceFromFormat(outFPA, "Stack", output->format)) {
         psError(PS_ERR_UNKNOWN, false, "Unable to generate output FPA.");
-        psFree(fpa);
         return false;
     }
@@ -294 +293 @@
         targetPSF->save = true;
     }
+
+#if 1
+    // Unconvolved stack
+    pmFPA *unconvFPA = pmFPAConstruct(config->camera, config->cameraName); // FPA to contain unconvolved output
+    if (!unconvFPA) {
+        psError(PS_ERR_UNEXPECTED_NULL, false, "Unable to construct an FPA from camera configuration.");
+        return false;
+    }
+    pmFPAfile *unConv = pmFPAfileDefineOutput(config, unconvFPA, "PPSTACK.UNCONV");
+    psFree(unconvFPA);                  // Drop reference
+    if (!unConv) {
+        psError(PS_ERR_IO, false, _("Unable to generate output file from PPSTACK.UNCONV"));
+        return false;
+    }
+    if (unConv->type != PM_FPA_FILE_IMAGE) {
+        psError(PS_ERR_IO, true, "PPSTACK.UNCONV is not of type IMAGE");
+        return false;
+    }
+    unConv->save = true;
+
+    if (!pmFPAAddSourceFromFormat(unconvFPA, "Stack", unConv->format)) {
+        psError(PS_ERR_UNKNOWN, false, "Unable to generate output FPA.");
+        return false;
+    }
+
+    // Unconvolved mask
+    pmFPAfile *unconvMask = pmFPAfileDefineOutput(config, unconvFPA, "PPSTACK.UNCONV.MASK");
+    if (!unconvMask) {
+        psError(PS_ERR_IO, false, _("Unable to generate output file from PPSTACK.UNCONV.MASK"));
+        return false;
+    }
+    if (unconvMask->type != PM_FPA_FILE_MASK) {
+        psError(PS_ERR_IO, true, "PPSTACK.UNCONV.MASK is not of type MASK");
+        return false;
+    }
+    unconvMask->save = true;
+
+    // Unconvolved variance
+    if (haveVariances) {
+        pmFPAfile *unconvVariance = pmFPAfileDefineOutput(config, unconvFPA, "PPSTACK.UNCONV.VARIANCE");
+        if (!unconvVariance) {
+            psError(PS_ERR_IO, false, _("Unable to generate output file from PPSTACK.UNCONV.VARIANCE"));
+            return false;
+        }
+        if (unconvVariance->type != PM_FPA_FILE_VARIANCE) {
+            psError(PS_ERR_IO, true, "PPSTACK.UNCONV.VARIANCE is not of type VARIANCE");
+            return false;
+        }
+        unconvVariance->save = true;
+    }
+#endif
 
     // Output JPEGs

trunk/ppStack/src/ppStackCleanup.c

@@ -74 +74 @@
 
         if (tempDelete) {
-            psString imageResolved = pmConfigConvertFilename(options->imageNames->data[i],
+            psString imageResolved = pmConfigConvertFilename(options->convImages->data[i],
                                                              config, false, false);
-            psString maskResolved = pmConfigConvertFilename(options->maskNames->data[i],
+            psString maskResolved = pmConfigConvertFilename(options->convMasks->data[i],
                                                             config, false, false);
-            psString varianceResolved = pmConfigConvertFilename(options->varianceNames->data[i],
+            psString varianceResolved = pmConfigConvertFilename(options->convVariances->data[i],
                                                                 config, false, false);
             if (unlink(imageResolved) == -1 || unlink(maskResolved) == -1 ||

trunk/ppStack/src/ppStackCombineFinal.c

@@ -10 +10 @@
 #include "ppStackLoop.h"
 
-bool ppStackCombineFinal(ppStackThreadData *stack, ppStackOptions *options, pmConfig *config)
+//#define TESTING                         // Enable test output
+
+bool ppStackCombineFinal(pmReadout *target, ppStackThreadData *stack, psArray *covariances,
+                         ppStackOptions *options, pmConfig *config, bool safe, bool normalise)
 {
     psAssert(stack, "Require stack");
     psAssert(options, "Require options");
     psAssert(config, "Require configuration");
+
+    if (!target->mask) {
+        target->mask = psImageAlloc(target->image->numCols, target->image->numRows, PS_TYPE_IMAGE_MASK);
+    }
 
     stack->lastScan = 0;            // Reset read
@@ -30 +37 @@
         }
 
-        // call: ppStackReadoutFinal(config, outRO, readouts, rejected)
+        // call: ppStackReadoutFinal(config, target, readouts, rejected)
         psThreadJob *job = psThreadJobAlloc("PPSTACK_FINAL_COMBINE"); // Job to start
+        psArrayAdd(job->args, 1, target);
         psArrayAdd(job->args, 1, thread);
         psArrayAdd(job->args, 1, options);
         psArrayAdd(job->args, 1, config);
+        PS_ARRAY_ADD_SCALAR(job->args, safe, PS_TYPE_U8);
+        PS_ARRAY_ADD_SCALAR(job->args, normalise, PS_TYPE_U8);
         if (!psThreadJobAddPending(job)) {
             psFree(job);
@@ -48 +58 @@
 
     // Sum covariance matrices
-    double sumWeights = 0.0;            // Sum of weights
-    for (int i = 0; i < options->num; i++) {
-        if (options->inputMask->data.U8[i]) {
-            psFree(options->covariances->data[i]);
-            options->covariances->data[i] = NULL;
-            continue;
+    if (covariances) {
+        double sumWeights = 0.0;            // Sum of weights
+        for (int i = 0; i < options->num; i++) {
+            if (options->inputMask->data.U8[i]) {
+                psFree(covariances->data[i]);
+                covariances->data[i] = NULL;
+                continue;
+            }
+            psKernel *covar = covariances->data[i]; // Covariance matrix
+            if (!covar) {
+                continue;
+            }
+            float weight = options->weightings->data.F32[i]; // Weight to apply
+            psBinaryOp(covar->image, covar->image, "*", psScalarAlloc(weight, PS_TYPE_F32));
+            sumWeights += weight;
         }
-        psKernel *covar = options->covariances->data[i]; // Covariance matrix
-        if (!covar) {
-            continue;
+        if (sumWeights > 0.0) {
+            target->covariance = psImageCovarianceSum(covariances);
+            psBinaryOp(target->covariance->image, target->covariance->image, "/",
+                       psScalarAlloc(sumWeights, PS_TYPE_F32));
+            psImageCovarianceTransfer(target->variance, target->covariance);
         }
-        float weight = options->weightings->data.F32[i]; // Weight to apply
-        psBinaryOp(covar->image, covar->image, "*", psScalarAlloc(weight, PS_TYPE_F32));
-        sumWeights += weight;
-    }
-    if (sumWeights > 0.0) {
-        pmReadout *outRO = options->outRO;  // Output readout
-        outRO->covariance = psImageCovarianceSum(options->covariances);
-        psBinaryOp(outRO->covariance->image, outRO->covariance->image, "/",
-                   psScalarAlloc(sumWeights, PS_TYPE_F32));
-        psFree(options->covariances); options->covariances = NULL;
-        psImageCovarianceTransfer(outRO->variance, outRO->covariance);
+    } else {
+        target->covariance = psImageCovarianceNone();
     }
 
 #ifdef TESTING
-    pmStackVisualPlotTestImage(outRO->image, "combined_initial.fits");
-    ppStackWriteImage("combined_final.fits", NULL, outRO->image, config);
+    static int pass = 0;                // Pass through
+    psString name = NULL;               // Name of file
+    psStringAppend(&name, "combined_image_final_%d.fits", pass);
+    pass++;
+    ppStackWriteImage(name, NULL, target->image, config);
+    psStringSubstitute(&name, "mask", "image");
+    ppStackWriteImage(name, NULL, target->mask, config);
+    psStringSubstitute(&name, "variance", "mask");
+    ppStackWriteImage(name, NULL, target->variance, config);
+    psFree(name);
+
+    pmStackVisualPlotTestImage(target->image, "combined_image_final.fits");
 #endif
 

trunk/ppStack/src/ppStackCombineInitial.c

@@ -9 +9 @@
 #include "ppStack.h"
 #include "ppStackLoop.h"
+
+//#define TESTING                         // Enable test output
 
 bool ppStackCombineInitial(ppStackThreadData *stack, ppStackOptions *options, pmConfig *config)
@@ -60 +62 @@
     // Harvest the jobs, gathering the inspection lists
     options->inspect = psArrayAlloc(options->num);
+    options->rejected = psArrayAlloc(options->num);
     for (int i = 0; i < options->num; i++) {
         if (options->inputMask->data.PS_TYPE_VECTOR_MASK_DATA[i]) {
@@ -65 +68 @@
         }
         options->inspect->data[i] = psArrayAllocEmpty(numChunk);
+        options->rejected->data[i] = psArrayAllocEmpty(numChunk);
     }
     psThreadJob *job;               // Completed job
@@ -71 +75 @@
                  "Job has incorrect type: %s", job->type);
         psArray *results = job->results; // Results of job
+        psAssert(results->n == 2, "Results array has wrong size!");
+        psArray *inspect = results->data[0]; // Pixels to inspect
+        psArray *reject = results->data[1];  // Pixels to reject
         for (int i = 0; i < options->num; i++) {
             if (options->inputMask->data.PS_TYPE_VECTOR_MASK_DATA[i]) {
                 continue;
             }
-            options->inspect->data[i] = psArrayAdd(options->inspect->data[i], 1, results->data[i]);
+            options->inspect->data[i] = psArrayAdd(options->inspect->data[i], 1, inspect->data[i]);
+            options->rejected->data[i] = psArrayAdd(options->rejected->data[i], 1, reject->data[i]);
         }
         psFree(job);
@@ -83 +91 @@
 
 #ifdef TESTING
-    ppStackWriteImage("combined_initial.fits", NULL, outRO->image, config);
-    pmStackVisualPlotTestImage(outRO->image, "combined_initial.fits");
+    ppStackWriteImage("combined_image_initial.fits", NULL, options->outRO->image, config);
+    ppStackWriteImage("combined_mask_initial.fits", NULL, options->outRO->mask, config);
+    ppStackWriteImage("combined_variance_initial.fits", NULL, options->outRO->variance, config);
+
+    pmStackVisualPlotTestImage(options->outRO->image, "combined_image_initial.fits");
 #endif
-
-    psFree(options->matchChi2); options->matchChi2 = NULL;
-
 
     if (options->stats) {

trunk/ppStack/src/ppStackCombinePrepare.c

@@ -33 +33 @@
     pmCell *outCell = pmFPAfileThisCell(config->files, view, "PPSTACK.OUTPUT"); // Output cell
     options->outRO = pmReadoutAlloc(outCell); // Output readout
+
+    pmCell *unconvCell = pmFPAfileThisCell(config->files, view, "PPSTACK.UNCONV"); // Unconvolved cell
+    options->unconvRO = pmReadoutAlloc(unconvCell); // Unconvolved readout
+
     psFree(view);
 
@@ -39 +43 @@
     psString maskBadStr = psMetadataLookupStr(NULL, recipe, "MASK.BAD"); // Name of bits to mask for bad
     psImageMaskType maskBad = pmConfigMaskGet(maskBadStr, config); // Bits to mask for bad pixels
+
     if (!pmReadoutStackDefineOutput(options->outRO, col0, row0, numCols, numRows, true, true, maskBad)) {
         psError(PS_ERR_UNKNOWN, false, "Unable to prepare output.");
@@ -44 +49 @@
     }
 
+    options->unconvRO->image = psImageCopy(NULL, options->outRO->image, PS_TYPE_F32);
+//    options->unconvRO->mask = psImageCopy(NULL, options->outRO->mask, PS_TYPE_IMAGE_MASK);
+    options->unconvRO->col0 = options->outRO->col0;
+    options->unconvRO->row0 = options->outRO->row0;
+
+
+
     return true;
 }

trunk/ppStack/src/ppStackConvolve.c

@@ -9 +9 @@
 #include "ppStack.h"
 #include "ppStackLoop.h"
+
+//#define TESTING
+
 
 // Update the value of a concept
@@ -39 +42 @@
     options->weightings = psVectorAlloc(num, PS_TYPE_F32); // Combination weightings for images (1/noise^2)
     psVectorInit(options->weightings, 0.0);
-    options->covariances = psArrayAlloc(num); // Covariance matrices
+    options->origCovars = psArrayAlloc(num);
+    options->convCovars = psArrayAlloc(num); // Covariance matrices
+
+    psVector *renorms = psVectorAlloc(num, PS_TYPE_F32); // Renormalisation values for variances
+    psVectorInit(renorms, NAN);
 
     psList *fpaList = psListAlloc(NULL); // List of input FPAs, for concept averaging
@@ -79 +86 @@
         // Background subtraction, scaling and normalisation is performed automatically by the image matching
         psTimerStart("PPSTACK_MATCH");
+        options->origCovars->data[i] = psMemIncrRefCounter(readout->covariance);
         if (!ppStackMatch(readout, options, i, config)) {
             psErrorStackPrint(stderr, "Unable to match image %d --- ignoring.", i);
@@ -85 +93 @@
             continue;
         }
-        options->covariances->data[i] = psMemIncrRefCounter(readout->covariance);
+        options->convCovars->data[i] = psMemIncrRefCounter(readout->covariance);
+
+        float renorm = psMetadataLookupF32(NULL, readout->analysis, PM_READOUT_ANALYSIS_RENORM);
+        if (!isfinite(renorm)) {
+            renorm = 1.0;
+        }
+        renorms->data.F32[i] = renorm;
 
         if (options->stats) {
@@ -114 +128 @@
         pmHDU *hdu = readout->parent->parent->parent->hdu; // HDU for convolved image
         assert(hdu);
-        ppStackWriteImage(options->imageNames->data[i], hdu->header, readout->image, config);
+        ppStackWriteImage(options->convImages->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);
+        ppStackWriteImage(options->convMasks->data[i], maskHeader, readout->mask, config);
         psFree(maskHeader);
-        psImageCovarianceTransfer(readout->variance, readout->covariance);
-        ppStackWriteImage(options->varianceNames->data[i], hdu->header, readout->variance, config);
+        ppStackWriteImage(options->convVariances->data[i], hdu->header, readout->variance, config);
 #ifdef TESTING
         {
@@ -129 +142 @@
             psFree(name);
         }
+        {
+            int numCols = readout->image->numCols, numRows = readout->image->numRows;
+            psImage *sn = psImageAlloc(numCols, numRows, PS_TYPE_F32);
+            for (int y = 0; y < numRows; y++) {
+                for (int x = 0; x < numCols; x++) {
+                    sn->data.F32[y][x] = readout->image->data.F32[y][x] /
+                        sqrtf(readout->variance->data.F32[y][x]);
+                }
+            }
+            psString name = NULL;
+            psStringAppend(&name, "signoise_%d.fits", i);
+            ppStackWriteImage(name, hdu->header, sn, config);
+            psFree(name);
+            psFree(sn);
+        }
 #endif
 
@@ -149 +177 @@
     psFree(rng);
 
-    psFree(options->norm); options->norm = NULL;
     psFree(options->sourceLists); options->sourceLists = NULL;
     psFree(options->psf); options->psf = NULL;
@@ -211 +238 @@
             numGood = 0;                    // Number of good images
             for (int i = 0; i < num; i++) {
-              if (options->inputMask->data.PS_TYPE_VECTOR_MASK_DATA[i] & PPSTACK_MASK_ALL) {
+                if (options->inputMask->data.PS_TYPE_VECTOR_MASK_DATA[i] & PPSTACK_MASK_ALL) {
                     continue;
                 }
@@ -233 +260 @@
     }
 
+    // Correct chi^2 for renormalisation
+    psBinaryOp(options->matchChi2, options->matchChi2, "/", renorms);
+    for (int i = 0; i < num; i++) {
+        psLogMsg("ppStack", PS_LOG_INFO, "Additional variance for image %d: %f\n",
+                 i, options->matchChi2->data.F32[i]);
+    }
+    psFree(renorms);
+
     return true;
 }

trunk/ppStack/src/ppStackFiles.c

@@ -22 +22 @@
 /// Output files for the combination
 static char *filesCombine[] = { "PPSTACK.OUTPUT", "PPSTACK.OUTPUT.MASK", "PPSTACK.OUTPUT.VARIANCE",
+                                "PPSTACK.UNCONV", "PPSTACK.UNCONV.MASK", "PPSTACK.UNCONV.VARIANCE",
                                 "PPSTACK.OUTPUT.JPEG1", "PPSTACK.OUTPUT.JPEG2", NULL };
 

trunk/ppStack/src/ppStackLoop.c

@@ -56 +56 @@
     // Start threading
     ppStackThreadInit();
-    ppStackThreadData *stack = ppStackThreadDataSetup(options, config);
+    ppStackThreadData *stack = ppStackThreadDataSetup(options, config, true);
     if (!stack) {
         psError(PS_ERR_IO, false, "Unable to initialise stack threads.");
@@ -62 +62 @@
         return false;
     }
-    psFree(options->cells); options->cells = NULL;
 
     // Prepare for combination
@@ -98 +97 @@
     // Final combination
     psTrace("ppStack", 2, "Final stack of convolved images....\n");
-    if (!ppStackCombineFinal(stack, options, config)) {
+    if (!ppStackCombineFinal(options->outRO, stack, options->convCovars, options, config, false, false)) {
         psError(PS_ERR_UNKNOWN, false, "Unable to perform final combination.");
         psFree(stack);
@@ -106 +105 @@
     psLogMsg("ppStack", PS_LOG_INFO, "Stage 5: Final Stack: %f sec", psTimerClear("PPSTACK_STEPS"));
     ppStackMemDump("final");
-
 
     // Clean up
@@ -121 +119 @@
     psFree(stack);
 
+#if 1
+    // Unconvolved stack --- it's cheap to calculate, compared to everything else!
+    if (options->convolve) {
+        // Start threading
+        ppStackThreadData *stack = ppStackThreadDataSetup(options, config, false);
+        if (!stack) {
+            psError(PS_ERR_IO, false, "Unable to initialise stack threads.");
+            psFree(options);
+            return false;
+        }
+        psTrace("ppStack", 2, "Stack of unconvolved images....\n");
+        if (!ppStackCombineFinal(options->unconvRO, stack, options->origCovars, options, config, false, true)) {
+            psError(PS_ERR_UNKNOWN, false, "Unable to perform unconvolved combination.");
+            psFree(stack);
+            psFree(options);
+            return false;
+        }
+        psLogMsg("ppStack", PS_LOG_INFO, "Stage 7: Unconvolved Stack: %f sec", psTimerClear("PPSTACK_STEPS"));
+        ppStackMemDump("unconv");
+
+        psFree(stack);
+    }
+    psFree(options->cells); options->cells = NULL;
+#endif
 
     // Photometry
@@ -129 +151 @@
         return false;
     }
-    psLogMsg("ppStack", PS_LOG_INFO, "Stage 7: Photometry Analysis: %f sec", psTimerClear("PPSTACK_STEPS"));
+    psLogMsg("ppStack", PS_LOG_INFO, "Stage 8: Photometry Analysis: %f sec", psTimerClear("PPSTACK_STEPS"));
     ppStackMemDump("photometry");
-
 
     // Finish up
@@ -140 +161 @@
         return false;
     }
-    psLogMsg("ppStack", PS_LOG_INFO, "Stage 8: Final output: %f sec", psTimerClear("PPSTACK_STEPS"));
+    psLogMsg("ppStack", PS_LOG_INFO, "Stage 9: Final output: %f sec", psTimerClear("PPSTACK_STEPS"));
     ppStackMemDump("finish");
 
     psFree(options);
+
     return true;
 }

trunk/ppStack/src/ppStackLoop.h

@@ -56 +56 @@
 // Final combination
 bool ppStackCombineFinal(
+    pmReadout *target,                  // Target readout
     ppStackThreadData *stack,           // Stack
+    psArray *covariances,               // Covariances
     ppStackOptions *options,            // Options
-    pmConfig *config                    // Configuration
+    pmConfig *config,                   // Configuration
+    bool safe,                          // Allow safe combination?
+    bool norm                           // Normalise images?
     );
 

trunk/ppStack/src/ppStackMatch.c

@@ -167 +167 @@
     )
 {
+#if 1
     bool mdok; // Status of metadata lookups
 
@@ -192 +193 @@
     psImageMaskType maskBad = pmConfigMaskGet("BLANK", config); // Bits to mask
 
+    psImageCovarianceTransfer(readout->variance, readout->covariance);
     return pmReadoutVarianceRenormalise(readout, maskBad, num, minValid, maxValid);
+#else
+    return true;
+#endif
 }
 
@@ -212 +217 @@
     int size = psMetadataLookupS32(NULL, ppsub, "KERNEL.SIZE"); // Kernel half-size
 
-    psString maskValStr = psMetadataLookupStr(NULL, recipe, "MASK.VAL"); // Name of bits to mask going in
+    psString maskValStr = psMetadataLookupStr(NULL, ppsub, "MASK.VAL"); // Name of bits to mask going in
     psImageMaskType maskVal = pmConfigMaskGet(maskValStr, config); // Bits to mask going in to pmSubtractionMatch
     psString maskPoorStr = psMetadataLookupStr(NULL, recipe, "MASK.POOR"); // Name of bits to mask for poor
@@ -257 +262 @@
             psFitsClose(fits);
 
-            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)) {
+            if (!readImage(&readout->image, options->convImages->data[index], config) ||
+                !readImage(&readout->mask, options->convMasks->data[index], config) ||
+                !readImage(&readout->variance, options->convVariances->data[index], config)) {
                 psError(PS_ERR_IO, false, "Unable to read previously produced image.");
                 return false;
@@ -378 +383 @@
             }
 #endif
+
+            fprintf(stderr, "vf = %f\n", psImageCovarianceFactor(readout->covariance));
+
 
             if (threads > 0) {
@@ -517 +525 @@
             psFree(iter);
             options->matchChi2->data.F32[index] = sum / (psImageCovarianceFactor(readout->covariance) * num);
+            fprintf(stderr, "chi2 = %f ; vf = %f\n", sum/num, psImageCovarianceFactor(readout->covariance));
         }
 
@@ -543 +552 @@
     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.");
-      psErrorClear();
+        psWarning("Can't measure background for image.");
+        psErrorClear();
     } else {
-      if (!psMetadataLookupBool(NULL, config->arguments, "PPSTACK.SKIP.BG.SUB")) {
-        psLogMsg("ppStack", PS_LOG_INFO, "Correcting convolved image background by %lf (+/- %lf)",
-                 psStatsGetValue(bg, PS_STAT_ROBUST_MEDIAN), psStatsGetValue(bg, PS_STAT_ROBUST_STDEV));
-        (void)psBinaryOp(readout->image, readout->image, "-",
-                         psScalarAlloc(psStatsGetValue(bg, PS_STAT_ROBUST_MEDIAN), PS_TYPE_F32));
-      }
+        if (!psMetadataLookupBool(NULL, config->arguments, "PPSTACK.SKIP.BG.SUB")) {
+            psLogMsg("ppStack", PS_LOG_INFO, "Correcting convolved image background by %lf (+/- %lf)",
+                     psStatsGetValue(bg, PS_STAT_ROBUST_MEDIAN), psStatsGetValue(bg, PS_STAT_ROBUST_STDEV));
+            (void)psBinaryOp(readout->image, readout->image, "-",
+                             psScalarAlloc(psStatsGetValue(bg, PS_STAT_ROBUST_MEDIAN), PS_TYPE_F32));
+        }
     }
 
@@ -569 +578 @@
     options->weightings->data.F32[index] = 1.0 / (psStatsGetValue(bg, PS_STAT_ROBUST_MEDIAN) *
                                                   psImageCovarianceFactor(readout->covariance));
+    psLogMsg("ppStack", PS_LOG_INFO, "Weighting for image %d is %f\n",
+             index, options->weightings->data.F32[index]);
     psMetadataAddF32(readout->analysis, PS_LIST_TAIL, "PPSTACK.WEIGHTING", 0,
                      "Weighting by 1/noise^2 for stack", options->weightings->data.F32[index]);

trunk/ppStack/src/ppStackOptions.c

@@ -12 +12 @@
         fclose(options->statsFile);
     }
-    psFree(options->imageNames);
-    psFree(options->maskNames);
-    psFree(options->varianceNames);
+    psFree(options->origImages);
+    psFree(options->origMasks);
+    psFree(options->origVariances);
+    psFree(options->origCovars);
+    psFree(options->convImages);
+    psFree(options->convMasks);
+    psFree(options->convVariances);
     psFree(options->psf);
     psFree(options->inputMask);
@@ -24 +28 @@
     psFree(options->matchChi2);
     psFree(options->weightings);
-    psFree(options->covariances);
+    psFree(options->convCovars);
     psFree(options->outRO);
+    psFree(options->unconvRO);
     psFree(options->inspect);
     psFree(options->rejected);
@@ -39 +44 @@
     options->stats = NULL;
     options->statsFile = NULL;
-    options->imageNames = NULL;
-    options->maskNames = NULL;
-    options->varianceNames = NULL;
+    options->origImages = NULL;
+    options->origMasks = NULL;
+    options->origVariances = NULL;
+    options->origCovars = NULL;
+    options->convImages = NULL;
+    options->convMasks = NULL;
+    options->convVariances = NULL;
     options->num = 0;
     options->psf = NULL;
@@ -54 +63 @@
     options->matchChi2 = NULL;
     options->weightings = NULL;
-    options->covariances = NULL;
+    options->convCovars = NULL;
     options->outRO = NULL;
+    options->unconvRO = NULL;
     options->inspect = NULL;
     options->rejected = NULL;

trunk/ppStack/src/ppStackOptions.h

@@ -11 +11 @@
     psMetadata *stats;                  // Statistics for output
     FILE *statsFile;                    // File to which to write statistics
-    psArray *imageNames, *maskNames, *varianceNames; // Filenames for the temporary convolved images
+    psArray *origImages, *origMasks, *origVariances; // Filenames of the original images
+    psArray *convImages, *convMasks, *convVariances; // Filenames for the temporary convolved images
+    psArray *origCovars;                // Original covariances matrices
     int num;                            // Number of inputs
     // Prepare
@@ -26 +28 @@
     psVector *matchChi2;                // chi^2 for stamps from matching
     psVector *weightings;               // Combination weightings for images (1/noise^2)
-    psArray *covariances;               // Covariance matrices
+    psArray *convCovars;                // Convolved covariance matrices
     // Combine initial
     pmReadout *outRO;                   // Output readout
+    pmReadout *unconvRO;                // Unconvolved readout
     psArray *inspect;                   // Array of arrays of pixels to inspect
     // Rejection

trunk/ppStack/src/ppStackReadout.c

@@ -25 +25 @@
     psVector *addVariance = options->matchChi2; // Additional variance when rejecting
 
-    psArray *inspect = ppStackReadoutInitial(config, outRO, thread->readouts, mask,
-                                             weightings, addVariance);
-
-    job->results = inspect;
+    job->results = ppStackReadoutInitial(config, outRO, thread->readouts, mask,
+                                         weightings, addVariance);
     thread->busy = false;
 
-    return inspect ? true : false;
+    return job->results ? true : false;
 }
 
@@ -39 +37 @@
 
     psArray *args = job->args;          // Arguments
-    ppStackThread *thread = args->data[0]; // Thread
-    ppStackOptions *options = args->data[1]; // Options
-    pmConfig *config = args->data[2];   // Configuration
-
-    pmReadout *outRO = options->outRO;  // Output readout
+    pmReadout *target = args->data[0];  // Output readout
+    ppStackThread *thread = args->data[1]; // Thread
+    ppStackOptions *options = args->data[2]; // Options
+    pmConfig *config = args->data[3];   // Configuration
+    bool safety = PS_SCALAR_VALUE(args->data[4], U8);    // Safety switch on?
+    bool normalise = PS_SCALAR_VALUE(args->data[5], U8); // Normalise images?
+
     psVector *mask = options->inputMask; // Mask for inputs
     psArray *rejected = options->rejected; // Rejected pixels
     psVector *weightings = options->weightings; // Weightings (1/noise^2) for each image
     psVector *addVariance = options->matchChi2; // Additional variance when rejecting
-
-    bool status = ppStackReadoutFinal(config, outRO, thread->readouts, mask, rejected,
-                                      weightings, addVariance); // Status of operation
+    psVector *norm = normalise ? options->norm : NULL; // Normalisations to apply to images
+
+    bool status = ppStackReadoutFinal(config, target, thread->readouts, mask, rejected,
+                                      weightings, addVariance, safety, norm); // Status of operation
 
     thread->busy = false;
@@ -63 +64 @@
 
     psArray *args = job->args;  // Input arguments
-    psArray *inspect = args->data[0]; // Array of pixel arrays
-    int index = PS_SCALAR_VALUE(args->data[1], S32); // Index of interest
-
-    psArray *inputs = inspect->data[index]; // Array of interest
-    psPixels *output = NULL;    // Output pixel list
-    for (int i = 0; i < inputs->n; i++) {
-        psPixels *input = inputs->data[i]; // Input pixel list
-        if (!input || input->n == 0) {
-            continue;
-        }
-        output = psPixelsConcatenate(output, input);
-    }
-
-    if (!output) {
-        // If there are no pixels to inspect, then just fake it
-        output = psPixelsAllocEmpty(0);
-    }
-
-    psFree(inputs);
-    inspect->data[index] = output;
+    psArray *inspects = args->data[0]; // Array of pixel arrays
+    psArray *rejects = args->data[1];  // Array of pixel arrays
+    int index = PS_SCALAR_VALUE(args->data[2], S32); // Index of interest
+
+    psArray *inInspects = inspects->data[index]; // Array of interest
+    psArray *inRejects = rejects->data[index]; // Array of interest
+    psAssert(inInspects->n == inRejects->n, "Size should be the same");
+    psPixels *outInspect = NULL, *outReject = NULL; // Output pixel lists
+    for (int i = 0; i < inInspects->n; i++) {
+        psPixels *inInspect = inInspects->data[i]; // Input pixel list
+        if (inInspect && inInspect->n > 0) {
+            outInspect = psPixelsConcatenate(outInspect, inInspect);
+        }
+        psPixels *inReject = inRejects->data[i]; // Input pixel list
+        if (inReject && inReject->n > 0) {
+            outReject = psPixelsConcatenate(outReject, inReject);
+        }
+    }
+
+    // If there are no pixels to inspect, then just fake it
+    if (!outInspect) {
+        outInspect = psPixelsAllocEmpty(0);
+    }
+    if (!outReject) {
+        outReject = psPixelsAllocEmpty(0);
+    }
+
+    psFree(inspects->data[index]);
+    inspects->data[index] = outInspect;
+    psFree(rejects->data[index]);
+    rejects->data[index] = outReject;
 
     return true;
@@ -104 +116 @@
 
     bool mdok;                          // Status of MD lookup
-    int iter = psMetadataLookupS32(NULL, recipe, "ITER"); // Rejection iterations
+    float iter = psMetadataLookupF32(NULL, recipe, "COMBINE.ITER"); // Rejection iterations
     float combineRej = psMetadataLookupF32(NULL, recipe, "COMBINE.REJ"); // Combination threshold
     float combineSys = psMetadataLookupF32(NULL, recipe, "COMBINE.SYS"); // Combination systematic error
@@ -114 +126 @@
     int kernelSize = psMetadataLookupS32(NULL, ppsub, "KERNEL.SIZE"); // Kernel half-size
 
-    psString maskValStr = psMetadataLookupStr(NULL, recipe, "MASK.VAL"); // Name of bits to mask going in
+    psString maskValStr = psMetadataLookupStr(NULL, recipe, "MASK.IN"); // Name of bits to mask going in
     psImageMaskType maskVal = pmConfigMaskGet(maskValStr, config); // Bits to mask going in to pmSubtractionMatch
+    psString maskSuspectStr = psMetadataLookupStr(NULL, recipe, "MASK.SUSPECT"); // Name of suspect mask bits
+    psImageMaskType maskSuspect = pmConfigMaskGet(maskSuspectStr, config); // Suspect bits
     psString maskBadStr = psMetadataLookupStr(NULL, recipe, "MASK.BAD"); // Name of bits to mask for bad
     psImageMaskType maskBad = pmConfigMaskGet(maskBadStr, config); // Bits to mask for bad pixels
@@ -149 +163 @@
     }
 
-    if (!pmStackCombine(outRO, stack, maskVal | maskBad, maskBad, kernelSize, iter,
-                        combineRej, combineSys, combineDiscard, true, useVariance, safe, false)) {
+    if (!pmStackCombine(outRO, stack, maskVal | maskBad, maskSuspect, maskBad, kernelSize, iter,
+                        combineRej, combineSys, combineDiscard, useVariance, safe, false)) {
         psError(PS_ERR_UNKNOWN, false, "Unable to combine input readouts with rejection.");
         psFree(stack);
@@ -156 +170 @@
     }
 
-    // Save list of pixels to inspect
+    // Save lists of pixels
     psArray *inspect = psArrayAlloc(num); // List of pixels to inspect
+    psArray *reject = psArrayAlloc(num);  // List of pixels rejected
     for (int i = 0; i < num; i++) {
         pmStackData *data = stack->data[i]; // Data for this image
@@ -168 +183 @@
         }
         inspect->data[i] = psMemIncrRefCounter(data->inspect);
+        reject->data[i] = psMemIncrRefCounter(data->reject);
     }
     psFree(stack);
@@ -173 +189 @@
     sectionNum++;
 
-    return inspect;
+    psArray *results = psArrayAlloc(2); // Array of results
+    results->data[0] = inspect;
+    results->data[1] = reject;
+
+    return results;
 }
 
@@ -180 +200 @@
 bool ppStackReadoutFinal(const pmConfig *config, pmReadout *outRO, const psArray *readouts,
                          const psVector *mask, const psArray *rejected, const psVector *weightings,
-                         const psVector *addVariance)
+                         const psVector *addVariance, bool safety, const psVector *norm)
 {
     assert(config);
@@ -196 +216 @@
 
     bool mdok;                          // Status of MD lookup
-    int iter = psMetadataLookupS32(NULL, recipe, "ITER"); // Rejection iterations
-    float combineRej = psMetadataLookupF32(NULL, recipe, "COMBINE.REJ"); // Combination threshold
-    float combineSys = psMetadataLookupF32(NULL, recipe, "COMBINE.SYS"); // Combination systematic error
-    float combineDiscard = psMetadataLookupF32(NULL, recipe, "COMBINE.DISCARD"); // Olympic discard fraction
     bool useVariance = psMetadataLookupBool(&mdok, recipe, "VARIANCE"); // Use variance for rejection?
     bool safe = psMetadataLookupBool(&mdok, recipe, "SAFE"); // Be safe when combining small numbers of pixels
 
-    psString maskValStr = psMetadataLookupStr(NULL, recipe, "MASK.VAL"); // Name of bits to mask going in
+    psString maskValStr = psMetadataLookupStr(NULL, recipe, "MASK.IN"); // Name of bits to mask going in
     psImageMaskType maskVal = pmConfigMaskGet(maskValStr, config); // Bits to mask going in to pmSubtractionMatch
+    psString maskSuspectStr = psMetadataLookupStr(NULL, recipe, "MASK.SUSPECT"); // Name of suspect mask bits
+    psImageMaskType maskSuspect = pmConfigMaskGet(maskSuspectStr, config); // Suspect bits
     psString maskBadStr = psMetadataLookupStr(NULL, recipe, "MASK.BAD"); // Name of bits to mask for bad
     psImageMaskType maskBad = pmConfigMaskGet(maskBadStr, config); // Bits to mask for bad pixels
@@ -211 +229 @@
     psArray *stack = psArrayAlloc(num); // Array for stacking
 
-    bool entire = true;                 // Combine entire image?
-    if (rejected) {
-        // We have rejection from a previous combination: combine without flagging pixels to inspect
-        entire = false;
-        safe = false;
-        iter = 0;
-        combineRej = NAN;
-        combineSys = NAN;
-    }
+    // We have rejection from a previous combination: combine without flagging pixels to inspect
+    safe &= safety;
+    int iter = 0;
+    float combineRej = NAN;
+    float combineSys = NAN;
+    float combineDiscard = NAN;
 
     for (int i = 0; i < num; i++) {
         pmReadout *ro = readouts->data[i];
-        if (mask->data.U8[i] & (PPSTACK_MASK_REJECT | PPSTACK_MASK_BAD)) {
-            // Image completely rejected since previous combination
-            entire = true;
-            continue;
-        } else if (mask->data.U8[i]) {
-            // Image completely rejected before original combination
+        if (mask->data.U8[i]) {
+            // Image completely rejected
             continue;
         }
@@ -252 +263 @@
         data->reject = rejected ? psMemIncrRefCounter(rejected->data[i]) : NULL;
         stack->data[i] = data;
-    }
-
-    if (!pmStackCombine(outRO, stack, maskVal | maskBad, maskBad, 0,
-                        iter, combineRej, combineSys, combineDiscard,
-                        entire, useVariance, safe, !rejected)) {
+
+        if (norm) {
+            float normalise = powf(10.0, -0.4 * norm->data.F32[i]); // Normalisation
+            psBinaryOp(ro->image, ro->image, "*", psScalarAlloc(normalise, PS_TYPE_F32));
+            psBinaryOp(ro->variance, ro->variance, "*", psScalarAlloc(PS_SQR(normalise), PS_TYPE_F32));
+        }
+    }
+
+    if (!pmStackCombine(outRO, stack, maskVal | maskBad, maskSuspect, maskBad, 0, iter, combineRej,
+                        combineSys, combineDiscard, useVariance, safe, true)) {
         psError(PS_ERR_UNKNOWN, false, "Unable to combine input readouts.");
         psFree(stack);

trunk/ppStack/src/ppStackReject.c

@@ -23 +23 @@
 
     int num = options->num;             // Number of inputs
-    options->rejected = psArrayAlloc(num);
 
     psMetadata *recipe = psMetadataLookupMetadata(NULL, config->recipes, PPSTACK_RECIPE); // ppStack recipe
@@ -53 +52 @@
         psThreadJob *job = psThreadJobAlloc("PPSTACK_INSPECT"); // Job to start
         psArrayAdd(job->args, 1, options->inspect);
+        psArrayAdd(job->args, 1, options->rejected);
         PS_ARRAY_ADD_SCALAR(job->args, i, PS_TYPE_S32);
         if (!psThreadJobAddPending(job)) {
@@ -92 +92 @@
 #endif
 
-        psPixels *reject = pmStackReject(options->inspect->data[i], options->numCols, options->numRows,
-                                         threshold, poorFrac, stride, options->regions->data[i],
-                                         options->kernels->data[i]); // Rejected pixels
-
 #ifdef TESTING
         {
-            psImage *mask = psPixelsToMask(NULL, reject,
+            psImage *mask = psPixelsToMask(NULL, options->rejected->data[i],
                                            psRegionSet(0, options->numCols - 1, 0, options->numRows - 1),
                                            0xff); // Mask image
             psString name = NULL;           // Name of image
-            psStringAppend(&name, "reject_%03d.fits", i);
+            psStringAppend(&name, "pre_reject_%03d.fits", i);
             pmStackVisualPlotTestImage(mask, name);
             psFits *fits = psFitsOpen(name, "w");
@@ -111 +107 @@
         }
 #endif
+
+        psPixels *reject = pmStackReject(options->inspect->data[i], options->numCols, options->numRows,
+                                         threshold, poorFrac, stride, options->regions->data[i],
+                                         options->kernels->data[i]); // Rejected pixels
 
         psFree(options->inspect->data[i]);
@@ -127 +127 @@
                           "exceeds limit (%.3f)", i, frac, imageRej);
                 psFree(reject);
-                // reject == NULL means reject image completely
                 reject = NULL;
                 options->inputMask->data.PS_TYPE_VECTOR_MASK_DATA[i] |= PPSTACK_MASK_BAD;
@@ -134 +133 @@
         }
 
-        // Images without a list of rejected pixels (the list may be empty) are rejected completely
-        options->rejected->data[i] = reject;
+        if (reject) {
+            // Add to list of pixels already rejected
+            reject = psPixelsConcatenate(reject, options->rejected->data[i]);
+            options->rejected->data[i] = psPixelsDuplicates(options->rejected->data[i], reject);
+        }
+
+#ifdef TESTING
+        {
+            psImage *mask = psPixelsToMask(NULL, options->rejected->data[i],
+                                           psRegionSet(0, options->numCols - 1, 0, options->numRows - 1),
+                                           0xff); // Mask image
+            psString name = NULL;           // Name of image
+            psStringAppend(&name, "reject_%03d.fits", i);
+            pmStackVisualPlotTestImage(mask, name);
+            psFits *fits = psFitsOpen(name, "w");
+            psFree(name);
+            psFitsWriteImage(fits, NULL, mask, 0, NULL);
+            psFree(mask);
+            psFitsClose(fits);
+        }
+#endif
 
         if (options->stats) {
@@ -143 +161 @@
                              "Number of pixels rejected", reject ? reject->n : 0);
         }
+
+        psFree(reject);
         psLogMsg("ppStack", PS_LOG_INFO, "Time to perform rejection on image %d: %f sec", i,
                  psTimerClear("PPSTACK_REJECT"));

trunk/ppStack/src/ppStackSetup.c

@@ -72 +72 @@
     }
 
-    options->imageNames = psArrayAlloc(num);
-    options->maskNames = psArrayAlloc(num);
-    options->varianceNames = psArrayAlloc(num);
+    options->convImages = psArrayAlloc(num);
+    options->convMasks = psArrayAlloc(num);
+    options->convVariances = psArrayAlloc(num);
     for (int i = 0; i < num; i++) {
         psString imageName = NULL, maskName = NULL, varianceName = NULL; // Names for convolved images
@@ -81 +81 @@
         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;
+        options->convImages->data[i] = imageName;
+        options->convMasks->data[i] = maskName;
+        options->convVariances->data[i] = varianceName;
     }
     psFree(outputName);
+
+    // Original images
+    options->origImages = psArrayAlloc(num);
+    options->origMasks = psArrayAlloc(num);
+    options->origVariances = psArrayAlloc(num);
+    pmFPAview *view = pmFPAviewAlloc(0);
+    for (int i = 0; i < num; i++) {
+        {
+            pmFPAfile *file = pmFPAfileSelectSingle(config->files, "PPSTACK.INPUT", i);
+            options->origImages->data[i] = pmFPAfileName(file, view, config);
+        }
+        {
+            // We want the convolved mask, since that defines the area that has been tested for outliers
+            options->origMasks->data[i] = psMemIncrRefCounter(options->convMasks->data[i]);
+        }
+        {
+            pmFPAfile *file = pmFPAfileSelectSingle(config->files, "PPSTACK.INPUT.VARIANCE", i);
+            options->origVariances->data[i] = pmFPAfileName(file, view, config);
+        }
+    }
+    psFree(view);
 
     if (!pmConfigMaskSetBits(NULL, NULL, config)) {

trunk/ppStack/src/ppStackSources.c

@@ -188 +188 @@
     }
 
+#if 0
+    // Position offsets
+    {
+        psArray *offsets = pmSourceMatchRelastro(matches, num, tol, iter1, starRej1,
+                                                  iter2, starRej2, starLimit); // Shifts for each image
+        if (!offsets) {
+            psError(PS_ERR_UNKNOWN, false, "Unable to measure offsets");
+            return false;
+        }
+        for (int i = 0; i < num; i++) {
+            psVector *offset = offsets->data[i];
+            fprintf(stderr, "Offset %d: %f,%f\n", i, offset->data.F32[0], offset->data.F32[1]);
+        }
+        psFree(offsets);
+    }
+#endif
+
 #ifdef TESTING
     dumpMatches("source_mags.dat", num, matches, zp, trans);

trunk/ppStack/src/ppStackThread.c

@@ -56 +56 @@
 }
 
-ppStackThreadData *ppStackThreadDataSetup(const ppStackOptions *options, const pmConfig *config)
+ppStackThreadData *ppStackThreadDataSetup(const ppStackOptions *options, const pmConfig *config, bool conv)
 {
     psAssert(options, "Require options");
@@ -62 +62 @@
 
     const psArray *cells = options->cells; // Array of input cells
-    const psArray *imageNames = options->imageNames; // Names of images to read
-    const psArray *maskNames = options->maskNames; // Names of masks to read
-    const psArray *varianceNames = options->varianceNames; // Names of variance maps to read
-    const psArray *covariances = options->covariances; // Covariance matrices (already read)
+    const psArray *imageNames = conv ? options->convImages : options->origImages; // Names of images to read
+    const psArray *maskNames = conv ? options->convMasks : options->origMasks; // Names of masks to read
+    const psArray *varianceNames = conv ? options->convVariances : options->origVariances; // Variance names
+    const psArray *covariances = conv ? options->convCovars : options->origCovars; // Covariance matrices
 
     PS_ASSERT_ARRAY_NON_NULL(cells, NULL);
-    PS_ASSERT_ARRAYS_SIZE_EQUAL(cells, imageNames, NULL);
-    PS_ASSERT_ARRAYS_SIZE_EQUAL(cells, maskNames, NULL);
-    PS_ASSERT_ARRAYS_SIZE_EQUAL(cells, varianceNames, NULL);
-    PS_ASSERT_ARRAYS_SIZE_EQUAL(cells, covariances, NULL);
+    if (imageNames) {
+        PS_ASSERT_ARRAYS_SIZE_EQUAL(cells, imageNames, NULL);
+    }
+    if (maskNames) {
+        PS_ASSERT_ARRAYS_SIZE_EQUAL(cells, maskNames, NULL);
+    }
+    if (varianceNames) {
+        PS_ASSERT_ARRAYS_SIZE_EQUAL(cells, varianceNames, NULL);
+    }
+    if (covariances) {
+        PS_ASSERT_ARRAYS_SIZE_EQUAL(cells, covariances, NULL);
+    }
 
     ppStackThreadData *stack = psAlloc(sizeof(ppStackThreadData)); // Thread data, to return
@@ -87 +95 @@
             continue;
         }
-        // Resolved names
-        psString imageResolved = pmConfigConvertFilename(imageNames->data[i], config, false, false);
-        psString maskResolved = pmConfigConvertFilename(maskNames->data[i], config, false, false);
-        psString varianceResolved = pmConfigConvertFilename(varianceNames->data[i], config, false, false);
-        stack->imageFits->data[i] = psFitsOpen(imageResolved, "r");
-        stack->maskFits->data[i] = psFitsOpen(maskResolved, "r");
-        stack->varianceFits->data[i] = psFitsOpen(varianceResolved, "r");
-        psFree(imageResolved);
-        psFree(maskResolved);
-        psFree(varianceResolved);
-        if (!stack->imageFits->data[i] || !stack->maskFits->data[i] || !stack->varianceFits->data[i]) {
-            psError(PS_ERR_UNKNOWN, false, "Unable to open convolved files %s, %s, %s",
-                    (char*)imageNames->data[i], (char*)maskNames->data[i], (char*)varianceNames->data[i]);
-            return NULL;
-        }
+
+// Open an image
+#define IMAGE_OPEN(NAMES, FITS, INDEX)          \
+        if (NAMES) { \
+            psString resolved = pmConfigConvertFilename((NAMES)->data[INDEX], config, false, false); \
+            (FITS)->data[INDEX] = psFitsOpen(resolved, "r");                            \
+            if (!(FITS)->data[INDEX]) { \
+                psError(PS_ERR_IO, false, "Unable to open file %s", (char*)(NAMES)->data[INDEX]); \
+                psFree(resolved); \
+                return NULL; \
+            } \
+            psFree(resolved); \
+        }
+
+        IMAGE_OPEN(imageNames, stack->imageFits, i);
+        IMAGE_OPEN(maskNames, stack->maskFits, i);
+        IMAGE_OPEN(varianceNames, stack->varianceFits, i);
     }
 
@@ -116 +126 @@
             }
             pmReadout *ro = pmReadoutAlloc(cell); // Readout for thread
-            ro->covariance = psMemIncrRefCounter(covariances->data[j]);
+            if (covariances) {
+                ro->covariance = psMemIncrRefCounter(covariances->data[j]);
+            }
             readouts->data[j] = ro;
         }
@@ -186 +198 @@
                 psFits *varianceFits = stack->varianceFits->data[i]; // FITS file for variance
 
-
-                int zMax = 0;
+                int zMax = 0;
                 bool keepReading = false;
-                if (pmReadoutMore(ro, imageFits, 0, &zMax, rows, config)) {
+
+                if (imageFits && pmReadoutMore(ro, imageFits, 0, &zMax, rows, config)) {
                     keepReading = true;
                     if (!pmReadoutReadChunk(ro, imageFits, 0, NULL, rows, overlap, config)) {
@@ -199 +211 @@
                 }
 
-                if (pmReadoutMoreMask(ro, maskFits, 0, &zMax, rows, config)) {
+                if (maskFits && pmReadoutMoreMask(ro, maskFits, 0, &zMax, rows, config)) {
                     keepReading = true;
                     if (!pmReadoutReadChunkMask(ro, maskFits, 0, NULL, rows, overlap, config)) {
@@ -209 +221 @@
                 }
 
-                if (pmReadoutMoreVariance(ro, varianceFits, 0, &zMax, rows, config)) {
+                if (varianceFits && pmReadoutMoreVariance(ro, varianceFits, 0, &zMax, rows, config)) {
                     keepReading = true;
                     if (!pmReadoutReadChunkVariance(ro, varianceFits, 0, NULL, rows, overlap, config)) {
@@ -263 +275 @@
 
     {
-        psThreadTask *task = psThreadTaskAlloc("PPSTACK_INSPECT", 2);
+        psThreadTask *task = psThreadTaskAlloc("PPSTACK_INSPECT", 3);
         task->function = &ppStackInspect;
         psThreadTaskAdd(task);
@@ -270 +282 @@
 
     {
-        psThreadTask *task = psThreadTaskAlloc("PPSTACK_FINAL_COMBINE", 3);
+        psThreadTask *task = psThreadTaskAlloc("PPSTACK_FINAL_COMBINE", 6);
         task->function = &ppStackReadoutFinalThread;
         psThreadTaskAdd(task);

trunk/ppStack/src/ppStackThread.h

@@ -35 +35 @@
 ppStackThreadData *ppStackThreadDataSetup(
     const ppStackOptions *options,      // Options
-    const pmConfig *config              // Configuration
+    const pmConfig *config,             // Configuration
+    bool conv                           // Use convolved products?
     );