Changeset 26076

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

Files:

• ippconfig/recipes (modified) (1 prop)
• ippconfig/recipes/filerules-mef.mdc (modified) (1 diff)
• ippconfig/recipes/filerules-simple.mdc (modified) (1 diff)
• ippconfig/recipes/filerules-split.mdc (modified) (1 diff)
• ippconfig/recipes/ppStack.config (modified) (2 diffs)
• ippconfig/recipes/ppSub.config (modified) (1 diff)
• ppStack (modified) (1 prop)
• ppStack/src/ppStack.h (modified) (2 diffs)
• ppStack/src/ppStackArguments.c (modified) (3 diffs)
• ppStack/src/ppStackCamera.c (modified) (3 diffs)
• ppStack/src/ppStackCleanup.c (modified) (1 diff)
• ppStack/src/ppStackCombineFinal.c (modified) (3 diffs)
• ppStack/src/ppStackCombineInitial.c (modified) (5 diffs)
• ppStack/src/ppStackCombinePrepare.c (modified) (3 diffs)
• ppStack/src/ppStackConvolve.c (modified) (9 diffs)
• ppStack/src/ppStackFiles.c (modified) (1 diff)
• ppStack/src/ppStackLoop.c (modified) (7 diffs)
• ppStack/src/ppStackLoop.h (modified) (1 diff)
• ppStack/src/ppStackMatch.c (modified) (8 diffs)
• ppStack/src/ppStackOptions.c (modified) (4 diffs)
• ppStack/src/ppStackOptions.h (modified) (2 diffs)
• ppStack/src/ppStackReadout.c (modified) (13 diffs)
• ppStack/src/ppStackReject.c (modified) (7 diffs)
• ppStack/src/ppStackSetup.c (modified) (2 diffs)
• ppStack/src/ppStackSources.c (modified) (1 diff)
• ppStack/src/ppStackThread.c (modified) (9 diffs)
• ppStack/src/ppStackThread.h (modified) (1 diff)
• ppSub/src/ppSubArguments.c (modified) (1 diff)
• psModules/src/camera (modified) (1 prop)
• psModules/src/camera/pmFPAMaskWeight.c (modified) (1 diff)
• psModules/src/camera/pmFPAMaskWeight.h (modified) (1 diff)
• psModules/src/camera/pmReadoutFake.c (modified) (1 prop)
• psModules/src/imcombine/pmStack.c (modified) (23 diffs, 1 prop)
• psModules/src/imcombine/pmStack.h (modified) (1 diff, 1 prop)
• psModules/src/objects (modified) (1 prop)
• psModules/src/objects/pmSourceMatch.c (modified) (18 diffs)
• psModules/src/objects/pmSourceMatch.h (modified) (5 diffs)

trunk/ippconfig/recipes/filerules-mef.mdc

@@ -263 +263 @@
 PPSTACK.OUTPUT.MASK     OUTPUT {OUTPUT}.mk.fits                  MASK      COMP_MASK  FPA        TRUE      NONE
 PPSTACK.OUTPUT.VARIANCE OUTPUT {OUTPUT}.wt.fits                  VARIANCE  COMP_WT    FPA        TRUE      NONE
+PPSTACK.UNCONV.MASK     OUTPUT {OUTPUT}.unconv.mask.fits         MASK      COMP_MASK  FPA        TRUE      NONE
+PPSTACK.UNCONV.VARIANCE OUTPUT {OUTPUT}.unconv.wt.fits           VARIANCE  COMP_WT    FPA        TRUE      NONE
 PPSTACK.TARGET.PSF      OUTPUT {OUTPUT}.target.psf               PSF       NONE       CHIP       TRUE      NONE
 PPSTACK.CONV.KERNEL     OUTPUT {OUTPUT}.{FILE.INDEX}.kernel      SUBKERNEL NONE       FPA        TRUE      NONE

trunk/ippconfig/recipes/filerules-simple.mdc

@@ -212 +212 @@
 PPSTACK.OUTPUT.MASK     OUTPUT {OUTPUT}.mask.fits    MASK      NONE       FPA        TRUE      NONE
 PPSTACK.OUTPUT.VARIANCE OUTPUT {OUTPUT}.weight.fits  VARIANCE  NONE       FPA        TRUE      NONE
+PPSTACK.UNCONV.MASK     OUTPUT {OUTPUT}.unconv.mask.fits MASK  COMP_MASK  FPA        TRUE      NONE
+PPSTACK.UNCONV.VARIANCE OUTPUT {OUTPUT}.unconv.wt.fits VARIANCE COMP_WT   FPA        TRUE      NONE
 PPSTACK.TARGET.PSF      OUTPUT {OUTPUT}.target.psf   PSF       NONE       CHIP       TRUE      NONE
 PPSTACK.CONV.KERNEL     OUTPUT {OUTPUT}.{FILE.INDEX}.kernel SUBKERNEL NONE FPA       TRUE      NONE

trunk/ippconfig/recipes/filerules-split.mdc

@@ -229 +229 @@
 PPSTACK.OUTPUT.MASK     OUTPUT {OUTPUT}.mask.fits                MASK      COMP_MASK  FPA        TRUE      NONE
 PPSTACK.OUTPUT.VARIANCE OUTPUT {OUTPUT}.wt.fits                  VARIANCE  COMP_WT    FPA        TRUE      NONE
+PPSTACK.UNCONV.MASK     OUTPUT {OUTPUT}.unconv.mask.fits         MASK      COMP_MASK  FPA        TRUE      NONE
+PPSTACK.UNCONV.VARIANCE OUTPUT {OUTPUT}.unconv.wt.fits           VARIANCE  COMP_WT    FPA        TRUE      NONE
 PPSTACK.TARGET.PSF      OUTPUT {OUTPUT}.target.psf               PSF       NONE       CHIP       TRUE      NONE
 PPSTACK.CONV.KERNEL     OUTPUT {OUTPUT}.{FILE.INDEX}.kernel      SUBKERNEL NONE       FPA        TRUE      NONE

trunk/ippconfig/recipes/ppStack.config

@@ -2 +2 @@
 
 CONVOLVE        BOOL    TRUE            # Convolve images when stacking?
-ITER            S32     1               # Number of rejection iterations
+COMBINE.ITER    F32     0.5             # Number of rejection iterations per input
 COMBINE.REJ     F32     2.5             # Rejection threshold in combination (sigma)
-COMBINE.SYS     F32     0.03            # Relative systematic error in combination
+COMBINE.SYS     F32     0.1             # Relative systematic error in combination
 COMBINE.DISCARD F32     0.2             # Discard fraction for Olympic weighted mean
-MASK.VAL        STR     MASK.VALUE,CONV.BAD,BURNTOOL    # Mask value of input bad pixels
+MASK.IN         STR     MASK.VALUE,CONV.BAD     # Mask value of input bad pixels
+MASK.SUSPECT    STR     SUSPECT,CONV.POOR       # Mask value of suspect pixels
 MASK.BAD        STR     BLANK           # Mask value to give bad pixels
 MASK.POOR       STR     CONV.POOR       # Mask value to give poor pixels
@@ -33 +34 @@
 PHOT.FRAC       F32     0.9             # Minimum fraction of good pixels
 
-ZP.RADIUS       F32     1.0             # Radius (pixels) for matching sources
+ZP.RADIUS       F32     4.0             # Radius (pixels) for matching sources
 ZP.ITER.1       S32     5               # Iterations for zero point; pass 1
 ZP.ITER.2       S32     1000            # Iterations for zero point; pass 2

trunk/ippconfig/recipes/ppSub.config

@@ -15 +15 @@
 SYS.ERR         F32     0.1             # Relative systematic error in images
 
-MASK.IN         STR     MASK.VALUE,CONV.BAD     # Mask value for input
+MASK.VAL        STR     MASK.VALUE,CONV.BAD     # Mask value for input
 MASK.BAD        STR     CONV.BAD        # Mask value to give bad pixels
 MASK.POOR       STR     CONV.POOR       # Mask value to give poor pixels

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?
     );
 

trunk/ppSub/src/ppSubArguments.c

@@ -86 +86 @@
     psMetadataAddS32(arguments, PS_LIST_TAIL, "-convolve", 0, "Image to convolve [1 or 2]", 0);
     psMetadataAddBool(arguments, PS_LIST_TAIL, "-photometry", 0, "Perform photometry?", NULL);
+    psMetadataAddF32(arguments, PS_LIST_TAIL, "-zp", 0, "Zero point for photometry", NAN);
     psMetadataAddBool(arguments, PS_LIST_TAIL, "-inverse", 0, "Generate inverse subtractions?", NULL);
     psMetadataAddBool(arguments, PS_LIST_TAIL, "-visual", 0, "Show diagnostic plots", NULL);

trunk/psModules/src/camera/pmFPAMaskWeight.c

@@ -444 +444 @@
     }
 
-    return true;
+    return psMetadataAddF32(readout->analysis, PS_LIST_TAIL, PM_READOUT_ANALYSIS_RENORM, 0,
+                            "Renormalisation of variance", PS_SQR(correction));
 }
 

trunk/psModules/src/camera/pmFPAMaskWeight.h

@@ -15 +15 @@
 /// @addtogroup Camera Camera Layout
 /// @{
+
+#define PM_READOUT_ANALYSIS_RENORM "READOUT.RENORM" // Name on analysis metadata for renormalisation
 
 /// Set a temporary readout mask using CELL.SATURATION and CELL.BAD

trunk/psModules/src/imcombine/pmStack.c

@@ -33 +33 @@
 #define NUM_DIRECT_STDEV 5              // For less than this number of values, measure stdev directly
 
+
 //#define TESTING                         // Enable test output
-//#define TEST_X 1085                     // x coordinate to examine
-//#define TEST_Y 3371                     // y coordinate to examine
+//#define TEST_X 3122-1                     // x coordinate to examine
+//#define TEST_Y 1028-1                     // y coordinate to examine
 
 
@@ -42 +43 @@
 typedef struct {
     psVector *pixels;                   // Pixel values
-    psVector *masks;                    // Pixel masks
     psVector *variances;                // Pixel variances
     psVector *weights;                  // Pixel weightings
     psVector *sources;                  // Pixel sources (which image did they come from?)
     psVector *limits;                   // Rejection limits
+    psVector *suspects;                 // Pixel is suspect?
     psVector *sort;                     // Buffer for sorting (to get a robust estimator of the standard dev)
 } combineBuffer;
