Changeset 42094 for trunk/ppStack/src/ppStackReadout.c

Timestamp:

Feb 28, 2022, 3:01:31 PM (4 years ago)

Author:

eugene

Message:

some cleanup of ppStackFiles: JPEG1 & JPEG2 are not part of the "convolved" output: they are activated explicitly in ppStackJPEGs.c; add option to generate stacks using the pmStackCombineByPercent method; trap failures in the ppStackReadout threads; fix normalization for both convolved and unconvolved versions (was inconsistent in ppStackLoop)

Location:

trunk/ppStack

Files:

• . (modified) (1 prop)
• src/ppStackReadout.c (modified) (12 diffs)

trunk/ppStack/src/ppStackReadout.c

@@ -19 +19 @@
                                          weightings, exposures, addVariance);
     thread->busy = false;
+
+    thread->status = job->results ? PPSTACK_THREAD_SUCCESS : PPSTACK_THREAD_FAILURE;
+
+    psAssert(job->results, "Stacking failed.");
 
     return job->results ? true : false;
@@ -48 +52 @@
     thread->busy = false;
 
+    thread->status = status ? PPSTACK_THREAD_SUCCESS : PPSTACK_THREAD_FAILURE;
+
     psAssert(status, "Stacking failed.");
 
@@ -53 +59 @@
 }
 
+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
+    pmReadout *expRO = options->expRO;  // Expmap 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
+
+    // fprintf (stderr, "start thread %d : scan range %d - %d\n", thread->entry, thread->firstScan, thread->lastScan);
+
+    bool status = ppStackReadoutPercent(config, outRO, expRO, thread->readouts, mask, weightings, exposures, addVariance);
+
+    thread->busy = false;
+
+    thread->status = status ? PPSTACK_THREAD_SUCCESS : PPSTACK_THREAD_FAILURE;
+
+    // fprintf (stderr, "finished thread %d, status %d : scan range %d - %d\n", thread->entry, thread->status, thread->firstScan, thread->lastScan);
+
+    psAssert(status, "Stacking failed.");
+
+    return status;
+}
+
+//////////////////////////////////// 
 
 bool ppStackInspect(psThreadJob *job)
@@ -145 +183 @@
 
     int num = readouts->n;              // Number of inputs
-    psArray *stack = psArrayAlloc(num); // Array for stacking
+    psArray *stackData = psArrayAlloc(num); // Array of data to be stacked
 
     for (int i = 0; i < num; i++) {
@@ -160 +198 @@
         }
 
-        stack->data[i] = pmStackDataAlloc(ro, weightings->data.F32[i], exposures->data.F32[i],
+        // stackData is an array of pmStackData structures
+        stackData->data[i] = pmStackDataAlloc(ro, weightings->data.F32[i], exposures->data.F32[i],
                                           addVariance->data.F32[i]);
     }
 
