Changeset 37067 for branches/eam_branches/ipp-ops-20130712/ppStack/src

Timestamp:

Jul 17, 2014, 12:32:26 PM (12 years ago)

Author:

eugene

Message:

merge changes (from past YEAR) into trunk

Location:

branches/eam_branches/ipp-ops-20130712/ppStack/src

Files:

• . (modified) (1 prop)
• ppStack.h (modified) (1 diff)
• ppStackCombineFinal.c (modified) (4 diffs)
• ppStackLoop.c (modified) (5 diffs)
• ppStackLoop.h (modified) (1 diff)
• ppStackMatch.c (modified) (2 diffs)
• ppStackOptions.c (modified) (2 diffs)
• ppStackOptions.h (modified) (1 diff)
• ppStackPrepare.c (modified) (2 diffs)
• ppStackReadout.c (modified) (10 diffs)
• ppStackThread.c (modified) (1 diff)

branches/eam_branches/ipp-ops-20130712/ppStack/src/ppStack.h

@@ -108 +108 @@
                          const psVector *addVariance, // Additional variance for rejection
                          bool safety,                 // Enable safety switch?
-                         const psVector *norm         // Normalisations to apply
+                         const psVector *norm,         // Normalisations to apply
+                         const psVector *bscaleApplyOffset  // hack for offset based on bscale
     );
 

branches/eam_branches/ipp-ops-20130712/ppStack/src/ppStackCombineFinal.c

@@ -4 +4 @@
 // #define TESTING                         // Enable test output
 