@@ -53 +54 @@
 {
     psFree(buffer->pixels);
-    psFree(buffer->masks);
     psFree(buffer->variances);
     psFree(buffer->weights);
     psFree(buffer->sources);
     psFree(buffer->limits);
+    psFree(buffer->suspects);
     psFree(buffer->sort);
     return;
@@ -69 +70 @@
 
     buffer->pixels = psVectorAlloc(numImages, PS_TYPE_F32);
-    buffer->masks = psVectorAlloc(numImages, PS_TYPE_VECTOR_MASK);
     buffer->variances = psVectorAlloc(numImages, PS_TYPE_F32);
     buffer->weights = psVectorAlloc(numImages, PS_TYPE_F32);
     buffer->sources = psVectorAlloc(numImages, PS_TYPE_U16);
     buffer->limits = psVectorAlloc(numImages, PS_TYPE_F32);
+    buffer->suspects = psVectorAlloc(numImages, PS_TYPE_U8);
     buffer->sort = psVectorAlloc(numImages, PS_TYPE_F32);
     return buffer;
@@ -143 +144 @@
                                    float *stdev, // Standard deviation value, to return
                                    const psVector *values, // Values to combine
-                                   const psVector *masks, // Mask to apply
                                    psVector *sortBuffer // Buffer for sorting
                                    )
 {
     assert(values);
-    assert(!masks || values->n == masks->n);
     assert(values->type.type == PS_TYPE_F32);
-    assert(!masks || masks->type.type == PS_TYPE_VECTOR_MASK);
     assert(sortBuffer && sortBuffer->nalloc >= values->n && sortBuffer->type.type == PS_TYPE_F32);
 
-    // Need to filter out clipped values
-    int num = 0;            // Number of valid values
-    for (int i = 0; i < values->n; i++) {
-        if (!masks || !masks->data.PS_TYPE_VECTOR_MASK_DATA[i]) {
-            sortBuffer->data.F32[num++] = values->data.F32[i];
-        }
-    }
-    sortBuffer->n = num;
-    if (!psVectorSortInPlace(sortBuffer)) {
+    int num = values->n;                // Number of values
+    sortBuffer = psVectorSortIndex(sortBuffer, values);
+    if (!sortBuffer) {
+        *median = NAN;
+        *stdev = NAN;
         return false;
     }
@@ -167 +161 @@
     if (num == 3) {
         // Attempt to measure standard deviation with only three values (and one of those possibly corrupted)
-        *median = sortBuffer->data.F32[1];
+        *median = values->data.F32[sortBuffer->data.S32[1]];
         if (stdev) {
-            float diff1 = sortBuffer->data.F32[0] - *median;
-            float diff2 = sortBuffer->data.F32[2] - *median;
+            float diff1 = values->data.F32[sortBuffer->data.S32[0]] - *median;
+            float diff2 = values->data.F32[sortBuffer->data.S32[2]] - *median;
             // This factor of sqrt(2) might not be exact, but it's about right
             *stdev = M_SQRT2 * PS_MIN(fabsf(diff1), fabsf(diff2));
         }
     } else {
-        *median = num % 2 ? sortBuffer->data.F32[num / 2] :
-            (sortBuffer->data.F32[num / 2 - 1] + sortBuffer->data.F32[num / 2]) / 2.0;
+        *median = num % 2 ? values->data.F32[sortBuffer->data.S32[num / 2]] :
+            (values->data.F32[sortBuffer->data.S32[num / 2 - 1]] +
+             values->data.F32[sortBuffer->data.S32[num / 2]]) / 2.0;
         if (stdev) {
             if (num <= NUM_DIRECT_STDEV) {
@@ -182 +177 @@
                 double sum = 0.0;
                 for (int i = 0; i < num; i++) {
-                    sum += PS_SQR(sortBuffer->data.F32[i] - *median);
+                    sum += PS_SQR(values->data.F32[sortBuffer->data.S32[i]] - *median);
                 }
                 *stdev = sqrt(sum / (double)(num - 1));
             } else {
                 // Standard deviation from the interquartile range
-                *stdev = 0.74 * (sortBuffer->data.F32[(int)(0.75 * num)] -
-                                 sortBuffer->data.F32[(int)(0.25 * num)]);
+                *stdev = 0.74 * (values->data.F32[sortBuffer->data.S32[(int)(0.75 * num)]] -
+                                 values->data.F32[sortBuffer->data.S32[(int)(0.25 * num)]]);
             }
         }
@@ -195 +190 @@
     return true;
 }
-
-#if 0
-// Return the weighted median for the pixels
-// This does not appear to produce as clean images as the weighted Olympic mean
-static float combinationWeightedMedian(const psVector *values, // Values to combine
-                                       const psVector *weights, // Weights to combine
-                                       const psVector *masks, // Mask to apply
-                                       psVector *sortBuffer // Buffer for sorting
-                                       )
-{
-    double sumWeight = 0.0;             // Sum of weights
-    for (int j = 0; j < values->n; j++) {
-        if (masks->data.PS_TYPE_VECTOR_MASK_DATA[j]) {
-            continue;
-        }
-        sumWeight += weights->data.F32[j];
-    }
-
-    sortBuffer = psVectorSortIndex(sortBuffer, values);
-    double target = sumWeight / 2.0;    // Target weight
-
-    int dominant = -1;                  // Index of dominant value, if any
-    double cumulativeWeight = 0.0;      // Sum of weights
-    for (int j = 0; j < values->n; j++) {
-        if (masks->data.PS_TYPE_VECTOR_MASK_DATA[j]) {
-            continue;
-        }
-        int index = sortBuffer->data.S32[j];  // Index of value of interest
-        float weight = weights->data.F32[index]; // Weight for value of interest
-        if (weight >= target) {
-            // Get the weighted median of the rest
-            dominant = index;
-            sumWeight -= weight;
-            target = sumWeight / 2.0;
-            continue;
-        }
-        cumulativeWeight += weight;
-        if (cumulativeWeight >= target) {
-            float median = values->data.F32[index]; // Weighted median median
-            if (dominant != -1) {
-                // In the case that a single value contains a disproportionate weight compared to the rest,
-                // we use a weighted mean between that dominant value and the weighted median of the rest.
-                return (values->data.F32[dominant] * weights->data.F32[dominant] + median * sumWeight) /
-                    (weights->data.F32[dominant] + sumWeight);
-            }
-            return median;
-        }
-    }
-
-    return NAN;
-}
-#endif
 
 // Return the weighted Olympic mean for the pixels
 static float combinationWeightedOlympic(const psVector *values, // Values to combine
                                         const psVector *weights, // Weights to combine
-                                        const psVector *masks, // Mask to apply
                                         float frac, // Fraction to discard
                                         psVector *sortBuffer // Buffer for sorting
                                         )
 {
-    int numGood = 0;                    // Number of good values
-    for (int i = 0; i < values->n; i++) {
-        if (masks->data.PS_TYPE_VECTOR_MASK_DATA[i]) {
-            continue;
-        }
-        numGood++;
-    }
+    int numGood = values->n;            // Number of good values
 
     int numBad = frac * numGood + 0.5;  // Number of bad values
@@ -272 +208 @@
     for (int i = 0, j = 0; i < values->n; i++) {
         int index = sortBuffer->data.S32[i]; // Index of interest
-        if (masks->data.PS_TYPE_VECTOR_MASK_DATA[index]) {
-            continue;
-        }
         j++;
         if (j > high) {
@@ -290 +223 @@
 
 // Mark a pixel for inspection
-static inline void combineInspect(const psArray *inputs, // Stack data
-                                  int x, int y, // Pixel
-                                  int source // Source image index
-                                  )
-{
+// Value in pixel doesn't seem to agree with the stack, so need to look closer
+static inline void combineMarkInspect(const psArray *inputs, // Stack data
+                                      int x, int y, // Pixel
+                                      int source // Source image index
+                                      )
+{
+#ifdef TESTING
+    if (x == TEST_X && y == TEST_Y) {
+        fprintf(stderr, "Marking image %d for inspection\n", source);
+    }
+#endif
     pmStackData *data = inputs->data[source]; // Stack data of interest
     if (!data) {
@@ -304 +243 @@
 }
 
-// Given a stack of images, combine with optional rejection.
-// Pixels in the stack that are rejected are marked for subsequent inspection
-static void combinePixels(psImage *image, // Combined image, for output
-                          psImage *mask, // Combined mask, for output
-                          psImage *variance, // Combined variance map, for output
-                          const psArray *inputs, // Stack data
-                          const psVector *weights, // Global (single value) weights for data, or NULL
-                          const psVector *addVariance, // Additional variance for data
-                          const psVector *reject, // Indices of pixels to reject, or NULL
-                          int x, int y, // Coordinates of interest; frame of output image
-                          psImageMaskType maskVal, // Value to mask
-                          psImageMaskType bad, // Value to give bad pixels
-                          int numIter, // Number of rejection iterations
-                          float rej, // Number of standard deviations at which to reject
-                          float sys,    // Relative systematic error
-                          float discard,// Fraction of values to discard (Olympic weighted mean)
-                          bool useVariance, // Use variance for rejection when combining?
-                          bool safe,    // Combine safely?
-                          bool rejectInspect, // Reject values marked for inspection from combination?
-                          combineBuffer *buffer // Buffer for combination; to avoid multiple allocations
-                         )
+// Mark a pixel for rejection
+// Cannot possibly inspect this pixel and confirm that it's good.
+// e.g., Only a single input
+static inline void combineMarkReject(const psArray *inputs, // Stack data
+                                     int x, int y, // Pixel
+                                     int source // Source image index
+                                     )
+{
+#ifdef TESTING
+    if (x == TEST_X && y == TEST_Y) {
+        fprintf(stderr, "Marking pixel image %d for rejection\n", source);
+    }
+#endif
+    pmStackData *data = inputs->data[source]; // Stack data of interest
+    if (!data) {
+        psWarning("Can't find input data for source %d", source);
+        return;
+    }
+    data->reject = psPixelsAdd(data->reject, data->reject->nalloc, x, y);
+    return;
+}
+
+
+// Extract vectors for simple combination/rejection operations
+static void combineExtract(int *num,                        // Number of good pixels
+                           bool *suspect,                   // Any suspect pixels?
+                           combineBuffer *buffer, // Buffer with vectors
+                           psImage *image, // Combined image, for output
+                           psImage *mask, // Combined mask, for output
+                           psImage *variance, // Combined variance map, for output
+                           const psArray *inputs, // Stack data
+                           const psVector *weights, // Global (single value) weights for data, or NULL
+                           const psVector *addVariance, // Additional variance for data
+                           const psVector *reject, // Indices of pixels to reject, or NULL
+                           int x, int y, // Coordinates of interest; frame of output image
+                           psImageMaskType maskVal, // Value to mask
+                           psImageMaskType maskSuspect // Value to suspect
+                           )
 {
     // Rudimentary error checking
+    assert(buffer);
     assert(image);
     assert(mask);
     assert(inputs);
-    assert(numIter >= 0);
-    assert(buffer);
-    assert(addVariance);
-    assert((useVariance && variance) || !useVariance);
 
     psVector *pixelData = buffer->pixels; // Values for the pixel of interest
-    psVector *pixelMasks = buffer->masks; // Masks for the pixel of interest
     psVector *pixelVariances = variance ? buffer->variances : NULL; // Variances for the pixel of interest
     psVector *pixelWeights = buffer->weights; // Image weights for the pixel of interest
     psVector *pixelSources = buffer->sources; // Sources for the pixel of interest
     psVector *pixelLimits = buffer->limits; // Limits for the pixel of interest
-    psVector *sort = buffer->sort;      // Sort buffer
+    psVector *pixelSuspects = buffer->suspects; // Is the pixel suspect?
+
+    if (suspect) {
+        *suspect = false;
+    }
 
     // Extract the pixel and mask data
-    int num = 0;                        // Number of good images
+    int numGood = 0;                    // Number of good pixels
     for (int i = 0, j = 0; i < inputs->n; i++) {
         // Check if this pixel has been rejected.  Assumes that the rejection pixel list is sorted --- it
@@ -364 +320 @@
         }
 
+        pixelSuspects->data.U8[numGood] = mask->data.PS_TYPE_IMAGE_MASK_DATA[yIn][xIn] & maskSuspect ?
+            true : false;
+
         psImage *image = data->readout->image; // Image of interest
         psImage *variance = data->readout->variance; // Variance map of interest
-        pixelData->data.F32[num] = image->data.F32[yIn][xIn];
+        pixelData->data.F32[numGood] = image->data.F32[yIn][xIn];
         if (variance) {
-            pixelVariances->data.F32[num] = variance->data.F32[yIn][xIn] * addVariance->data.F32[i];
-        }
-        pixelWeights->data.F32[num] = data->weight;
-        pixelSources->data.U16[num] = i;
-        num++;
-    }
-    pixelData->n = num;
+            pixelVariances->data.F32[numGood] = variance->data.F32[yIn][xIn] * addVariance->data.F32[i];
+        }
+        pixelWeights->data.F32[numGood] = data->weight;
+        pixelSources->data.U16[numGood] = i;
+        numGood++;
+    }
+    pixelData->n = numGood;
     if (variance) {
-        pixelVariances->n = num;
-    }
-    pixelWeights->n = num;
-    pixelSources->n = num;
-    pixelLimits->n = num;
+        pixelVariances->n = numGood;
+    }
+    pixelWeights->n = numGood;
+    pixelSources->n = numGood;
+    pixelLimits->n = numGood;
+    pixelSuspects->n = numGood;
+    *num = numGood;
 
 #ifdef TESTING
     if (x == TEST_X && y == TEST_Y) {
-        for (int i = 0; i < num; i++) {
-            fprintf(stderr, "Input %d (%" PRIu16 "): %f %f %f\n",
+        for (int i = 0; i < numGood; i++) {
+            fprintf(stderr, "Input %d (%" PRIu16 "): %f %f (%f) %f %d\n",
                     i, pixelSources->data.U16[i], pixelData->data.F32[i], pixelVariances->data.F32[i],
-                    pixelWeights->data.F32[i]);
-        }
-    }
-#endif
-
-    // The sensible thing to do varies according to how many good pixels there are.
+                    addVariance->data.F32[i], pixelWeights->data.F32[i], pixelSuspects->data.U8[i]);
+        }
+    }
+#endif
+
+    return;
+}
+
+
+// Combine pixels
+static void combinePixels(psImage *image, // Combined image, for output
+                          psImage *mask, // Combined mask, for output
+                          psImage *variance, // Combined variance map, for output
+                          int num,      // Number of good pixels
+                          combineBuffer *buffer, // Buffer with vectors
+                          int x, int y, // Coordinates of interest; frame of output image
+                          psImageMaskType bad, // Value for bad pixels
+                          bool safe             // Safe combination?
+                          )
+{
+    psVector *pixelData = buffer->pixels; // Values for the pixel of interest
+    psVector *pixelVariances = variance ? buffer->variances : NULL; // Variances for the pixel of interest
+    psVector *pixelWeights = buffer->weights; // Image weights for the pixel of interest
+
     // Default option is that the pixel is bad
     float imageValue = NAN, varianceValue = NAN; // Value for combined image and variance map
     psImageMaskType maskValue = bad;    // Value for combined mask
     switch (num) {
-      case 0:
-        // Nothing to combine: it's bad
-#ifdef TESTING
-    if (x == TEST_X && y == TEST_Y) {
-        fprintf(stderr, "No inputs to combine, pixel is bad.\n");
-    }
-#endif
-        break;
+      case 0: {
+          // Nothing to combine: it's bad
+#ifdef TESTING
+          if (x == TEST_X && y == TEST_Y) {
+              fprintf(stderr, "No inputs to combine, pixel is bad.\n");
+          }
+#endif
+          break;
+      }
       case 1: {
           // Accept the single pixel unless we have to be safe
@@ -427 +407 @@
       }
       case 2: {
-          // Accept the mean of the pixels only if we're going to reject based on the variance, or we're not
-          // playing safe
-          if (useVariance || !safe) {
+          // Automatically accept the mean of the pixels only if we're not playing safe
+          if (!safe) {
               float mean, var;   // Mean and variance from combination
               if (combinationMeanVariance(&mean, &var, pixelData, pixelVariances, pixelWeights)) {
@@ -437 +416 @@
 #ifdef TESTING
                   if (x == TEST_X && y == TEST_Y) {
-                      fprintf(stderr, "Two inputs to combine using variance/unsafe --> %f %f\n",
-                              mean, var);
+                      fprintf(stderr, "Two inputs to combine using unsafe --> %f %f\n", mean, var);
                   }
 #endif
               }
           }
-          if (useVariance && numIter > 0) {
-              // Use variance to check that the two are consistent
-              float diff = 0.5 * (pixelData->data.F32[0] - pixelData->data.F32[1]); // Mean flux difference
-              float var1 = pixelVariances->data.F32[0]; // Variance of first
-              float var2 = pixelVariances->data.F32[1]; // Variance of second
-              // Systematic error contributes to the rejection level
-              var1 += PS_SQR(sys * pixelData->data.F32[0]);
-              var2 += PS_SQR(sys * pixelData->data.F32[1]);
-
-              float sigma2 = var1 + var2; // Combined variance
-              if (PS_SQR(diff) > PS_SQR(rej) * sigma2) {
-                  // Not consistent: mark both for inspection
-                  if (rejectInspect) {
-                      imageValue = NAN;
-                      varianceValue = NAN;
-                      maskValue = bad;
-                  } else {
-                      combineInspect(inputs, x, y, pixelSources->data.U16[0]);
-                      combineInspect(inputs, x, y, pixelSources->data.U16[1]);
-                  }
-#ifdef TESTING
-                  if (x == TEST_X && y == TEST_Y) {
-                      fprintf(stderr, "Both pixels marked for inspection (%f > %f x %f\n)",
-                              diff, rej, sqrtf(sigma2));
-                  }
-#endif
+#ifdef TESTING
+          else {
+              if (x == TEST_X && y == TEST_Y) {
+                  fprintf(stderr, "Two inputs to combine, safety on, pixel is bad\n");
               }
           }
+#endif
           break;
       }
       default: {
-          // Record the value derived with no clipping, because pixels rejected using the harsh clipping
-          // applied in the first pass might later be accepted.
+          // Can combine without too much worrying
           float mean, var;           // Mean and variance of the combination
           if (!combinationMeanVariance(&mean, &var, pixelData, pixelVariances, pixelWeights)) {
@@ -488 +444 @@
           }
 #endif
-
-          // Prepare for clipping iteration
-          if (numIter > 0) {
-              pixelMasks->n = num;
-              psVectorInit(pixelMasks, 0);
-              if (useVariance) {
-                  // Convert to rejection limits --- saves doing it later.
-                  // Using squared rejection limit because it's cheaper than sqrts
-                  float rej2 = PS_SQR(rej); // Rejection level squared
-                  double sumWeights = 0.0;
-                  for (int i = 0; i < num; i++) {
-                      sumWeights += pixelWeights->data.F32[i];
-                  }
-                  for (int i = 0; i < num; i++) {
-                      // Systematic error contributes to the rejection level
-                      float sysVar = PS_SQR(sys * pixelData->data.F32[i]);
-#ifdef TESTING
-                      // Correct variance for comparison against weighted mean including itself
-                      float compare = 1.0 - pixelWeights->data.F32[i] / sumWeights;
-                      if (x == TEST_X && y == TEST_Y) {
-                          fprintf(stderr, "Variance %d (%d): %f %f %f\n", i, pixelSources->data.U16[i],
-                                  pixelVariances->data.F32[i], sysVar, compare);
-                      }
-#endif
-
-                      pixelLimits->data.F32[i] = rej2 * (pixelVariances->data.F32[i] + sysVar);
-                  }
-              }
-          }
-
-          // The clipping that follows is solely to identify suspect pixels.
-          // These suspect pixels will be inspected in more detail by other functions.
-          int numClipped = INT_MAX;     // Number of pixels clipped per iteration
-          int totalClipped = 0;         // Total number of pixels clipped
-          for (int i = 0; i < numIter && numClipped > 0 && num - totalClipped > 2; i++) {
-              numClipped = 0;
-              float median = NAN;       // Middle of distribution
-              float limit = NAN;        // Rejection limit
-              if (!useVariance) {
-                  float stdev;  // Median and stdev of the combination, for rejection
-                  if (!combinationMedianStdev(&median, useVariance ? NULL : &stdev,
-                                              pixelData, pixelMasks, sort)) {
-                      psWarning("Bad median/stdev at %d,%d", x, y);
-                      break;
-                  }
-                  limit = rej * stdev;
-#ifdef TESTING
-                  if (x == TEST_X && y == TEST_Y) {
-                      fprintf(stderr, "Rejecting without variance; rejection limit: %f\n", limit);
-                  }
-#endif
-              } else {
-#ifdef TESTING
-                  if (x == TEST_X && y == TEST_Y) {
-                      fprintf(stderr, "Rejecting with variance...\n");
-                  }
-#endif
-                  median = combinationWeightedOlympic(pixelData, pixelWeights, pixelMasks, discard, sort);
-              }
-
-#ifdef TESTING
-              if (x == TEST_X && y == TEST_Y) {
-                  fprintf(stderr, "Median: %f\n", median);
-              }
-#endif
-
-
-// Mask a pixel for inspection
-#define MASK_PIXEL_FOR_INSPECTION() \
-    pixelMasks->data.PS_TYPE_VECTOR_MASK_DATA[j] = 0xff; \
-    if (!rejectInspect) { \
-        combineInspect(inputs, x, y, pixelSources->data.U16[j]); \
-    } \
-    numClipped++; \
-    totalClipped++;
-
-              for (int j = 0; j < num; j++) {
-                  if (pixelMasks->data.PS_TYPE_VECTOR_MASK_DATA[j]) {
-                      continue;
-                  }
-                  float diff = pixelData->data.F32[j] - median; // Difference from expected
-                  if (useVariance) {
-                      // Comparing squares --- cheaper than lots of sqrts
-                      // pixelVariances includes the rejection limit, from above
-                      if (PS_SQR(diff) > pixelLimits->data.F32[j]) {
-                          MASK_PIXEL_FOR_INSPECTION();
-#ifdef TESTING
-                          if (x == TEST_X && y == TEST_Y) {
-                              fprintf(stderr, "Rejecting input %d based on variance: %f > %f\n",
-                                      j, diff, sqrtf(pixelLimits->data.F32[j]));
-                          }
-#endif
-                      }
-                  } else if (fabsf(diff) > limit) {
-                      MASK_PIXEL_FOR_INSPECTION();
-#ifdef TESTING
-                      if (x == TEST_X && y == TEST_Y) {
-                          fprintf(stderr, "Rejecting input %d based on distribution: %f > %f\n",
-                                  j, diff, limit);
-                      }
-#endif
-                  }
-              }
-          }
-
-          if (rejectInspect && totalClipped > 0) {
-              // Get rid of the masked values
-              // The alternative to this is to make combinationMeanVariance() accept a mask
-              int good = 0;            // Index of good value
-              for (int i = 0; i < num; i++) {
-                  if (pixelMasks->data.PS_TYPE_VECTOR_MASK_DATA[i]) {
-                      continue;
-                  }
-                  if (i != good) {
-                      pixelData->data.F32[good] = pixelData->data.F32[i];
-                      pixelVariances->data.F32[good] = pixelVariances->data.F32[i];
-                      pixelWeights->data.F32[good] = pixelWeights->data.F32[i];
-                      pixelData->data.F32[good] = pixelData->data.F32[i];
-                  }
-                  good++;
-              }
-              pixelData->n = good;
-              pixelVariances->n = good;
-              pixelWeights->n = good;
-              if (combinationMeanVariance(&mean, &var, pixelData, pixelVariances, pixelWeights)) {
-                  imageValue = mean;
-                  varianceValue = var;
-                  maskValue = 0;
-              } else {
-                  imageValue = NAN;
-                  varianceValue = NAN;
-                  maskValue = bad;
-              }
-          }
-
           break;
       }
@@ -633 +454 @@
     }
 
+#ifdef TESTING
+    mask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] = num;
+#endif
+
+
     return;
+}
+
+
+// Test pixels to be combined
+// Returns false to repeat without suspect pixels
+static bool combineTest(int num,      // Number of good pixels
+                        bool suspect, // Does the stack contain suspect pixels?
+                        psArray *inputs,       // Original inputs (for flagging)
+                        combineBuffer *buffer, // Buffer with vectors
+                        int x, int y, // Coordinates of interest; frame of output image
+                        float iter, // Number of rejection iterations per input
+                        float rej, // Number of standard deviations at which to reject
+                        float sys,    // Relative systematic error
+                        float olympic,// Fraction of values to discard (Olympic weighted mean)
+                        bool useVariance, // Use variance for rejection when combining?
+                        bool safe    // Combine safely?
+                        )
+{
+    if (iter <= 0) {
+        return true;
+    }
+
+    int numIter = PS_MAX(iter * num, 1); // Number of iterations
+
+#ifdef TESTING
+    if (x == TEST_X && y == TEST_Y) {
+        fprintf(stderr, "Testing pixel %d,%d: %d %f %f %f %d %d\n",
+                x, y, numIter, rej, sys, olympic, useVariance, safe);
+    }
+#endif
+
+    psVector *pixelData = buffer->pixels; // Values for the pixel of interest
+    psVector *pixelWeights = buffer->weights; // Is the pixel suspect?
+    psVector *pixelVariances = buffer->variances; // Variances for the pixel of interest
+    psVector *pixelSources = buffer->sources; // Sources for the pixel of interest
+    psVector *pixelSuspects = buffer->suspects; // Is the pixel suspect?
+    psVector *pixelLimits = buffer->limits; // Is the pixel suspect?
+
+    // Set up rejection limits
+    float rej2 = PS_SQR(rej); // Rejection level squared
+    if (num > 2 && useVariance) {
+        // Convert rejection limits --- saves doing it later multiple times
+        // Using squared rejection limit because it's cheaper than sqrts
+        double sumWeights = 0.0;
+        for (int i = 0; i < num; i++) {
+            sumWeights += pixelWeights->data.F32[i];
+        }
+        for (int i = 0; i < num; i++) {
+            // Systematic error contributes to the rejection level
+            float sysVar = PS_SQR(sys * pixelData->data.F32[i]);
+#ifdef TESTING
+            // Correct variance for comparison against weighted mean including itself
+            float compare = 1.0 - pixelWeights->data.F32[i] / sumWeights;
+            if (x == TEST_X && y == TEST_Y) {
+                fprintf(stderr, "Variance %d (%d): %f %f %f\n", i, pixelSources->data.U16[i],
+                        pixelVariances->data.F32[i], sysVar, compare);
+            }
+#endif
+            pixelLimits->data.F32[i] = rej2 * (pixelVariances->data.F32[i] + sysVar);
+        }
+    }
+
+    int maskIndex = 0;                  // Index of pixel to mask
+    int totalClipped = 0;               // Total number of pixels clipped
+    for (int i = 0; i < numIter && maskIndex >= 0; i++) {
+        maskIndex = -1;                 // Nothing to reject
+
+        switch (num) {
+          case 0:
+            break;
+          case 1:
+            if (i == 0 && safe) {
+                combineMarkReject(inputs, x, y, pixelSources->data.U16[0]);
+            }
+            break;
+          case 2: {
+              if (useVariance) {
+                  // Use variance to check that the two are consistent
+                  float diff = 0.5 * (pixelData->data.F32[0] - pixelData->data.F32[1]); // Mean flux difference
+                  float var1 = pixelVariances->data.F32[0]; // Variance of first
+                  float var2 = pixelVariances->data.F32[1]; // Variance of second
+                  // Systematic error contributes to the rejection level
+                  var1 += PS_SQR(sys * pixelData->data.F32[0]);
+                  var2 += PS_SQR(sys * pixelData->data.F32[1]);
+
+                  float sigma2 = var1 + var2; // Combined variance
+                  if (PS_SQR(diff) > rej2 * sigma2) {
+                      // Not consistent: don't believe either!
+                      if (i == 0 && suspect) {
+                          combineMarkReject(inputs, x, y, pixelSources->data.U16[0]);
+                          combineMarkReject(inputs, x, y, pixelSources->data.U16[1]);
+                      } else {
+                          combineMarkInspect(inputs, x, y, pixelSources->data.U16[0]);
+                          combineMarkInspect(inputs, x, y, pixelSources->data.U16[1]);
+                      }
+#ifdef TESTING
+                      if (x == TEST_X && y == TEST_Y) {
+                          fprintf(stderr, "Flagged both pixels (%f > %f x %f\n)",
+                                  diff, rej, sqrtf(sigma2));
+                      }
+#endif
+                  }
+              } else if (i == 0 && safe) {
+                  // Can't test them, and we want to be safe, so reject
+                  combineMarkReject(inputs, x, y, pixelSources->data.U16[0]);
+                  combineMarkReject(inputs, x, y, pixelSources->data.U16[1]);
+              }
+              break;
+          }
+#if 0
+          case 3: {
+              // Want to be a bit careful on the rejection than for a larger number of inputs
+              if (!useVariance) {
+                  return combineTestGeneral(num, suspect, inputs, buffer, x, y, numIter, rej, sys,
+                                            olympic, useVariance, safe, allowSuspect);
+              }
+
+              // Differences between pixel values
+              float diff01 = pixelData->data.F32[0] - pixelData->data.F32[1];
+              float diff12 = pixelData->data.F32[1] - pixelData->data.F32[2];
+              float diff20 = pixelData->data.F32[2] - pixelData->data.F32[0];
+              // Variance for each pixel
+              float var0 = pixelVariances->data.F32[0] + PS_SQR(sys * pixelData->data.F32[0]);
+              float var1 = pixelVariances->data.F32[1] + PS_SQR(sys * pixelData->data.F32[1]);
+              float var2 = pixelVariances->data.F32[2] + PS_SQR(sys * pixelData->data.F32[2]);
+              // Errors in pixel differences
+              float err01 = var0 + var1;
+              float err12 = var1 + var2;
+              float err20 = var2 + var0;
+
+#ifdef TESTING
+              if (x == TEST_X && y == TEST_Y) {
+                  fprintf(stderr, "Diff 0-1: %f %f\n", diff01, err01);
+                  fprintf(stderr, "Diff 1-2: %f %f\n", diff12, err12);
+                  fprintf(stderr, "Diff 2-0: %f %f\n", diff20, err20);
+              }
+#endif
+
+              int badPairs = 0;         // Number of bad pairs
+              bool bad01 = false, bad12 = false, bad20 = false; // Pair is bad?
+              if (PS_SQR(diff01) > rej2 * err01) {
+                  bad01 = true;
+                  badPairs++;
+              }
+              if (PS_SQR(diff12) > rej2 * err12) {
+                  bad12 = true;
+                  badPairs++;
+              }
+              if (PS_SQR(diff20) > rej2 * err20) {
+                  bad20 = true;
+                  badPairs++;
+              }
+
+              if (badPairs > 0 && allowSuspect && suspect) {
+                  return false;
+              }
+
+              switch (badPairs) {
+                case 0:
+                  // Nothing to worry about!
+                  break;
+                case 1:
+                  // Can't tell which image is bad, so be sure to get it
+                  if (bad01) {
+                      combineMarkInspect(inputs, x, y, pixelSources->data.U16[0]);
+                      combineMarkInspect(inputs, x, y, pixelSources->data.U16[1]);
+                      break;
+                  }
+                  if (bad12) {
+                      combineMarkInspect(inputs, x, y, pixelSources->data.U16[1]);
+                      combineMarkInspect(inputs, x, y, pixelSources->data.U16[2]);
+                      break;
+                  }
+                  if (bad20) {
+                      combineMarkInspect(inputs, x, y, pixelSources->data.U16[2]);
+                      combineMarkInspect(inputs, x, y, pixelSources->data.U16[0]);
+                      break;
+                  }
+                  psAbort("Should never get here");
+                case 2:
+                  if (bad01 && bad12) {
+                      // 2 and 0 are good
+                      combineMarkInspect(inputs, x, y, pixelSources->data.U16[1]);
+                      break;
+                  }
+                  if (bad12 && bad20) {
+                      // 0 and 1 are good
+                      combineMarkInspect(inputs, x, y, pixelSources->data.U16[2]);
+                      break;
+                  }
+                  if (bad20 && bad01) {
+                      // 1 and 2 are good
+                      combineMarkInspect(inputs, x, y, pixelSources->data.U16[0]);
+                      break;
+                  }
+                  psAbort("Should never get here");
+                case 3:
+                  // Everything's bad
+                  combineMarkInspect(inputs, x, y, pixelSources->data.U16[0]);
+                  combineMarkInspect(inputs, x, y, pixelSources->data.U16[1]);
+                  combineMarkInspect(inputs, x, y, pixelSources->data.U16[2]);
+                  break;
+              }
+              break;
+          }
+#endif
+          default: {
+              if (useVariance) {
+                  float median = combinationWeightedOlympic(pixelData, pixelWeights,
+                                                            olympic, buffer->sort); // Median for stack
+#ifdef TESTING
+                  if (x == TEST_X && y == TEST_Y) {
+                      fprintf(stderr, "Flag with variance, median = %f\n", median);
+                  }
+#endif
+                  float worst = -INFINITY; // Largest deviation
+                  for (int j = 0; j < num; j++) {
+                      float diff = pixelData->data.F32[j] - median; // Difference from expected
+#ifdef TESTING
+                      if (x == TEST_X && y == TEST_Y) {
+                          fprintf(stderr, "Testing input %d: %f\n", j, diff);
+                      }
+#endif
+
+                      // Comparing squares --- cheaper than lots of sqrts
+                      // pixelVariances includes the rejection limit, from above
+                      float diff2 = PS_SQR(diff); // Square difference
+                      if (diff2 > pixelLimits->data.F32[j]) {
+                          float dev = diff2 / pixelLimits->data.F32[j]; // Deviation
+                          if (dev > worst) {
+                              worst = dev;
+                              maskIndex = j;
+                          }
+                      }
+                  }
+              } else {
+                  float median = NAN, stdev = NAN;  // Median and stdev of the combination, for rejection
+                  combinationMedianStdev(&median, &stdev, pixelData, buffer->sort);
+                  float limit = rej * stdev; // Rejection limit
+#ifdef TESTING
+                  if (x == TEST_X && y == TEST_Y) {
+                      fprintf(stderr, "Flag without variance; median = %f, stdev = %f, limit = %f\n",
+                              median, stdev, limit);
+                  }
+#endif
+                  float worst = -INFINITY; // Largest deviation
+                  for (int j = 0; j < num; j++) {
+                      float diff = fabsf(pixelData->data.F32[j] - median); // Difference from expected
+
+                      if (diff > limit) {
+                          float dev = diff / limit; // Deviation
+                          if (dev > worst) {
+                              worst = dev;
+                              maskIndex = j;
+                          }
+                      }
+                  }
+              }
+          }
+        }
+
+        // Do the actual rejection of the pixel
+        if (maskIndex >= 0) {
+            if (suspect) {
+#ifdef TESTING
+                if (x == TEST_X && y == TEST_Y) {
+                    fprintf(stderr, "Throwing out all suspect pixels\n");
+                }
+#endif
+                // Throw out all suspect pixels
+                int numGood = 0;        // Number of good pixels
+                for (int j = 0; j < num; j++) {
+                    if (pixelSuspects->data.U8[j]) {
+                        combineMarkReject(inputs, x, y, pixelSources->data.U16[j]);
+                        continue;
+                    }
+                    if (numGood == j) {
+                        numGood++;
+                        continue;
+                    }
+                    pixelData->data.F32[numGood] = pixelData->data.F32[j];
+                    pixelWeights->data.F32[numGood] = pixelWeights->data.F32[j];
+                    pixelSources->data.U16[numGood] = pixelSources->data.U16[j];
+                    pixelLimits->data.F32[numGood] = pixelLimits->data.F32[j];
+                    pixelVariances->data.F32[numGood] = pixelVariances->data.F32[j];
+                    numGood++;
+                }
+                pixelData->n = numGood;
+                pixelWeights->n = numGood;
+                pixelSources->n = numGood;
+                pixelLimits->n = numGood;
+                pixelVariances->n = numGood;
+                totalClipped += num - numGood;
+                num = numGood;
+                suspect = false;
+            } else {
+                // Throw out masked pixel
+#ifdef TESTING
+                if (x == TEST_X && y == TEST_Y) {
+                    fprintf(stderr, "Throwing out input %d\n", maskIndex);
+                }
+#endif
+                combineMarkInspect(inputs, x, y, pixelSources->data.U16[maskIndex]);
+                int numGood = 0;        // Number of good pixels
+                for (int j = 0; j < num; j++) {
+                    if (j == maskIndex) {
+                        continue;
+                    }
+                    if (numGood == j) {
+                        numGood++;
+                        continue;
+                    }
+                    pixelData->data.F32[numGood] = pixelData->data.F32[j];
+                    pixelWeights->data.F32[numGood] = pixelWeights->data.F32[j];
+                    pixelSources->data.U16[numGood] = pixelSources->data.U16[j];
+                    pixelLimits->data.F32[numGood] = pixelLimits->data.F32[j];
+                    pixelVariances->data.F32[numGood] = pixelVariances->data.F32[j];
+                    numGood++;
+                }
+                pixelData->n = numGood;
+                pixelWeights->n = numGood;
+                pixelSources->n = numGood;
+                pixelLimits->n = numGood;
+                pixelVariances->n = numGood;
+                totalClipped++;
+                num--;
+            }
+        }
+    }
+
+    return true;
 }
 
@@ -639 +796 @@
 // Ensure the input array of pmStackData is valid, and get some details out of it
 static bool validateInputData(bool *haveVariances, // Do we have variance maps in the sky images?
-                              bool *haveRejects, // Do we have lists of rejected pixels?
                               int *num,    // Number of inputs
                               int *numCols, int *numRows, // Size of (sky) images
-                              psArray *input // Input array of pmStackData to validate
+                              const psArray *input, // Input array of pmStackData to validate
+                              const pmReadout *output // Output readout
     )
 {
@@ -668 +825 @@
         PS_ASSERT_IMAGE_TYPE(data->readout->variance, PS_TYPE_F32, false);
     }
-    *haveRejects = (data->reject != NULL);
+    bool haveRejects = (data->reject != NULL); // Do we have rejected pixels?
 
     // Make sure the rest correspond with the first
@@ -681 +838 @@
             return false;
         }
-        if ((*haveRejects && !data->reject) || (data->reject && !*haveRejects)) {
+        if ((haveRejects && !data->reject) || (data->reject && !haveRejects)) {
             psError(PS_ERR_UNEXPECTED_NULL, true,
                     "The rejected pixels are specified in some but not all inputs.");
@@ -697 +854 @@
             PS_ASSERT_IMAGE_TYPE(data->readout->variance, PS_TYPE_F32, false);
         }
+    }
+
+    PM_ASSERT_READOUT_NON_NULL(output, false);
+    if (output->image) {
+        PS_ASSERT_IMAGE_NON_NULL(output->image, false);
+        PS_ASSERT_IMAGE_TYPE(output->image, PS_TYPE_F32, false);
+        PS_ASSERT_IMAGE_NON_NULL(output->mask, false);
+        PS_ASSERT_IMAGE_TYPE(output->mask, PS_TYPE_IMAGE_MASK, false);
+        PS_ASSERT_IMAGES_SIZE_EQUAL(output->image, output->mask, false);
     }
 
@@ -782 +948 @@
 
 /// Stack input images
-bool pmStackCombine(pmReadout *combined, psArray *input, psImageMaskType maskVal, psImageMaskType bad,
-                    int kernelSize, int numIter, float rej, float sys, float discard,
-                    bool entire, bool useVariance, bool safe, bool rejectInspect)
-{
-    PS_ASSERT_PTR_NON_NULL(combined, false);
+bool pmStackCombine(pmReadout *combined, psArray *input, psImageMaskType maskVal, psImageMaskType maskSuspect,
+                    psImageMaskType bad, int kernelSize, float iter, float rej, float sys, float olympic,
+                    bool useVariance, bool safe, bool rejection)
+{
     bool haveVariances;                 // Do we have the variance maps?
-    bool haveRejects;                   // Do we have lists of rejected pixels?
     int num;                            // Number of inputs
     int numCols, numRows;               // Size of (sky) images
-    if (!validateInputData(&haveVariances, &haveRejects, &num, &numCols, &numRows, input)) {
+    if (!validateInputData(&haveVariances, &num, &numCols, &numRows, input, combined)) {
         return false;
     }
     PS_ASSERT_INT_NONNEGATIVE(kernelSize, false);
     PS_ASSERT_INT_POSITIVE(bad, false);
-    PS_ASSERT_INT_NONNEGATIVE(numIter, false);
     if (isnan(rej)) {
-        PS_ASSERT_INT_EQUAL(numIter, 0, false);
+        PS_ASSERT_FLOAT_EQUAL(iter, 0, false);
     } else {
+        PS_ASSERT_FLOAT_LARGER_THAN(iter, 0, false);
         PS_ASSERT_FLOAT_LARGER_THAN(rej, 0.0, false);
-    }
-    if (haveRejects) {
-        // This is a subsequent combination, so expect that the image and mask already exist
-        PS_ASSERT_IMAGE_NON_NULL(combined->image, false);
-        PS_ASSERT_IMAGE_TYPE(combined->image, PS_TYPE_F32, false);
-        PS_ASSERT_IMAGE_NON_NULL(combined->mask, false);
-        PS_ASSERT_IMAGE_TYPE(combined->mask, PS_TYPE_IMAGE_MASK, false);
-        PS_ASSERT_IMAGES_SIZE_EQUAL(combined->image, combined->mask, false);
     }
     if (useVariance && !haveVariances) {
@@ -833 +989 @@
         }
 #endif
-        if (!haveRejects && !data->inspect) {
-            data->inspect = psPixelsAllocEmpty(PIXEL_LIST_BUFFER);
+        if (!rejection) {
+            // Ensure pixels can be put on the appropriate list
+            if (!data->inspect) {
+                data->inspect = psPixelsAllocEmpty(PIXEL_LIST_BUFFER);
+            }
+            if (!data->reject) {
+                data->reject = psPixelsAllocEmpty(PIXEL_LIST_BUFFER);
+            }
         }
     }
@@ -863 +1025 @@
     combineBuffer *buffer = combineBufferAlloc(num);
 
-    if (haveRejects) {
-        psImage *combinedImage = combined->image; // Combined image
-        psImage *combinedMask = combined->mask; // Combined mask
-        psImage *combinedVariance = combined->variance; // Combined variance map
-
-        psArray *pixelMap = pixelMapGenerate(input, minInputCols, maxInputCols,
-                                             minInputRows, maxInputRows); // Map of pixels to source
-        psPixels *pixels = NULL;            // Total list of pixels, with no duplicates
+    // Pull the products out, allocate if necessary
+    psImage *combinedImage = combined->image; // Combined image
+    if (!combinedImage) {
+        combined->image = psImageAlloc(numCols, numRows, PS_TYPE_F32);
+        combinedImage = combined->image;
+    }
+    psImage *combinedMask = combined->mask; // Combined mask
+    if (!combinedMask) {
+        combined->mask = psImageAlloc(numCols, numRows, PS_TYPE_IMAGE_MASK);
+        combinedMask = combined->mask;
+    }
+
+    psImage *combinedVariance = combined->variance; // Combined variance map
+    if (haveVariances && !combinedVariance) {
+        combined->variance = psImageAlloc(numCols, numRows, PS_TYPE_F32);
+        combinedVariance = combined->variance;
+    }
+
+    // Set up rejection list
+    psArray *pixelMap = NULL;           // Map of pixels to source
+    if (rejection) {
+        pixelMap = pixelMapGenerate(input, minInputCols, maxInputCols, minInputRows, maxInputRows);
+    }
+
+    // Combine each pixel
+    for (int y = minInputRows; y < maxInputRows; y++) {
+        for (int x = minInputCols; x < maxInputCols; x++) {
+#ifdef TESTING
+            if (x == TEST_X && y == TEST_Y) {
+                fprintf(stderr, "Combining pixel %d,%d: %x %x %f %f %f %f %d %d %d\n",
+                        x, y, maskVal, bad, iter, rej, sys, olympic, useVariance, safe, rejection);
+            }
+#endif
+            psVector *reject = NULL; // Images to reject for this pixel
+            if (rejection) {
+                reject = pixelMapQuery(pixelMap, minInputCols, minInputRows, x, y);
+#ifdef TESTING
+                if (x == TEST_X && y == TEST_Y) {
+                    fprintf(stderr, "Rejected inputs: ");
+                    if (!reject) {
+                        fprintf(stderr, "<none>\n");
+                    } else {
+                        for (int i = 0; i < reject->n; i++) {
+                            fprintf(stderr, "%d ", reject->data.U16[i]);
+                        }
+                        fprintf(stderr, "\n");
+                    }
+                }
+#endif
+            }
+
+            int num;                    // Number of good pixels
+            bool suspect;               // Suspect pixels in stack?
+            combineExtract(&num, &suspect, buffer, combinedImage, combinedMask, combinedVariance,
+                           input, weights, addVariance, reject, x, y, maskVal, maskSuspect);
+            combinePixels(combinedImage, combinedMask, combinedVariance, num, buffer, x, y, bad, safe);
+
+            if (iter > 0) {
+                combineTest(num, suspect, input, buffer, x, y, iter, rej, sys, olympic,
+                            useVariance, safe);
+            }
+        }
+    }
+
+    psFree(pixelMap);
+    psFree(weights);
+    psFree(buffer);
+
+
+
+#ifndef PS_NO_TRACE
+    if (!rejection && psTraceGetLevel("psModules.imcombine") >= 5) {
         for (int i = 0; i < num; i++) {
-            pmStackData *data = input->data[i]; // Stacking data; contains the list of pixels
-            if (!data) {
+            pmStackData *data = input->data[i]; // Stack data for this input
+            if (!data || !data->inspect) {
                 continue;
             }
-            pixels = psPixelsConcatenate(pixels, data->reject);
-        }
-        pixels = psPixelsDuplicates(pixels, pixels);
-
-        if (entire) {
-            // Combine entire image
-            for (int y = minInputRows; y < maxInputRows; y++) {
-                for (int x = minInputCols; x < maxInputCols; x++) {
-                    psVector *reject = pixelMapQuery(pixelMap, minInputCols, minInputRows,
-                                                     x, y); // Reject these images
-                    combinePixels(combinedImage, combinedMask, combinedVariance, input, weights,
-                                  addVariance, reject, x, y, maskVal, bad, numIter, rej, sys, discard,
-                                  useVariance, safe, rejectInspect, buffer);
-                }
-            }
-        } else {
-            // Only combine previously rejected pixels
-            for (int i = 0; i < pixels->n; i++) {
-                // Pixel coordinates are in the frame of the original image
-                int x = pixels->data[i].x, y = pixels->data[i].y; // Coordinates of interest
-                if (x < minInputCols || x >= maxInputCols || y < minInputRows || y >= maxInputRows) {
-                    continue;
-                }
-                psVector *reject = pixelMapQuery(pixelMap, minInputCols, minInputRows,
-                                                 x, y); // Reject these images
-                combinePixels(combinedImage, combinedMask, combinedVariance, input, weights,
-                              addVariance, reject, x, y, maskVal, bad, numIter, rej, sys, discard,
-                              useVariance, safe, rejectInspect, buffer);
-            }
-        }
-        psFree(pixels);
-        psFree(pixelMap);
-    } else {
-        // Pull the products out, allocate if necessary
-        psImage *combinedImage = combined->image; // Combined image
-        if (!combinedImage) {
-            combined->image = psImageAlloc(numCols, numRows, PS_TYPE_F32);
-            combinedImage = combined->image;
-        }
-        psImage *combinedMask = combined->mask; // Combined mask
-        if (!combinedMask) {
-            combined->mask = psImageAlloc(numCols, numRows, PS_TYPE_IMAGE_MASK);
-            combinedMask = combined->mask;
-        }
-
-        psImage *combinedVariance = combined->variance; // Combined variance map
-        if (haveVariances && !combinedVariance) {
-            combined->variance = psImageAlloc(numCols, numRows, PS_TYPE_F32);
-            combinedVariance = combined->variance;
-        }
-
-        for (int y = minInputRows; y < maxInputRows; y++) {
-            for (int x = minInputCols; x < maxInputCols; x++) {
-                combinePixels(combinedImage, combinedMask, combinedVariance, input, weights,
-                              addVariance, NULL, x, y, maskVal, bad, numIter, rej, sys, discard,
-                              useVariance, safe, rejectInspect, buffer);
-            }
-        }
-
-#ifndef PS_NO_TRACE
-        if (psTraceGetLevel("psModules.imcombine") >= 5) {
-            for (int i = 0; i < num; i++) {
-                pmStackData *data = input->data[i]; // Stack data for this input
-                if (!data || !data->inspect) {
-                    continue;
-                }
-                psTrace("psModules.imcombine", 5, "Image %d: %ld pixels to inspect.\n", i, data->inspect->n);
-            }
-        }
-#endif
-    }
-
-    psFree(weights);
-    psFree(buffer);
+            psTrace("psModules.imcombine", 5, "Image %d: %ld pixels to inspect.\n", i, data->inspect->n);
+        }
+    }
+#endif
 
     return true;

trunk/psModules/src/imcombine/pmStack.h

@@ -44 +44 @@
                     psArray *input,     ///< Input array of pmStackData
                     psImageMaskType maskVal, ///< Mask value of bad pixels
+                    psImageMaskType suspect, ///< Mask value of suspect pixels
                     psImageMaskType bad,     ///< Mask value to give rejected pixels
                     int kernelSize,     ///< Half-size of the convolution kernel
-                    int numIter,        ///< Number of iterations
+                    float iter,         ///< Number of iterations per input
                     float rej,          ///< Rejection limit (standard deviations)
                     float sys,          ///< Relative systematic error
                     float discard,      ///< Fraction of values to discard for Olympic weighted mean
-                    bool entire,        ///< Combine entire image even if rejection lists provided?
                     bool useVariance,   ///< Use variance values for rejection?
                     bool safe,          ///< Play safe with small numbers of input pixels (mask if N <= 2)?

trunk/psModules/src/objects/pmSourceMatch.c

@@ -17 +17 @@
 #define SOURCES_MAX_LEAF 2              // Maximum number of points on a tree leaf
 #define ARRAY_BUFFER 16                 // Buffer for array
+
 
 //////////////////////////////////////////////////////////////////////////////////////////////////////////////
@@ -124 +125 @@
     match->mag = psVectorAllocEmpty(ARRAY_BUFFER, PS_TYPE_F32);
     match->magErr = psVectorAllocEmpty(ARRAY_BUFFER, PS_TYPE_F32);
+    match->x = psVectorAllocEmpty(ARRAY_BUFFER, PS_TYPE_F32);
+    match->y = psVectorAllocEmpty(ARRAY_BUFFER, PS_TYPE_F32);
     match->image = psVectorAllocEmpty(ARRAY_BUFFER, PS_TYPE_U32);
     match->index = psVectorAllocEmpty(ARRAY_BUFFER, PS_TYPE_U32);
@@ -133 +136 @@
 void pmSourceMatchAdd(pmSourceMatch *match, // Match data
                       float mag, float magErr, // Magnitude and error
+                      float x, float y,        // Position
                       int image, // Image index
                       int index // Source index
@@ -141 +145 @@
     match->mag = psVectorExtend(match->mag, match->mag->nalloc, 1);
     match->magErr = psVectorExtend(match->magErr, match->magErr->nalloc, 1);
+    match->x = psVectorExtend(match->x, match->x->nalloc, 1);
+    match->y = psVectorExtend(match->y, match->y->nalloc, 1);
     match->image = psVectorExtend(match->image, match->image->nalloc, 1);
     match->index = psVectorExtend(match->index, match->index->nalloc, 1);
@@ -147 +153 @@
     match->mag->data.F32[num] = mag;
     match->magErr->data.F32[num] = magErr;
+    match->x->data.F32[num] = x;
+    match->y->data.F32[num] = y;
     match->image->data.S32[num] = image;
     match->index->data.S32[num] = index;
@@ -192 +200 @@
                 pmSourceMatch *match = pmSourceMatchAlloc(); // Match data
                 pmSourceMatchAdd(match, magImage->data.F32[j], magErrImage->data.F32[j],
-                                 i, indices->data.S32[j]);
+                                 xImage->data.F32[j], yImage->data.F32[j], i, indices->data.S32[j]);
                 matches->data[j] = match;
             }
@@ -212 +220 @@
                 pmSourceMatch *match = pmSourceMatchAlloc(); // Match data
                 pmSourceMatchAdd(match, magImage->data.F32[j], magErrImage->data.F32[j],
-                                 i, indices->data.S32[j]);
+                                 xImage->data.F32[j], yImage->data.F32[j], i, indices->data.S32[j]);
                 matches->data[k] = match;
             }
@@ -238 +246 @@
                     pmSourceMatch *match = matches->data[index]; // Match data
                     pmSourceMatchAdd(match, magImage->data.F32[j], magErrImage->data.F32[j],
-                                     i, indices->data.S32[j]);
+                                     xImage->data.F32[j], yImage->data.F32[j], i, indices->data.S32[j]);
                     numMatch++;
                 } else {
@@ -244 +252 @@
                     pmSourceMatch *match = pmSourceMatchAlloc(); // Match data
                     pmSourceMatchAdd(match, magImage->data.F32[j], magErrImage->data.F32[j],
-                                     i, indices->data.S32[j]);
+                                     xImage->data.F32[j], yImage->data.F32[j], i, indices->data.S32[j]);
                     xMaster->data.F32[numMaster] = xImage->data.F32[j];
                     yMaster->data.F32[numMaster] = yImage->data.F32[j];
@@ -304 +312 @@
         pmSourceMatch *match = matches->data[i]; // Match of interest
         for (int j = 0; j < match->num && !source; j++) {
-            if (!isfinite(match->mag->data.F32[j]) || !isfinite(match->magErr->data.F32[j])) {
+            if (!isfinite(match->mag->data.F32[j]) || !isfinite(match->magErr->data.F32[j]) ||
+                !isfinite(match->x->data.F32[j]) || !isfinite(match->y->data.F32[j])) {
                 continue;
             }
@@ -365 +374 @@
         double star = 0.0, starErr = 0.0; // Accumulators for star
         for (int j = 0; j < match->num; j++) {
-            if (match->mask->data.PS_TYPE_VECTOR_MASK_DATA[j]) {
+            if (match->mask->data.PS_TYPE_VECTOR_MASK_DATA[j] & PM_SOURCE_MATCH_MASK_PHOT) {
                 continue;
             }
@@ -396 +405 @@
         pmSourceMatch *match = matches->data[i]; // Matched stars
         for (int j = 0; j < match->num; j++) {
-            if (match->mask->data.PS_TYPE_VECTOR_MASK_DATA[j]) {
+            if (match->mask->data.PS_TYPE_VECTOR_MASK_DATA[j] & PM_SOURCE_MATCH_MASK_PHOT) {
                 continue;
             }
@@ -424 +433 @@
         pmSourceMatch *match = matches->data[i]; // Matched stars
         for (int j = 0; j < match->num; j++) {
-            if (match->mask->data.PS_TYPE_VECTOR_MASK_DATA[j]) {
+            if (match->mask->data.PS_TYPE_VECTOR_MASK_DATA[j] & PM_SOURCE_MATCH_MASK_PHOT) {
                 continue;
             }
@@ -543 +552 @@
         pmSourceMatch *match = matches->data[i]; // Matched stars
         for (int j = 0; j < match->num; j++) {
-            if (match->mask->data.PS_TYPE_VECTOR_MASK_DATA[j]) {
+            if (match->mask->data.PS_TYPE_VECTOR_MASK_DATA[j] & PM_SOURCE_MATCH_MASK_PHOT) {
                 continue;
             }
@@ -564 +573 @@
                 if (PS_SQR(dev) > starClip * (PS_SQR(magErr) + sysErr2)) {
                     numRejected++;
-                    match->mask->data.PS_TYPE_VECTOR_MASK_DATA[j] = 0xFF;
+                    match->mask->data.PS_TYPE_VECTOR_MASK_DATA[j] |= PM_SOURCE_MATCH_MASK_PHOT;
                 }
             }
@@ -627 +636 @@
             return NULL;
         }
-        psTrace("psModules.objects", 3, "Pass 1: Determined %d/%d are photometric", numPhoto, numImages);
+        psTrace("psModules.objects", 3, "Pass 1: Determined %d/%d are photometric\n", numPhoto, numImages);
 
         fracRej = sourceMatchRelphotReject(trans, stars, matches, zp, photo, badImage, rej1, sys1);
-        psTrace("psModules.objects", 3, "Pass 1: %f%% of measurements rejected", fracRej * 100);
+        psTrace("psModules.objects", 3, "Pass 1: %f%% of measurements rejected\n", fracRej * 100);
 
         chi2 = sourceMatchRelphotIterate(trans, stars, badImage, matches, zp, photo, sys1);
@@ -649 +658 @@
             return NULL;
         }
-        psTrace("psModules.objects", 3, "Pass 2: Determined %d/%d are photometric", numPhoto, numImages);
+        psTrace("psModules.objects", 3, "Pass 2: Determined %d/%d are photometric\n", numPhoto, numImages);
 
         fracRej = sourceMatchRelphotReject(trans, stars, matches, zp, photo, badImage, rej2, sys2);
-        psTrace("psModules.objects", 3, "Pass 2: %f%% of measurements rejected", fracRej * 100);
+        psTrace("psModules.objects", 3, "Pass 2: %f%% of measurements rejected\n", fracRej * 100);
 
         chi2 = sourceMatchRelphotIterate(trans, stars, badImage, matches, zp, photo, sys2);
@@ -668 +677 @@
     return trans;
 }
+
+
+// Iterate on the star positions and image shifts
+// Returns the solution chi^2
+static float sourceMatchRelastroIterate(psVector *xShift, psVector *yShift, // Shift for image
+                                        psVector *xStar, psVector *yStar,   // Position for star
+                                        const psArray *matches // Array of matches
+    )
+{
+    psAssert(matches, "Need list of matches");
+
+    int numImages = xShift->n;          // Number of images
+    int numStars = matches->n;          // Number of stars
+
+    psAssert(xShift && xShift->type.type == PS_TYPE_F32 && yShift && yShift->type.type == PS_TYPE_F32,
+             "Need shifts");
+    psAssert(yShift->n == numImages, "Not enough shifts: %ld\n", yShift->n);
+    psAssert(xStar && xStar->type.type == PS_TYPE_F32 && yStar && yStar->type.type == PS_TYPE_F32,
+             "Need star positions");
+    psAssert(xStar->n == numStars && yStar->n == numStars, "Not enough stars\n");
+
+    // Solve the star positions
+    psVectorInit(xStar, NAN);
+    psVectorInit(yStar, NAN);
+    psVector *starMask = psVectorAlloc(numStars, PS_TYPE_U8); // Mask for stars
+    psVectorInit(starMask, 0xFF);
+    int numGoodStars = 0;               // Number of stars with good measurements
+    for (int i = 0; i < numStars; i++) {
+        pmSourceMatch *match = matches->data[i]; // Matched stars
+        int numMeasurements = 0;        // Number of unmasked measurements for star
+        double xSum = 0.0, ySum = 0.0;  // Accumulators for star
+        for (int j = 0; j < match->num; j++) {
+            if (match->mask->data.PS_TYPE_VECTOR_MASK_DATA[j] & PM_SOURCE_MATCH_MASK_ASTRO) {
+                continue;
+            }
+            numMeasurements++;
+            int index = match->image->data.U32[j]; // Image index
+
+            xSum += match->x->data.F32[j] - xShift->data.F32[index];
+            ySum += match->y->data.F32[j] - yShift->data.F32[index];
+        }
+        if (numMeasurements > 1) {
+            // It's only a good star (contributing to the chi^2) if there's more than 1 measurement
+            numGoodStars++;
+            xStar->data.F32[i] = xSum / numMeasurements;
+            yStar->data.F32[i] = ySum / numMeasurements;
+            starMask->data.U8[i] = 0;
+        }
+    }
+
+    // Solve for the shifts
+    psVectorInit(xShift, 0.0);
+    psVectorInit(yShift, 0.0);
+    psVector *num = psVectorAlloc(numImages, PS_TYPE_S32);    // Number of stars
+    psVectorInit(num, 0);
+    for (int i = 0; i < numStars; i++) {
+        if (starMask->data.U8[i]) {
+            continue;
+        }
+        pmSourceMatch *match = matches->data[i]; // Matched stars
+        for (int j = 0; j < match->num; j++) {
+            if (match->mask->data.PS_TYPE_VECTOR_MASK_DATA[j] & PM_SOURCE_MATCH_MASK_ASTRO) {
+                continue;
+            }
+            int index = match->image->data.U32[j]; // Image index
+
+            xShift->data.F32[index] += match->x->data.F32[j] - xStar->data.F32[i];
+            yShift->data.F32[index] += match->y->data.F32[j] - yStar->data.F32[i];
+            num->data.S32[index]++;
+        }
+    }
+    for (int i = 0; i < numImages; i++) {
+        xShift->data.F32[i] /= num->data.S32[i];
+        yShift->data.F32[i] /= num->data.S32[i];
+        psTrace("psModules.objects", 3, "Shift for image %d: %f,%f\n",
+                i, xShift->data.F32[i], yShift->data.F32[i]);
+    }
+    psFree(num);
+
+    // Once more through to evaluate chi^2
+    float chi2 = 0.0;                   // chi^2 for iteration
+    int dof = 0;                        // Degrees of freedom
+    for (int i = 0; i < numStars; i++) {
+        pmSourceMatch *match = matches->data[i]; // Matched stars
+        if (starMask->data.U8[i]) {
+            continue;
+        }
+        for (int j = 0; j < match->num; j++) {
+            if (match->mask->data.PS_TYPE_VECTOR_MASK_DATA[j]) {
+                continue;
+            }
+
+            int index = match->image->data.U32[j]; // Image index
+            float dx = match->x->data.F32[j] - xShift->data.F32[index] - xStar->data.F32[i];
+            float dy = match->y->data.F32[j] - yShift->data.F32[index] - yStar->data.F32[i];
+
+            chi2 += PS_SQR(dx) + PS_SQR(dy);
+            dof++;
+        }
+    }
+    dof -= numGoodStars + numImages;
+    chi2 /= dof;
+
+    return chi2;
+}
+
+// Reject star measurements
+// Returns the fraction of measurements that were rejected
+static float sourceMatchRelastroReject(const psVector *xShift, const psVector *yShift, // Shifts for each image
+                                       const psVector *xStar, const psVector *yStar, // Positions for each star
+                                       const psArray *matches, // Array of matches
+                                       float chi2,             // chi^2 from fit
+                                       float rej               // Rejection threshold
+                                )
+{
+    psAssert(matches, "Need list of matches");
+
+    int numImages = xShift->n;          // Number of images
+    int numStars = matches->n;          // Number of stars
+
+    psAssert(xShift && xShift->type.type == PS_TYPE_F32 && yShift && yShift->type.type == PS_TYPE_F32,
+             "Need shifts");
+    psAssert(yShift->n == numImages, "Not enough shifts: %ld\n", yShift->n);
+    psAssert(xStar && xStar->type.type == PS_TYPE_F32 && yStar && yStar->type.type == PS_TYPE_F32,
+             "Need star positions");
+    psAssert(xStar->n == numStars && yStar->n == numStars, "Not enough stars\n");
+
+    int numRejected = 0;                // Number rejected
+    int numMeasurements = 0;            // Number of measurements
+
+    float thresh = PS_SQR(rej) * chi2;    // Threshold for rejection
+
+    for (int i = 0; i < numStars; i++) {
+        pmSourceMatch *match = matches->data[i]; // Matched stars
+        for (int j = 0; j < match->num; j++) {
+            if (match->mask->data.PS_TYPE_VECTOR_MASK_DATA[j] & PM_SOURCE_MATCH_MASK_ASTRO) {
+                continue;
+            }
+            numMeasurements++;
+            int index = match->image->data.U32[j]; // Image index
+
+            float dx = match->x->data.F32[j] - xShift->data.F32[index] - xStar->data.F32[i];
+            float dy = match->y->data.F32[j] - yShift->data.F32[index] - yStar->data.F32[i];
+
+            if (PS_SQR(dx) + PS_SQR(dy) > thresh) {
+                numRejected++;
+                match->mask->data.PS_TYPE_VECTOR_MASK_DATA[j] |= PM_SOURCE_MATCH_MASK_ASTRO;
+            }
+        }
+    }
+
+    return (float)numRejected / (float)numMeasurements;
+}
+
+psArray *pmSourceMatchRelastro(const psArray *matches, // Array of matches
+                               int numImages,          // Number of images
+                               float tol, // Relative tolerance for convergence
+                               int iter1, // Number of iterations for pass 1
+                               float rej1, // Limit on rejection between iterations for pass 1
+                               int iter2, // Number of iterations for pass 2
+                               float rej2, // Limit on rejection between iterations for pass 2
+                               float rejLimit // Limit on rejection between iterations
+    )
+{
+    PS_ASSERT_ARRAY_NON_NULL(matches, NULL);
+
+    int numStars = matches->n;          // Number of stars
+    psVector *xShift = psVectorAlloc(numImages, PS_TYPE_F32); // x shift for each image
+    psVector *yShift = psVectorAlloc(numImages, PS_TYPE_F32); // y shift for each image
+    psVectorInit(xShift, 0.0);
+    psVectorInit(yShift, 0.0);
+    psVector *xStar = psVectorAlloc(numStars, PS_TYPE_F32); // x position for each star
+    psVector *yStar = psVectorAlloc(numStars, PS_TYPE_F32); // y position for each star
+
+    float chi2 = sourceMatchRelastroIterate(xShift, yShift, xStar, yStar, matches); // chi^2 for solution
+    psTrace("psModules.objects", 1, "Initial: chi^2 = %f\n", chi2);
+    float lastChi2 = INFINITY;          // chi^2 on last iteration
+    float fracRej = INFINITY;           // Fraction of measurements rejected
+
+    // In the first passes, the shifts are not well deteremined: use high systematic error and
+    // rejection thresholds
+    for (int i = 0; i < iter1; i++) {
+        fracRej = sourceMatchRelastroReject(xShift, yShift, xStar, yStar, matches, chi2, rej1);
+        psTrace("psModules.objects", 3, "Pass 1: %f%% of measurements rejected\n", fracRej * 100);
+
+        chi2 = sourceMatchRelastroIterate(xShift, yShift, xStar, yStar, matches);
+        psTrace("psModules.objects", 1, "Pass 1: iter = %d: chi^2 = %f rejected = %f\n", i, chi2, fracRej);
+    }
+
+    for (int i = 0; i < iter2 && (fabsf(lastChi2 - chi2) > tol * chi2 || fracRej > rejLimit); i++) {
+        lastChi2 = chi2;
+
+        fracRej = sourceMatchRelastroReject(xShift, yShift, xStar, yStar, matches, chi2, rej2);
+        psTrace("psModules.objects", 3, "Pass 2: %f%% of measurements rejected\n", fracRej * 100);
+
+        chi2 = sourceMatchRelastroIterate(xShift, yShift, xStar, yStar, matches);
+        psTrace("psModules.objects", 1, "Pass 2: iter = %d: chi^2 = %f rejected = %f\n", i, chi2, fracRej);
+    }
+
+    psFree(xStar);
+    psFree(yStar);
+
+    if (fabsf(lastChi2 - chi2) > tol * chi2 || fracRej > rejLimit) {
+        psWarning("Unable to converge to relphot solution (%f,%f)", (lastChi2 - chi2) / chi2, fracRej);
+    }
+
+    psArray *results = psArrayAlloc(numImages); // Array of results
+    for (int i = 0; i < numImages; i++) {
+        psVector *offset = results->data[i] = psVectorAlloc(2, PS_TYPE_F32); // Offset for image
+        offset->data.F32[0] = xShift->data.F32[i];
+        offset->data.F32[1] = yShift->data.F32[i];
+    }
+    psFree(xShift);
+    psFree(yShift);
+
+    return results;
+}
+

trunk/psModules/src/objects/pmSourceMatch.h

@@ -1 +1 @@
 #ifndef PM_SOURCE_MATCH_H
 #define PM_SOURCE_MATCH_H
+
+#include <pslib.h>
+
+/// Mask values for matched sources
+typedef enum {
+    PM_SOURCE_MATCH_MASK_PHOT = 0x01,   // Source was rejected from photometry fit
+    PM_SOURCE_MATCH_MASK_ASTRO = 0x02,     // Source was rejected from astrometry fit
+} pmSourceMatchMask;
 
 /// Matched sources
@@ -12 +20 @@
     psVector *mag;                      // Magnitudes
     psVector *magErr;                   // Magnitude errors
+    psVector *x, *y;                    // Positions
     psVector *mask;                     // Mask for measurements
 } pmSourceMatch;
@@ -26 +35 @@
     PS_ASSERT_VECTOR_NON_NULL((MATCH)->mag, RVAL); \
     PS_ASSERT_VECTOR_NON_NULL((MATCH)->magErr, RVAL); \
+    PS_ASSERT_VECTOR_NON_NULL((MATCH)->x, RVAL); \
+    PS_ASSERT_VECTOR_NON_NULL((MATCH)->y, RVAL); \
     PS_ASSERT_VECTOR_SIZE((MATCH)->image, (MATCH)->num, RVAL); \
     PS_ASSERT_VECTOR_SIZE((MATCH)->index, (MATCH)->num, RVAL); \
     PS_ASSERT_VECTOR_SIZE((MATCH)->mag, (MATCH)->num, RVAL); \
     PS_ASSERT_VECTOR_SIZE((MATCH)->magErr, (MATCH)->num, RVAL); \
+    PS_ASSERT_VECTOR_SIZE((MATCH)->x, (MATCH)->num, RVAL); \
+    PS_ASSERT_VECTOR_SIZE((MATCH)->y, (MATCH)->num, RVAL); \
 }
 
@@ -38 +51 @@
 void pmSourceMatchAdd(pmSourceMatch *match, // Match data
                       float mag, float magErr, // Magnitude and error
+                      float x, float y,        // Position
                       int image, // Image index
                       int index // Source index
@@ -75 +89 @@
     );
 
+/// Perform relative astrometry to calibrate images
+psArray *pmSourceMatchRelastro(const psArray *matches, // Array of matches
+                               int numImages,          // Number of images
+                               float tol, // Relative tolerance for convergence
+                               int iter1, // Number of iterations for pass 1
+                               float rej1, // Limit on rejection between iterations for pass 1
+                               int iter2, // Number of iterations for pass 2
+                               float rej2, // Limit on rejection between iterations for pass 2
+                               float rejLimit // Limit on rejection between iterations
+    );
+
 #endif