Changeset 36855

Timestamp:

Jun 12, 2014, 11:52:06 PM (12 years ago)

Author:

mhuber

Message:

adding BSCALEOFFSET option hack to apply -0.5*bscale offset to each input warp, default is -0.0 offset

Location:

trunk

Files:

• ippconfig/gpc1/ppStack.config (modified) (3 diffs)
• ippconfig/recipes/ppStack.config (modified) (1 diff)
• ppStack/src/ppStack.h (modified) (1 diff)
• ppStack/src/ppStackCombineFinal.c (modified) (2 diffs)
• ppStack/src/ppStackLoop.c (modified) (3 diffs)
• ppStack/src/ppStackLoop.h (modified) (1 diff)
• ppStack/src/ppStackOptions.c (modified) (2 diffs)
• ppStack/src/ppStackOptions.h (modified) (1 diff)
• ppStack/src/ppStackPrepare.c (modified) (2 diffs)
• ppStack/src/ppStackReadout.c (modified) (4 diffs)
• ppStack/src/ppStackThread.c (modified) (1 diff)

trunk/ippconfig/gpc1/ppStack.config

@@ -160 +160 @@
     STACK.TYPE              STR   DEEP_STACK  ## needed?
     PSF.INPUT.CLIP.SIMPLE   BOOL  True
-    PSF.INPUT.MAX           F32   6.0   # never really want >6 pixels in refstack, even on edges 
+    PSF.INPUT.MAX           F32   7.0   # center/edge can vary 1-2 pix, prefer <6 pix but needs to be <7 pix for edges
     PSF.INPUT.CLIP.NSIGMA   F32   1.0   # sample typically on rising side of distribution (set even smaller?)
     PSF.INPUT.THRESH        F32   NAN
@@ -168 +168 @@
     PSF.TARGET.AS.MAX       BOOL  TRUE      # Set the target PSF FWHM as the maximum of accepted input FWHM values.
     PSF.TARGET.AS.MAX.EPSILON F32 0.1       # Amount to set the target PSF FWHM larger than the maximum input. Target = eps + max(input)
+    THRESHOLD.MASK  F32     0.1             # Threshold for mask deconvolution (0..1)
+    COMBINE.ITER    F32     0.5             # Number of rejection iterations per input
+    COMBINE.REJ     F32     3.0             # Rejection threshold in combination (sigma)
+    COMBINE.SYS     F32     0.1             # Relative systematic error in combination
+    COMBINE.DISCARD F32     0.2             # Discard fraction for Olympic weighted mean
+    IMAGE.REJ       F32     0.1             # Rejected pixel fraction threshold for rejecting entire image
+    NMINPIX         S32     4               # Minimum input per pixel
+    MASK.VAL        STR     SUSPECT,MASK.VALUE,CONV.BAD,GHOST # Mask value of input bad pixels
+    MASK.BLANKBORDER S32    20              # Mask blank border in final stack output
+    BSCALEOFFSET    BOOL    TRUE            # HACK for removing bscale offset from input warps because of compression (0.5*BSCALE)
 END
 
@@ -193 +203 @@
     MASK.VAL        STR     SUSPECT,MASK.VALUE,CONV.BAD,GHOST # Mask value of input bad pixels
     MASK.BLANKBORDER S32    20              # Mask blank border in final stack output
+    BSCALEOFFSET    BOOL    TRUE            # HACK for removing bscale offset from input warps because of compression (0.5*BSCALE)
 END
 

trunk/ippconfig/recipes/ppStack.config

@@ -34 +34 @@
 TARGET.MINMAG   F32     3.0             # Minimum magnitude difference to tolerate in stamp
 NMINPIX         S32     1               # Minimum input per pixel 
+BSCALEOFFSET    BOOL    FALSE           # HACK for removing offset from input warps because of compression (0.5*BSCALE)
 
 BACKGROUND.MODEL    BOOL   FALSE        # Construct a stacked version of the warp stage background

trunk/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
     );
 

trunk/ppStack/src/ppStackCombineFinal.c

@@ -44 +44 @@
 
 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");
@@ -107 +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);

trunk/ppStack/src/ppStackLoop.c

@@ -120 +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);
@@ -128 +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);
@@ -207 +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);

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

trunk/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;

trunk/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

trunk/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

trunk/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;
@@ -201 +203 @@
                          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);
@@ -271 +273 @@
         stack->data[i] = data;
 
+        //MEH -- apply bscale offset before norm   
+        if (bscaleApplyOffset) {
+            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

trunk/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);