Index: trunk/magic/remove/src/streaksremove.c =================================================================== --- trunk/magic/remove/src/streaksremove.c (revision 26408) +++ trunk/magic/remove/src/streaksremove.c (revision 26438) @@ -12,12 +12,16 @@ static pmConfig *parseArguments(int argc, char **argv); static bool readAndCopyToOutput(streakFiles *sf, bool exciseAll); -static void exciseNonWarpedPixels(streakFiles *sfiles, psImageMaskType newMaskValue); -static bool warpedPixel(streakFiles *sfiles, int x, int y); +static void exciseNonDiffedPixels(streakFiles *sfiles, psImageMaskType newMaskValue); +static bool diffedPixel(streakFiles *sfiles, int x, int y); static void excisePixel(streakFiles *sfiles, unsigned int x, unsigned int y, bool streak, psImageMaskType newMaskValue); static void writeImages(streakFiles *sf, bool exciseImageCube); static void updateAstrometry(streakFiles *sfiles); static void censorSources(streakFiles *sfiles, psImageMaskType maskStreak); -static long censorPixels(streakFiles *sfiles, psImage * pixels, bool checkNonWarpedPixels, psU16 maskStreak); - +static long censorPixels(streakFiles *sfiles, psImage * pixels, bool checkNonDiffedPixels, psU16 maskStreak); + +// Note: For clarity the flow of this program is in main(). +// There is not a lot of error checking is done in main. +// Until the end, where we might be doing Nebulous operations, called functions exit when an error +// is encountered. int main(int argc, char *argv[]) @@ -26,5 +30,5 @@ psLibInit(NULL); - psTimerStart("STREAKSREMOVE"); + psTimerStart("TOTAL_TIME"); pmConfig *config = parseArguments(argc, argv); @@ -46,4 +50,5 @@ streakFiles *sfiles = openFiles(config, true, argv[0]); setupAstrometry(sfiles); + sfiles->stats = psMetadataAlloc(); // Optionally we can set pixels that are masked to NAN since they couldn't have been @@ -57,12 +62,12 @@ bool exciseAll = false; - // --keepnonwarped is a test and debug mode - bool keepNonWarpedPixels = psMetadataLookupBool(&status, config->arguments, "KEEP_NON_WARPED"); - - // we need to check for non warped pixels unless we've been asked not to or the stage is diff - // (By definition pixels in diff images were included in a diff) - bool checkNonWarpedPixels = ! (keepNonWarpedPixels || (sfiles->stage == IPP_STAGE_DIFF) ); - - if (checkNonWarpedPixels ) { + // --keepnondiffed is a test and debug mode + bool keepNonDiffedPixels = psMetadataLookupBool(&status, config->arguments, "KEEP_NON_DIFFED"); + + // we need to check for non diffed pixels unless we've been asked not to or the stage is diff + // (By definition pixels in diff images were included in the difference images) + bool checkNonDiffedPixels = ! (keepNonDiffedPixels || (sfiles->stage == IPP_STAGE_DIFF) ); + + if (checkNonDiffedPixels ) { // From magic ICD: // In the raw and detrended images, the pixels which were not @@ -72,11 +77,12 @@ // if no skycells are provided sfiles->exciseAll is set to true - psTimerStart("COMPUTE_WARPED_PIXELS"); - if (! computeWarpedPixels(sfiles) ) { + psTimerStart("COMPUTE_DIFFED_PIXELS"); + if (! computeDiffedPixels(sfiles) ) { // we have no choice to excise all pixels exciseAll = true; } - psF64 cwp_t = psTimerClear("COMPUTE_WARPED_PIXELS"); - psLogMsg("streaksremove", PS_LOG_INFO, "time to compute warped pixels: %f\n", cwp_t); + psF64 cdp_t = psTimerClear("COMPUTE_DIFFED_PIXELS"); + psMetadataAddF64(sfiles->stats, PS_LIST_TAIL, "COMPUTE_UNDIFFED_PIXELS", PS_META_REPLACE, "time to compute non-diffedpixels", cdp_t); + psLogMsg("streaksremove", PS_LOG_INFO, "time to compute diffed pixels: %f\n", cdp_t); } @@ -97,4 +103,10 @@ long totalStreakPixels = 0; + // accumulators for the various timers + psF64 gsp_t = 0; + psF64 enw_t = 0; + psF64 rms_t = 0; + psF64 cs_t = 0; + psF64 wi_t = 0; // Iterate through each component of the input (except for raw images there is only one) do { @@ -126,4 +138,6 @@ StreakPixels *pixels = streak_on_component(streaks, sfiles->astrom, sfiles->inImage->numCols, sfiles->inImage->numRows); + gsp_t += psTimerClear("GET_STREAK_PIXELS"); + psMetadataAddF64(sfiles->stats, PS_LIST_TAIL, "GET_STREAK_PIXELS", PS_META_REPLACE, "", gsp_t); psLogMsg("streaksremove", PS_LOG_INFO, "time to get streak pixels: %f\n", psTimerClear("GET_STREAK_PIXELS")); @@ -131,19 +145,23 @@ // otherwise it contained an image cube (video cell) which is handled in the if block if (sfiles->inImage->image) { - if (checkNonWarpedPixels) { - psTimerStart("EXCISE_NON_WARPED"); - - // set non-warped pixels and variance to NAN, mask to maskStreak (since the pixel + if (checkNonDiffedPixels) { + psTimerStart("EXCISE_NON_DIFFED"); + + // set non-diffed pixels and variance to NAN, mask to maskStreak (since the pixel // is excised as part of the destreaking process) - exciseNonWarpedPixels(sfiles, sfiles->maskStreak); - - psLogMsg("streaksremove", PS_LOG_INFO, "time to excise non warped pixels: %f\n", psTimerClear("EXCISE_NON_WARPED")); + exciseNonDiffedPixels(sfiles, sfiles->maskStreak); + + enw_t += psTimerClear("EXCISE_NON_DIFFED"); + psMetadataAddF64(sfiles->stats, PS_LIST_TAIL, "EXCISE_NON_DIFFED", PS_META_REPLACE, "", enw_t); + psLogMsg("streaksremove", PS_LOG_INFO, "time to excise non diffed pixels: %f\n", enw_t); } psTimerStart("REMOVE_STREAKS"); - totalStreakPixels += censorPixels(sfiles, pixels, checkNonWarpedPixels, sfiles->maskStreak); - - psLogMsg("streaksremove", PS_LOG_INFO, "time to remove streak pixels: %f\n", psTimerClear("REMOVE_STREAKS")); + totalStreakPixels += censorPixels(sfiles, pixels, checkNonDiffedPixels, sfiles->maskStreak); + + rms_t += psTimerClear("REMOVE_STREAKS"); + psMetadataAddF64(sfiles->stats, PS_LIST_TAIL, "REMOVE_STREAKS", PS_META_REPLACE, "", enw_t); + psLogMsg("streaksremove", PS_LOG_INFO, "time to remove streak pixels: %f\n", rms_t); if (nanForRelease) { @@ -167,9 +185,15 @@ } + psTimerStart("CENSOR_SOURCES"); censorSources(sfiles, sfiles->maskStreak); + cs_t += psTimerClear("CENSOR_SOURCES"); + psMetadataAddF64(sfiles->stats, PS_LIST_TAIL, "CENSOR_SOURCES", PS_META_REPLACE, "", cs_t); + // write the destreaked "temporary" images and the recovery images + psTimerStart("WRITE_IMAGES"); writeImages(sfiles, exciseImageCube); - + wi_t += psTimerClear("WRITE_IMAGES"); + psMetadataAddF64(sfiles->stats, PS_LIST_TAIL, "WRITE_IMAGES", PS_META_REPLACE, "", wi_t); psLogMsg("streaksremove", PS_LOG_INFO, "time to process component %d: %f\n", sfiles->extnum, psTimerClear("PROCESS_COMPONENT")); @@ -179,5 +203,11 @@ psFree(streaks); - psLogMsg("streaksremove", PS_LOG_INFO, "pixels: %ld streak pixels: %ld %4.2f%%\n", totalPixels, totalStreakPixels, 100. * totalStreakPixels / totalPixels); + if (exciseAll) { + totalStreakPixels = totalPixels; + } + + psF64 maskedFraction = totalStreakPixels / totalPixels; + psLogMsg("streaksremove", PS_LOG_INFO, "pixels: %ld streak pixels: %ld %4.2f%%\n", totalPixels, totalStreakPixels, maskedFraction * 100); + psMetadataAddF64(sfiles->stats, PS_LIST_TAIL, "MASKED_FRACTION", PS_META_REPLACE, "", maskedFraction); // check the weight and mask files for extra extensions that might be in files @@ -188,10 +218,11 @@ psTimerStart("CLOSE_IMAGES"); - closeImages(sfiles); - - psLogMsg("streaksremove", PS_LOG_INFO, "time to close images: %f\n", psTimerClear("CLOSE_IMAGES")); - - + psF64 ci_t = psTimerClear("CLOSE_IMAGES"); + psMetadataAddF64(sfiles->stats, PS_LIST_TAIL, "CLOSE_IMAGES", PS_META_REPLACE, "", ci_t); + + psLogMsg("streaksremove", PS_LOG_INFO, "time to close images: %f\n", ci_t); + + psTimerStart("REPLICATE_OUTPUTS"); if (!replicateOutputs(sfiles)) { psErrorStackPrint(stderr, "failed to replicate output files"); @@ -199,4 +230,6 @@ exit(PS_EXIT_UNKNOWN_ERROR); } + psF64 ro_t = psTimerClear("REPLICATE_OUTPUTS"); + psMetadataAddF64(sfiles->stats, PS_LIST_TAIL, "REPLICATE_OUTPUTS", PS_META_REPLACE, "", ro_t); // NOTE: from here on we can't just quit if something goes wrong. @@ -208,9 +241,9 @@ // swap the instances for the input and output // Note this is a nebulous database operation. No file I/O is performed + psTimerStart("SWAP_INSTANCES"); if (!swapOutputsToInputs(sfiles)) { - // XXX: Now what? // It is up to the program that reverts failed destreak runs to insure that - // any input files that have been swapped are restored and that the de-streaked - // versions are deleted + // any original non-destreaked input files that have been swapped are restored and that the de-streaked + // versions are deleted. psErrorStackPrint(stderr, "failed to swap files"); @@ -218,4 +251,24 @@ // XXX: pick a specific error code for this failure exit(PS_EXIT_UNKNOWN_ERROR); + } + psF64 si_t = psTimerClear("SWAP_INSTANCES"); + psMetadataAddF64(sfiles->stats, PS_LIST_TAIL, "SWAP_INSTANCES", PS_META_REPLACE, "", si_t); + } + + psF64 total_time = psTimerClear("TOTAL_TIME"); + psMetadataAddF64(sfiles->stats, PS_LIST_TAIL, "TOTAL_TIME", PS_META_REPLACE, "", total_time); + psLogMsg("streaksremove", PS_LOG_INFO, "time to run streaksremove: %f\n", total_time); + + if (sfiles->statsFile) { + const char *statsMDC = psMetadataConfigFormat(sfiles->stats); + if (!statsMDC || strlen(statsMDC) == 0) { + psError(PS_ERR_IO, false, "Unable to get statistics MDC file.\n"); + } else { + fprintf(sfiles->statsFile, "%s", statsMDC); + psFree((void *)statsMDC); + fclose(sfiles->statsFile); + sfiles->statsFile = NULL; + psFree(sfiles->stats); + sfiles->stats = NULL; } } @@ -228,5 +281,4 @@ streaksNebulousCleanup(); pmConfigDone(); - psLogMsg("streaksremove", PS_LOG_INFO, "time to run streaksremove: %f\n", psTimerClear("STREAKSREMOVE")); psLibFinalize(); @@ -237,5 +289,5 @@ static long -censorPixels(streakFiles *sfiles, psImage *pixels, bool checkNonWarpedPixels, psU16 maskStreak) +censorPixels(streakFiles *sfiles, psImage *pixels, bool checkNonDiffedPixels, psU16 maskStreak) { long streakPixels = 0; @@ -245,5 +297,5 @@ if (psImageGet(pixels, x, y)) { ++streakPixels; - if (!checkNonWarpedPixels || warpedPixel(sfiles, x, y)) { + if (!checkNonDiffedPixels || diffedPixel(sfiles, x, y)) { excisePixel(sfiles, x, y, true, maskStreak); @@ -251,5 +303,5 @@ } else { // This pixel was not included in any warp and has thus already excised - // by exciseNonWarpedPixels + // by exciseNonDiffedPixels } } @@ -268,5 +320,5 @@ fprintf(stderr, "\t-weight WEIGHT.fits: weight file to de-streak\n"); fprintf(stderr, "\t-replace: replace the input images with the output\n"); - fprintf(stderr, "\t-keepnonwarped: do not exise pixels that were not part of difference processing\n"); + fprintf(stderr, "\t-keepnondiffed: do not exise pixels that were not part of difference processing\n"); fprintf(stderr, "\t-transparent val: instead of setting excicsed pixel to NAN add val\n"); fprintf(stderr, "\t-chip_mask MASK.fits: name of mask for chip stage (camera stage mask is passed tih -mask)\n"); @@ -359,8 +411,8 @@ } - if ((argnum = psArgumentGet(argc, argv, "-keepnonwarped"))) { - psArgumentRemove(argnum, &argc, argv); - psMetadataAddBool(config->arguments, PS_LIST_TAIL, "KEEP_NON_WARPED", 0, - "skip excising of non warped pixels", true); + if ((argnum = psArgumentGet(argc, argv, "-keepnondiffed"))) { + psArgumentRemove(argnum, &argc, argv); + psMetadataAddBool(config->arguments, PS_LIST_TAIL, "KEEP_NON_DIFFED", 0, + "skip excising of non diffed pixels", true); } @@ -377,5 +429,5 @@ } - // if skycells are not provided then we have to execise all pixels unless -keepnonwarped + // if skycells are not provided then we have to execise all pixels unless -keepnondiffed pmConfigFileSetsMD(config->arguments, &argc, argv, "SKYCELLS", "-skycell", "-skycelllist"); @@ -440,4 +492,11 @@ psMetadataAddStr(config->arguments, PS_LIST_TAIL, "SOURCES", 0, "name of input sources file", argv[argnum]); + psArgumentRemove(argnum, &argc, argv); + } + + if ((argnum = psArgumentGet(argc, argv, "-stats"))) { + psArgumentRemove(argnum, &argc, argv); + psMetadataAddStr(config->arguments, PS_LIST_TAIL, "STATS", 0, + "name of input stats file", argv[argnum]); psArgumentRemove(argnum, &argc, argv); } @@ -589,4 +648,5 @@ strkGetMaskValues(sf); + // add the STREAK bit to the mask image pixels setStreakBits(sf->inMask->image, sf->maskStreak); } @@ -721,11 +781,10 @@ sfiles->inMask->image->data.PS_TYPE_IMAGE_MASK_DATA[y][x]; } - sfiles->outMask->image->data.PS_TYPE_IMAGE_MASK_DATA[y][x] = - sfiles->inMask->image->data.PS_TYPE_IMAGE_MASK_DATA[y][x] | newMaskValue; + sfiles->outMask->image->data.PS_TYPE_IMAGE_MASK_DATA[y][x] |= newMaskValue; } } static void -exciseNonWarpedPixels(streakFiles *sfiles, psImageMaskType newMaskValue) +exciseNonDiffedPixels(streakFiles *sfiles, psImageMaskType newMaskValue) { int cell_x0 = sfiles->astrom->cell_x0; @@ -749,5 +808,5 @@ } - psU8 *pixels = sfiles->warpedPixels->data.U8[yChip]; + psU8 *pixels = sfiles->diffedPixels->data.U8[yChip]; if (xParity == 1) { @@ -772,10 +831,10 @@ static bool -warpedPixel(streakFiles *sfiles, int x, int y) +diffedPixel(streakFiles *sfiles, int x, int y) { PixelPos chipCoord; if (!CHIP_LEVEL_INPUT(sfiles->stage)) { - // if we're here on a skycell image by definition this pixel was warped + // if we're here on a skycell image by definition this pixel was diffed return true; } @@ -790,12 +849,12 @@ cellToChipInt(&chipCoord.x, &chipCoord.y, sfiles->astrom, x, y); - if (chipCoord.x < 0 || chipCoord.x >= sfiles->warpedPixels->numCols) { + if (chipCoord.x < 0 || chipCoord.x >= sfiles->diffedPixels->numCols) { return false; } - if (chipCoord.y < 0 || chipCoord.y >= sfiles->warpedPixels->numRows) { + if (chipCoord.y < 0 || chipCoord.y >= sfiles->diffedPixels->numRows) { return false; } - return psImageGet(sfiles->warpedPixels, chipCoord.x, chipCoord.y) ? true : false; + return psImageGet(sfiles->diffedPixels, chipCoord.x, chipCoord.y) ? true : false; }