Changeset 41841

Timestamp:

Oct 15, 2021, 6:45:21 AM (5 years ago)

Author:

eugene

Message:

adding simple percentile stack

Location:

branches/eam_branches/ppStack.20211015/src

Files:

• ppStackCombineInitial.c (modified) (1 diff)
• ppStackCombinePercent.c (added)
• ppStackCombinePrepare.c (modified) (4 diffs)
• ppStackLoop.c (modified) (5 diffs)
• ppStackLoop.h (modified) (3 diffs)
• ppStackLoopByPercent.c (added)
• ppStackMatch.c (modified) (1 diff)
• ppStackOptions.c (modified) (1 diff)
• ppStackOptions.h (modified) (1 diff)
• ppStackReadout.c (modified) (2 diffs)
• ppStackSetup.c (modified) (1 diff)
• ppStackThread.c (modified) (1 diff)

branches/eam_branches/ppStack.20211015/src/ppStackCombineInitial.c

@@ -9 +9 @@
     psAssert(options, "Require options");
     psAssert(config, "Require configuration");
-
-/*     if (!options->convolve) { */
-/*         // No need to do initial combination when we haven't convolved */
-/*         // XXX either allocate inspect and rejected here, or do not require them downstream */
-/*         return true; */
-/*     } */
 
     psTimerStart("PPSTACK_INITIAL");

branches/eam_branches/ppStack.20211015/src/ppStackCombinePrepare.c

@@ -31 +31 @@
       options->expRO = pmReadoutAlloc(expCell); //Output readout
     }
-/*     else { */
-/*       options->expRO = NULL; */
-/*     } */
+
     int bkg_r0 = 0,bkg_c0 = 0;
     int bkg_nC = 0,bkg_nR = 0;
@@ -53 +51 @@
       
       options->bkgRO = pmReadoutAlloc(bkgCell); // BKG readout
-      //      if (!pmHDUGenerateForFPA(options->bkgRO->parent->parent->parent)) {
       options->bkgRO->parent->parent->parent->hdu = pmHDUAlloc(NULL);
+
       if (!options->bkgRO->parent->parent->parent->hdu) {
         fprintf(stderr,"failed to generate a HDU for this thing.\n");
@@ -68 +66 @@
                                                                         ro->parent->parent->parent->concepts);
 
-    }
-    else {
+    } else {
       options->bkgRO = NULL;
     }
@@ -97 +94 @@
         return false;
       }
-      
-
-      
     }
 

branches/eam_branches/ppStack.20211015/src/ppStackLoop.c

@@ -31 +31 @@
     ppStackMemDump("prepare");
     if (options->quality) return true; // Can't do anything else
+
+    // if we want to skip the rejection steps, and instead use a percentile clipping:
+    if (options->clipPercent) {
+      bool result = ppStackLoopByPercent (config, options);
+      return result;
+    }
 
     // Convolve inputs
@@ -88 +94 @@
         pmCellFreeData(options->cells->data[i]);
     }
-    //MEH -- must uncomment back out -- unclear is should be moved after pixel rejection
+    // MEH -- must uncomment back out -- unclear is should be moved after pixel rejection
     psFree(stack);
 
@@ -126 +132 @@
         return false;
       }
-    }
-    else {
+    } else {
       // Since we haven't convolved, I believe we do need to normalize here.
-      //MEH -- see below for comment on ppStackCombineFinal and bscaleOffset
+      // MEH -- see below for comment on ppStackCombineFinal and bscaleOffset
       if (!ppStackCombineFinal(stack, options->convCovars, options, config, false, true, true, true)) {
         psError(psErrorCodeLast(), false, "Unable to perform final combination.");
@@ -154 +159 @@
       return false;
     }
-    ppStackFileActivation(config,PPSTACK_FILES_BKG ,false);    
+    ppStackFileActivation(config, PPSTACK_FILES_BKG, false);    
+
     // Photometry
     psTrace("ppStack", 1, "Photometering stacked image....\n");
@@ -184 +190 @@
         return false;
     }
-    //MEH -- also must uncomment back out..
+    // MEH -- also must uncomment back out..
     psFree(stack);
     psLogMsg("ppStack", PS_LOG_INFO, "Stage 7: Cleanup, WCS & JPEGS: %f sec", psTimerClear("PPSTACK_STEPS"));

branches/eam_branches/ppStack.20211015/src/ppStackLoop.h

@@ -7 +7 @@
     ppStackOptions *options             // Options for stacking
     );
+
 // Median only loop.
 bool ppStackMedianLoop(
@@ -12 +13 @@
     ppStackOptions *options             // Options for stacking
     );
+
+// stack by percentile range
+bool ppStackLoopByPercent(
+  pmConfig *config,
+  ppStackOptions *options
+  );
 
 // Setup
@@ -92 +99 @@
 
 
+bool ppStackCombinePercent(ppStackThreadData *stack, ppStackOptions *options, pmConfig *config);
+
 bool ppStackUpdateHeader(ppStackThreadData *stack, ppStackOptions *options, pmConfig *config);
 bool ppStackJPEGs(ppStackThreadData *stack, ppStackOptions *options, pmConfig *config);

branches/eam_branches/ppStack.20211015/src/ppStackMatch.c

@@ -524 +524 @@
         psFree(conv);
     } else {
+        // this branch is done if options->convolve is FALSE
         // Match the normalisation
         float norm = powf(10.0, -0.4 * options->norm->data.F32[index]); // Normalisation

branches/eam_branches/ppStack.20211015/src/ppStackOptions.c

@@ -56 +56 @@
     options->num = 0;
     options->quality = 0;
+
+    options->clipPercent = false;
+    options->clipMinRange = 0.16; // 1-sigma default
+    options->clipMaxRange = 0.84; // 1-sigma default
+
     options->psf = NULL;
     options->sumExposure = NAN;

branches/eam_branches/ppStack.20211015/src/ppStackOptions.h

@@ -16 +16 @@
     int num;                            // Number of inputs
     int quality;                        // Bad data quality flag
+
+    bool clipPercent;                   // use percentile range to clip?
+    psF32 clipMinRange;                 // min percentile
+    psF32 clipMaxRange;                 // max percentile
+
     // Prepare
     pmPSF *psf;                         // Target PSF

branches/eam_branches/ppStack.20211015/src/ppStackReadout.c

@@ -53 +53 @@
 }
 
+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 +339 @@
     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;
+}
+
+

branches/eam_branches/ppStack.20211015/src/ppStackSetup.c

@@ -18 +18 @@
     options->convolve = psMetadataLookupBool(NULL, recipe, "CONVOLVE"); // Convolve images?
     options->doBackground = psMetadataLookupBool(NULL, recipe, "BACKGROUND.MODEL");
+
+    options->clipPercent  = psMetadataLookupBool(NULL, recipe, "COMBINE.PERCENT"); // use percentile range to clip?
+    options->clipMinRange = psMetadataLookupF32 (NULL, recipe, "COMBINE.MIN.RANGE"); // min percentile
+    options->clipMaxRange = psMetadataLookupF32 (NULL, recipe, "COMBINE.MAX.RANGE"); // max percentile
+
 
     if (!psMetadataLookupBool(NULL, config->arguments, "HAVE.PSF")) {

branches/eam_branches/ppStack.20211015/src/ppStackThread.c

@@ -283 +283 @@
     }
 
+    {
+        psThreadTask *task = psThreadTaskAlloc("PPSTACK_PERCENT_COMBINE", 3);
+        task->function = &ppStackReadoutPercentThread;
+        psThreadTaskAdd(task);
+        psFree(task);
+    }
     return;
 }