Index: branches/eam_branches/ipp-ops-20130712/ppStack/src/ppStack.h =================================================================== --- branches/eam_branches/ipp-ops-20130712/ppStack/src/ppStack.h (revision 37064) +++ branches/eam_branches/ipp-ops-20130712/ppStack/src/ppStack.h (revision 37067) @@ -108,5 +108,6 @@ 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 ); Index: branches/eam_branches/ipp-ops-20130712/ppStack/src/ppStackCombineFinal.c =================================================================== --- branches/eam_branches/ipp-ops-20130712/ppStack/src/ppStackCombineFinal.c (revision 37064) +++ branches/eam_branches/ipp-ops-20130712/ppStack/src/ppStackCombineFinal.c (revision 37067) @@ -4,6 +4,45 @@ // #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,4 +59,7 @@ 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,4 +107,5 @@ 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,4 +152,10 @@ #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) Index: branches/eam_branches/ipp-ops-20130712/ppStack/src/ppStackLoop.c =================================================================== --- branches/eam_branches/ipp-ops-20130712/ppStack/src/ppStackLoop.c (revision 37064) +++ branches/eam_branches/ipp-ops-20130712/ppStack/src/ppStackLoop.c (revision 37067) @@ -87,5 +87,6 @@ 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,5 +120,5 @@ 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,5 +128,6 @@ 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,5 +183,6 @@ 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,5 +208,8 @@ // 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); Index: branches/eam_branches/ipp-ops-20130712/ppStack/src/ppStackLoop.h =================================================================== --- branches/eam_branches/ipp-ops-20130712/ppStack/src/ppStackLoop.h (revision 37064) +++ branches/eam_branches/ipp-ops-20130712/ppStack/src/ppStackLoop.h (revision 37067) @@ -65,5 +65,6 @@ bool safe, // Allow safe combination? bool norm, // Normalise images? - bool grow // Grow rejection masks? + bool grow, // Grow rejection masks? + bool bscaleoffset // Apply bscale offset? ); Index: branches/eam_branches/ipp-ops-20130712/ppStack/src/ppStackMatch.c =================================================================== --- branches/eam_branches/ipp-ops-20130712/ppStack/src/ppStackMatch.c (revision 37064) +++ branches/eam_branches/ipp-ops-20130712/ppStack/src/ppStackMatch.c (revision 37067) @@ -6,7 +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,4 +155,9 @@ 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. Index: branches/eam_branches/ipp-ops-20130712/ppStack/src/ppStackOptions.c =================================================================== --- branches/eam_branches/ipp-ops-20130712/ppStack/src/ppStackOptions.c (revision 37064) +++ branches/eam_branches/ipp-ops-20130712/ppStack/src/ppStackOptions.c (revision 37067) @@ -19,4 +19,5 @@ psFree(options->inputMask); psFree(options->sourceLists); + psFree(options->bscaleOffset); psFree(options->norm); psFree(options->sources); @@ -65,4 +66,5 @@ options->inputMask = NULL; options->sourceLists = NULL; + options->bscaleOffset = NULL; options->norm = NULL; options->sources = NULL; Index: branches/eam_branches/ipp-ops-20130712/ppStack/src/ppStackOptions.h =================================================================== --- branches/eam_branches/ipp-ops-20130712/ppStack/src/ppStackOptions.h (revision 37064) +++ branches/eam_branches/ipp-ops-20130712/ppStack/src/ppStackOptions.h (revision 37067) @@ -35,4 +35,5 @@ 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 Index: branches/eam_branches/ipp-ops-20130712/ppStack/src/ppStackPrepare.c =================================================================== --- branches/eam_branches/ipp-ops-20130712/ppStack/src/ppStackPrepare.c (revision 37064) +++ branches/eam_branches/ipp-ops-20130712/ppStack/src/ppStackPrepare.c (revision 37067) @@ -148,4 +148,9 @@ } + //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,4 +167,13 @@ } 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 Index: branches/eam_branches/ipp-ops-20130712/ppStack/src/ppStackReadout.c =================================================================== --- branches/eam_branches/ipp-ops-20130712/ppStack/src/ppStackReadout.c (revision 37064) +++ branches/eam_branches/ipp-ops-20130712/ppStack/src/ppStackReadout.c (revision 37067) @@ -34,4 +34,5 @@ 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,7 +41,8 @@ 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,4 +118,6 @@ 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,5 +165,5 @@ 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,4 +188,7 @@ psFree(stack); + //MEH change to trace + //psLogMsg("ppStack", PS_LOG_INFO, "initial stack image sectionNum %d", sectionNum); + sectionNum++; @@ -197,5 +204,5 @@ 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,4 +223,6 @@ 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,4 +274,11 @@ 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,5 +289,5 @@ 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,4 +309,7 @@ psFree(stack); + //MEH change to trace + //psLogMsg("ppStack", PS_LOG_INFO, "final stack image sectionNum %d", sectionNum); + sectionNum++; Index: branches/eam_branches/ipp-ops-20130712/ppStack/src/ppStackThread.c =================================================================== --- branches/eam_branches/ipp-ops-20130712/ppStack/src/ppStackThread.c (revision 37064) +++ branches/eam_branches/ipp-ops-20130712/ppStack/src/ppStackThread.c (revision 37067) @@ -277,5 +277,5 @@ { - psThreadTask *task = psThreadTaskAlloc("PPSTACK_FINAL_COMBINE", 6); + psThreadTask *task = psThreadTaskAlloc("PPSTACK_FINAL_COMBINE", 7); task->function = &ppStackReadoutFinalThread; psThreadTaskAdd(task);