+//MEH -- adhoc addition to blank mask border of final stack since rejection different/none on order of KERNEL.SIZE with overlap in CombineInitial 
+static void stackBorderMask(psImage *image, // Image to mark as blank
+                            psImage *mask, // Mask to mark as blank (or NULL)
+                            psImage *variance, // Weight map to mark as blank (or NULL)
+                            int numCols, int numRows, // Size of image
+                            int size, // Size to mark blank
+                            psImageMaskType blank // Blank mask value
+    )
+{
+    for (int y = size; y < numRows - size; y++) {
+        for (int x = 0; x < size; x++) {
+            image->data.F32[y][x] = NAN;
+            mask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] = blank;
+            variance->data.F32[y][x] = NAN;
+        }
+        for (int x = numCols - size; x < numCols; x++) {
+            image->data.F32[y][x] = NAN;
+            mask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] = blank;
+            variance->data.F32[y][x] = NAN;
+        }
+    }
+    for (int y = 0; y < size; y++) {
+        for (int x = 0; x < numCols; x++) {
+            image->data.F32[y][x] = NAN;
+            mask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] = blank;
+            variance->data.F32[y][x] = NAN;
+        }
+    }
+    for (int y = numRows - size; y < numRows; y++) {
+        for (int x = 0; x < numCols; x++) {
+            image->data.F32[y][x] = NAN;
+            mask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] = blank;
+            variance->data.F32[y][x] = NAN;
+        }
+    }
+    return;
+}
+
+
 bool ppStackCombineFinal(ppStackThreadData *stack, psArray *covariances, ppStackOptions *options,
-                         pmConfig *config, bool safe, bool normalise, bool grow)
+                         pmConfig *config, bool safe, bool normalise, bool grow, bool bscaleoffset)
 {
     psAssert(stack, "Require stack");
@@ -20 +59 @@
     psAssert(recipe, "We've thrown an error on this before.");
     float poorFrac = psMetadataLookupF32(NULL, recipe, "POOR.FRACTION"); // Fraction for "poor"
+
+    int sizeBlank = psMetadataLookupS32(NULL, recipe, "MASK.BLANKBORDER"); // Pixels to mask BLANK from edge 
+    psImageMaskType maskBlank = pmConfigMaskGet("BLANK", config); // Bits to mask for bad pixels
 
     // Grow the list of rejected pixels, if desired
@@ -65 +107 @@
         PS_ARRAY_ADD_SCALAR(job->args, safe, PS_TYPE_U8);
         PS_ARRAY_ADD_SCALAR(job->args, normalise, PS_TYPE_U8);
+        PS_ARRAY_ADD_SCALAR(job->args, bscaleoffset, PS_TYPE_U8);
         if (!psThreadJobAddPending(job)) {
             psFree(reject);
@@ -109 +152 @@
 #endif
 
+    //MEH blank mask/manual reject border on final stack -- 
+    if (sizeBlank > 0) {
+        stackBorderMask(outRO->image,outRO->mask,outRO->variance,numCols,numRows,sizeBlank,maskBlank);
+        stackBorderMask(expRO->image,expRO->mask,expRO->variance,numCols,numRows,sizeBlank,0);
+    }
+
     if (options->stats) {
         // only add the timer if it has not been set (convolved stack)

branches/eam_branches/ipp-ops-20130712/ppStack/src/ppStackLoop.c

@@ -87 +87 @@
         pmCellFreeData(options->cells->data[i]);
     }
-    //    psFree(stack);
+    //MEH -- must uncomment back out -- unclear is should be moved after pixel rejection
+    psFree(stack);
 
     // Pixel rejection
@@ -119 +120 @@
     psTrace("ppStack", 2, "Final stack of convolved images....\n");
     if (options->convolve) {
-      if (!ppStackCombineFinal(stack, options->convCovars, options, config, false, false, true)) {
+      if (!ppStackCombineFinal(stack, options->convCovars, options, config, false, false, true, false)) {
         psError(psErrorCodeLast(), false, "Unable to perform final combination.");
         psFree(stack);
@@ -127 +128 @@
     else {
       // Since we haven't convolved, I believe we do need to normalize here.
-      if (!ppStackCombineFinal(stack, options->convCovars, options, config, false, true, true)) {
+      //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.");
         psFree(stack);
@@ -181 +183 @@
         return false;
     }
-    //    psFree(stack);
+    //MEH -- also must uncomment back out..
+    psFree(stack);
     psLogMsg("ppStack", PS_LOG_INFO, "Stage 7: Cleanup, WCS & JPEGS: %f sec", psTimerClear("PPSTACK_STEPS"));
     ppStackMemDump("cleanup");
@@ -205 +208 @@
         // generate the unconvolved stack. NOTE: this one must be normalized since the inputs have not been
         psTrace("ppStack", 2, "Stack of unconvolved images....\n");
-        if (!ppStackCombineFinal(stack, options->origCovars, options, config, false, true, false)) {
+        //MEH -- terrible hack for bscale offset to input warps -- treat like normalization 
+        // -- if BSCALEOFFSET TRUE, 0.5*bscale for offset values, otherwise value 0.0
+        // -- bscaleOffset set when files read in and applied in CombineFinal (if arg true, only valid for unconv case)
+        if (!ppStackCombineFinal(stack, options->origCovars, options, config, false, true, false, true)) {
             psError(psErrorCodeLast(), false, "Unable to perform unconvolved combination.");
             psFree(stack);

branches/eam_branches/ipp-ops-20130712/ppStack/src/ppStackLoop.h

@@ -65 +65 @@
     bool safe,                          // Allow safe combination?
     bool norm,                          // Normalise images?
-    bool grow                           // Grow rejection masks?
+    bool grow,                           // Grow rejection masks?
+    bool bscaleoffset                   // Apply bscale offset?
     );
 

branches/eam_branches/ipp-ops-20130712/ppStack/src/ppStackMatch.c

@@ -6 +6 @@
 #define COVAR_FRAC 0.01                 // Truncation fraction for covariance matrix
 
-// #define TESTING                         // Enable debugging output
+// #define TESTING                      // Enable debugging output
 // #define TESTING_REUSE_CONV           // Enable debugging for re-using convolved outputs from previous run
-
 // #define TESTING_CZW 
 
@@ -156 +155 @@
     if (options->convolve) {
         pmReadout *conv = pmReadoutAlloc(NULL); // Conv readout, for holding results temporarily
+
+        // MEH -- earlier defn needed with new simple mode modifications to  use TESTING_REUSE_CONV...
+        //const char *typeStr = psMetadataLookupStr(NULL, ppsub, "KERNEL.TYPE"); // Kernel type
+        //pmSubtractionKernelsType type = pmSubtractionKernelsTypeFromString(typeStr); // Kernel type
+
 #ifdef TESTING_REUSE_CONV
         // This is a hack to use the temporary convolved images and kernel generated previously.

branches/eam_branches/ipp-ops-20130712/ppStack/src/ppStackOptions.c

@@ -19 +19 @@
     psFree(options->inputMask);
     psFree(options->sourceLists);
+    psFree(options->bscaleOffset);
     psFree(options->norm);
     psFree(options->sources);
@@ -65 +66 @@
     options->inputMask = NULL;
     options->sourceLists = NULL;
+    options->bscaleOffset = NULL;
     options->norm = NULL;
     options->sources = NULL;

branches/eam_branches/ipp-ops-20130712/ppStack/src/ppStackOptions.h

@@ -35 +35 @@
     float clippedMean;                  // clipped mean of input fwhm
     float clippedStdev;                 // clipped stdev of input fwhm 
+    psVector *bscaleOffset;             // hack offset of bzero using 0.5*bscale if option, 0.0 if not
     // Convolve
     psArray *cells;                     // Cells for convolved images --- a handle for reading again

branches/eam_branches/ipp-ops-20130712/ppStack/src/ppStackPrepare.c

@@ -148 +148 @@
     }
 
+    //MEH -- bscale offset hack for warps 
+    options->bscaleOffset = psVectorAlloc(options->num, PS_TYPE_F32);
+    psVectorInit(options->bscaleOffset, 0.0);
+    bool setbscaleoffset = psMetadataLookupBool(NULL,recipe, "BSCALEOFFSET");
+
     psArray *psfs = psArrayAlloc(num); // PSFs for PSF envelope
     int numCols = 0, numRows = 0;   // Size of image
@@ -162 +167 @@
         }
         options->sumExposure += options->exposures->data.F32[i];
+
+        //MEH -- hdu get redefn again below if conv (lots of this happening..) 
+        if (setbscaleoffset) {
+            pmHDU *hdu = pmHDUFromCell(cell);
+            assert(hdu && hdu->header);
+            float bscale = psMetadataLookupF32(NULL, hdu->header, "BSCALE");
+            options->bscaleOffset->data.F32[i] = 0.5*bscale;
+        }
+        psLogMsg("ppStack", PS_LOG_INFO, "bscaleOffset: %d %f", i,options->bscaleOffset->data.F32[i]);
 
         // Get list of PSFs, to determine target PSF

branches/eam_branches/ipp-ops-20130712/ppStack/src/ppStackReadout.c

@@ -34 +34 @@
     bool safety = PS_SCALAR_VALUE(args->data[4], U8);    // Safety switch on?
     bool normalise = PS_SCALAR_VALUE(args->data[5], U8); // Normalise images?
+    bool bscaleoffset = PS_SCALAR_VALUE(args->data[6], U8); // Apply bscale offset?
 
     psVector *mask = options->inputMask; // Mask for inputs
@@ -40 +41 @@
     psVector *addVariance = options->matchChi2; // Additional variance when rejecting
     psVector *norm = normalise ? options->norm : NULL; // Normalisations to apply to images
+    psVector *bscaleApplyOffset = bscaleoffset ? options->bscaleOffset : NULL; // BSCALE offset to apply to images
 
     bool status = ppStackReadoutFinal(config, options->outRO, options->expRO, thread->readouts, mask, reject,
-                                      weightings, exposures, addVariance, safety, norm); // Status of operation
+                                      weightings, exposures, addVariance, safety, norm, bscaleApplyOffset); // Status of operation
 
     thread->busy = false;
@@ -116 +118 @@
     bool useVariance = psMetadataLookupBool(&mdok, recipe, "VARIANCE"); // Use variance for rejection?
     bool safe = psMetadataLookupBool(&mdok, recipe, "SAFE"); // Be safe when combining small numbers of pixels
+   
+    int nminpix = psMetadataLookupS32(&mdok, recipe, "NMINPIX"); // Minimum input per pixel to combine with
 
     psMetadata *ppsub = psMetadataLookupMetadata(NULL, config->recipes, "PPSUB"); // PPSUB recipe
@@ -161 +165 @@
 
     if (!pmStackCombine(outRO, NULL, stack, maskBad, maskSuspect, maskBlank, kernelSize, iter,
-                        combineRej, combineSys, combineDiscard, useVariance, safe, false)) {
+                        combineRej, combineSys, combineDiscard, useVariance, safe, nminpix, false)) {
         psError(PS_ERR_UNKNOWN, false, "Unable to combine input readouts with rejection.");
         psFree(stack);
@@ -184 +188 @@
     psFree(stack);
 
+    //MEH change to trace
+    //psLogMsg("ppStack", PS_LOG_INFO, "initial stack image sectionNum %d", sectionNum);
+
     sectionNum++;
 
@@ -197 +204 @@
                          const psVector *mask, const psArray *rejected, const psVector *weightings,
                          const psVector *exposures, const psVector *addVariance, bool safety,
-                         const psVector *norm)
+                         const psVector *norm, const psVector *bscaleApplyOffset)
 {
     assert(config);
@@ -216 +223 @@
     bool useVariance = psMetadataLookupBool(&mdok, recipe, "VARIANCE"); // Use variance for rejection?
     bool safe = psMetadataLookupBool(&mdok, recipe, "SAFE"); // Be safe when combining small numbers of pixels
+
+    int nminpix = psMetadataLookupS32(&mdok, recipe, "NMINPIX"); // Minimum input per pixel to combine with
 
     psString maskBadStr = psMetadataLookupStr(NULL, recipe, "MASK.VAL"); // Name of bits to mask for bad
@@ -265 +274 @@
         stack->data[i] = data;
 
+        //MEH -- apply bscale offset before norm   
+        if (bscaleApplyOffset) {
+            //MEH change to trace
+            //psLogMsg("ppStack", PS_LOG_INFO, "bscaleApplyOffset: %d %f", i, bscaleApplyOffset->data.F32[i]);
+            psBinaryOp(ro->image, ro->image, "-", psScalarAlloc(bscaleApplyOffset->data.F32[i], PS_TYPE_F32));
+        }
+
         if (norm) {
             float normalise = powf(10.0, -0.4 * norm->data.F32[i]); // Normalisation
@@ -273 +289 @@
 
     if (!pmStackCombine(outRO, expRO, stack, maskBad, maskSuspect, maskBlank, 0, iter, combineRej,
-                        combineSys, combineDiscard, useVariance, safe, rejected)) {
+                        combineSys, combineDiscard, useVariance, safe, nminpix, rejected)) {
         psError(PS_ERR_UNKNOWN, false, "Unable to combine input readouts.");
         psFree(stack);
@@ -293 +309 @@
     psFree(stack);
 
+    //MEH change to trace
+    //psLogMsg("ppStack", PS_LOG_INFO, "final stack image sectionNum %d", sectionNum);
+
     sectionNum++;
 

branches/eam_branches/ipp-ops-20130712/ppStack/src/ppStackThread.c

@@ -277 +277 @@
 
     {
-        psThreadTask *task = psThreadTaskAlloc("PPSTACK_FINAL_COMBINE", 6);
+        psThreadTask *task = psThreadTaskAlloc("PPSTACK_FINAL_COMBINE", 7);
         task->function = &ppStackReadoutFinalThread;
         psThreadTaskAdd(task);