-    if (!pmStackCombine(outRO, NULL, stack, maskBad, maskSuspect, maskBlank, kernelSize, iter,
+    if (!pmStackCombine(outRO, NULL, stackData, maskBad, maskSuspect, maskBlank, kernelSize, iter,
                         combineRej, combineSys, combineDiscard, useVariance, safe, nminpix, false)) {
         psError(PS_ERR_UNKNOWN, false, "Unable to combine input readouts with rejection.");
-        psFree(stack);
+        psFree(stackData);
+        
+        // XXX EAM : if pmStackCombine fails above this return will cause the thread to block.
+        // The failure status results in a psThread.fault.  the psThreadLauncher function, which
+        // each background thread runs as a loop, will sleep until the fault is cleared externally
+        // by the handler.  psThreadPoolWait() catches, counts, and clears these faults, but
+        // ppStackThreadRead does not.  
         return false;
     }
@@ -175 +220 @@
     psArray *reject = psArrayAlloc(num);  // List of pixels rejected
     for (int i = 0; i < num; i++) {
-        pmStackData *data = stack->data[i]; // Data for this image
+        pmStackData *data = stackData->data[i]; // Data for this image
         if (!data) {
             continue;
@@ -186 +231 @@
         reject->data[i] = psMemIncrRefCounter(data->reject);
     }
-    psFree(stack);
+    psFree(stackData);
 
     //MEH change to trace
@@ -247 +292 @@
 
     int num = readouts->n;              // Number of inputs
-    psArray *stack = psArrayAlloc(num); // Array for stacking
+    psArray *stackData = psArrayAlloc(num); // Array of data to be stacked
 
     // We have rejection from a previous combination: combine without flagging pixels to inspect
@@ -272 +317 @@
                                              addVariance ? addVariance->data.F32[i] : NAN);
         data->reject = rejected ? psMemIncrRefCounter(rejected->data[i]) : NULL;
-        stack->data[i] = data;
+        stackData->data[i] = data;
 
         //MEH -- apply bscale offset before norm   
@@ -288 +333 @@
     }
 
-    if (!pmStackCombine(outRO, expRO, stack, maskBad, maskSuspect, maskBlank, 0, iter, combineRej,
+    if (!pmStackCombine(outRO, expRO, stackData, maskBad, maskSuspect, maskBlank, 0, iter, combineRej,
                         combineSys, combineDiscard, useVariance, safe, nminpix, rejected)) {
         psError(PS_ERR_UNKNOWN, false, "Unable to combine input readouts.");
-        psFree(stack);
+        psFree(stackData);
         return false;
     }
@@ -307 +352 @@
     expChip->data_exists = true;
 
-    psFree(stack);
+    psFree(stackData);
 
     //MEH change to trace
@@ -316 +361 @@
     return true;
 }
+
+// NOTE: expRO is allowed to be NULL
+bool ppStackReadoutPercent(const pmConfig *config, pmReadout *outRO, pmReadout *expRO, 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
+
+    // fprintf (stderr, "starting ReadoutPercent, %d\n", sectionNum);
+
+    // 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
+    psF64 rejectFraction = psMetadataLookupF32(&mdok, recipe, "COMBINE.REJECT.FRACTION"); // fraction of outliers to reject
+
+    char defaultBlankStr[16] = "BLANK";
+
+    psString maskBadStr     = psMetadataLookupStr(&mdok, recipe, "MASK.VAL"); // Name of bits to mask for bad
+    psString maskSuspectStr = psMetadataLookupStr(&mdok, recipe, "MASK.SUSPECT"); // Name of suspect mask bits
+    psString maskBlankStr   = psMetadataLookupStr(&mdok, recipe, "MASK.BLANK"); // Name of bits to set for empty pixels
+    if (!maskBlankStr) {
+      maskBlankStr = psMetadataLookupStr(&mdok, recipe, "MASK.BAD"); // Old name for MASK.BLANK
+    }
+    if (!maskBlankStr) {
+      maskBlankStr = defaultBlankStr; // this is statically allocated above
+    }
+
+    psImageMaskType maskBad     = pmConfigMaskGet(maskBadStr, config); // Bits to mask for bad pixels
+    psImageMaskType maskSuspect = pmConfigMaskGet(maskSuspectStr, config); // Suspect bits
+    psImageMaskType maskBlank   = pmConfigMaskGet(maskBlankStr, config); // Bits to mask for bad pixels
+
+    int num = readouts->n;              // Number of inputs
+    psArray *stackData = psArrayAlloc(num); // Array of data to be stacked
+
+    for (int i = 0; i < num; i++) {
+        stackData->data[i] = NULL;
+
+        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);
+        }
+
+        stackData->data[i] = pmStackDataAlloc(ro, weightings->data.F32[i], exposures->data.F32[i],
+                                              addVariance->data.F32[i]);
+    }
+
+    if (!pmStackCombineByPercentile(outRO, expRO, stackData, rejectFraction, maskBad, maskSuspect, maskBlank)) {
+        psError(PS_ERR_UNKNOWN, false, "Unable to combine input readouts with rejection.");
+        psFree(stackData);
+        return false;
+    }
+
+    outRO->data_exists = true;                 // output readout
+    outRO->parent->data_exists = true;         // output cell
+    outRO->parent->parent->data_exists = true; // output chip
+
+    if (expRO) {
+      expRO->data_exists = true;                 // expmap readout 
+      expRO->parent->data_exists = true;               // expmap cell   
+      expRO->parent->parent->data_exists = true; // expmap chip 
+    }
+      
+    psFree(stackData);
+    sectionNum++;
+
+    return true;
+}
+
+