Index: branches/eam_branches/ppStack.20211015/src/ppStackReadout.c
===================================================================
--- branches/eam_branches/ppStack.20211015/src/ppStackReadout.c	(revision 41840)
+++ branches/eam_branches/ppStack.20211015/src/ppStackReadout.c	(revision 41841)
@@ -53,4 +53,27 @@
 }
 
+bool ppStackReadoutPercentThread(psThreadJob *job)
+{
+    PS_ASSERT_THREAD_JOB_NON_NULL(job, false);
+
+    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
+    psVector *mask = options->inputMask; // Mask for inputs
+    psVector *weightings = options->weightings; // Weightings (1/noise^2) for each image
+    psVector *exposures = options->exposures;   // Exposure times for each image
+    psVector *addVariance = options->matchChi2; // Additional variance when rejecting
+
+    job->results = ppStackReadoutInitial(config, outRO, thread->readouts, mask,
+                                         weightings, exposures, addVariance);
+    thread->busy = false;
+
+    return job->results ? true : false;
+}
+
+//////////////////////////////////// 
 
 bool ppStackInspect(psThreadJob *job)
@@ -316,2 +339,83 @@
     return true;
 }
+
+psArray *ppStackReadoutPercent(const pmConfig *config, pmReadout *outRO, const psArray *readouts,
+			       const psVector *mask, const psVector *weightings, const psVector *exposures,
+			       const psVector *addVariance)
+{
+    assert(config);
+    assert(outRO);
+    assert(readouts);
+    assert(mask && mask->n == readouts->n && mask->type.type == PS_TYPE_VECTOR_MASK);
+    assert(weightings && weightings->n == readouts->n && weightings->type.type == PS_TYPE_F32);
+    assert(addVariance && addVariance->n == readouts->n && addVariance->type.type == PS_TYPE_F32);
+    static int sectionNum = 0;          // Section number; for debugging outputs
+
+    // Get the recipe values
+    psMetadata *recipe = psMetadataLookupMetadata(NULL, config->recipes, PPSTACK_RECIPE); // ppStack recipe
+    psAssert(recipe, "We've thrown an error on this before.");
+
+    bool mdok;                          // Status of MD lookup
+    float iter = psMetadataLookupF32(NULL, recipe, "COMBINE.ITER"); // Rejection iterations
+    float minRange = psMetadataLookupF32(NULL, recipe, "COMBINE.MIN.RANGE"); // Combination threshold
+    float maxRange = psMetadataLookupF32(NULL, recipe, "COMBINE.MAX.RANGE"); // Combination threshold
+
+    psMetadata *ppsub = psMetadataLookupMetadata(NULL, config->recipes, "PPSUB"); // PPSUB recipe
+    int kernelSize = psMetadataLookupS32(NULL, ppsub, "KERNEL.SIZE"); // Kernel half-size
+
+    psString maskBadStr = psMetadataLookupStr(NULL, recipe, "MASK.VAL"); // Name of bits to mask for bad
+    psImageMaskType maskBad = pmConfigMaskGet(maskBadStr, config); // Bits to mask for bad pixels
+
+    psString maskSuspectStr = psMetadataLookupStr(NULL, recipe, "MASK.SUSPECT"); // Name of suspect mask bits
+    psImageMaskType maskSuspect = pmConfigMaskGet(maskSuspectStr, config); // Suspect bits
+
+    bool status = false;
+    psImageMaskType maskBlank;
+    psString maskBlankStr = psMetadataLookupStr(&status, recipe, "MASK.BLANK"); // Name of bits to set for empty pixels
+    if (maskBlankStr) {
+      maskBlank = pmConfigMaskGet(maskBlankStr, config); // Bits to mask for bad pixels
+    } else {
+      maskBlankStr = psMetadataLookupStr(&status, recipe, "MASK.BAD"); // Old name for MASK.BLANK
+      if (maskBlankStr) {
+	maskBlank = pmConfigMaskGet(maskBlankStr, config); // Bits to mask for bad pixels
+      } else {
+	maskBlank = pmConfigMaskGet("BLANK", config); // Bits to mask for bad pixels
+      }
+    }
+
+    int num = readouts->n;              // Number of inputs
+    psArray *stack = psArrayAlloc(num); // Array for stacking
+
+    for (int i = 0; i < num; i++) {
+        pmReadout *ro = readouts->data[i];
+        if (!ro || mask->data.PS_TYPE_VECTOR_MASK_DATA[i]) {
+            // Bad image
+            continue;
+        }
+
+        // Ensure there is a mask, or pmStackCombine will complain
+        if (!ro->mask) {
+            ro->mask = psImageAlloc(ro->image->numCols, ro->image->numRows, PS_TYPE_IMAGE_MASK);
+            psImageInit(ro->mask, 0);
+        }
+
+        stack->data[i] = pmStackDataAlloc(ro, weightings->data.F32[i], exposures->data.F32[i],
+                                          addVariance->data.F32[i]);
+    }
+
+    if (!pmStackCombineByPercentile(outRO, stack, minRange, maxRange, maskBad, maskSuspect, maskBlank)) {
+        psError(PS_ERR_UNKNOWN, false, "Unable to combine input readouts with rejection.");
+        psFree(stack);
+        return false;
+    }
+
+    sectionNum++;
+
+    psArray *results = psArrayAlloc(2); // Array of results
+    results->data[0] = inspect;
+    results->data[1] = reject;
+
+    return results;
+}
+
+