Index: /branches/eam_branches/20091201/magic/Makefile =================================================================== --- /branches/eam_branches/20091201/magic/Makefile (revision 26869) +++ /branches/eam_branches/20091201/magic/Makefile (revision 26870) @@ -2,12 +2,11 @@ default: all -all: ssa-core-cpp magic remove +all: ssa-core-cpp magic remove verify -install: magic.install remove.install +install: magic.install remove.install verify.install clean: if [ -d magic ]; then (cd magic && make -f Makefile.magic clean); fi if [ -d ssa-core-cpp ]; then (cd ssa-core-cpp && make clean); fi - cd $(REMOVE) && make -f Makefile.simple clean update: Index: /branches/eam_branches/20091201/magic/censorObjects/src/censor.c =================================================================== --- /branches/eam_branches/20091201/magic/censorObjects/src/censor.c (revision 26869) +++ /branches/eam_branches/20091201/magic/censorObjects/src/censor.c (revision 26870) @@ -14,5 +14,8 @@ static void usage(void) { - fprintf(stderr, "USAGE: censor [-file image(s)] [-list imagelist] (output)\n"); + fprintf(stderr, "USAGE: censorObjects -file (sources) -masklist (masklist) (output)\n"); + fprintf(stderr, " censorObjects -file (sources) -mask (mask.fits) (output)\n"); + fprintf(stderr, " censorObjects -list (sourceslist) -masklist (masklist) (output)\n"); + fprintf(stderr, " censorObjects -list (sourceslist) -mask (mask.fits) (output)\n"); exit(PS_EXIT_CONFIG_ERROR); } Index: /branches/eam_branches/20091201/magic/remove/src/Makefile.simple =================================================================== --- /branches/eam_branches/20091201/magic/remove/src/Makefile.simple (revision 26869) +++ /branches/eam_branches/20091201/magic/remove/src/Makefile.simple (revision 26870) @@ -12,5 +12,5 @@ streaksremove.o \ streaksextern.o \ - warpedpixels.o \ + diffedpixels.o \ Line.o Index: /branches/eam_branches/20091201/magic/remove/src/isdestreaked.c =================================================================== --- /branches/eam_branches/20091201/magic/remove/src/isdestreaked.c (revision 26869) +++ /branches/eam_branches/20091201/magic/remove/src/isdestreaked.c (revision 26870) @@ -1,3 +1,5 @@ #include "streaksremove.h" + +char *streaksProgram = "isdestreaked"; int Index: /branches/eam_branches/20091201/magic/remove/src/streaksastrom.c =================================================================== --- /branches/eam_branches/20091201/magic/remove/src/streaksastrom.c (revision 26869) +++ /branches/eam_branches/20091201/magic/remove/src/streaksastrom.c (revision 26870) @@ -379,6 +379,5 @@ { if (!pmAstromLinearizeTransforms((pmFPA *) astrom->fpa, (pmChip *) astrom->chip)) { - // streaksExit("linear fit to astrometry failed\n", PS_EXIT_UNKNOWN_ERROR); - psError(PS_ERR_UNKNOWN, false, "linear fit to astrometry failed ignoring"); + psErrorStackPrint(stderr, "linear fit to astrometry failed. ignoring\n"); return false; } Index: /branches/eam_branches/20091201/magic/remove/src/streakscompare.c =================================================================== --- /branches/eam_branches/20091201/magic/remove/src/streakscompare.c (revision 26869) +++ /branches/eam_branches/20091201/magic/remove/src/streakscompare.c (revision 26870) @@ -1,3 +1,5 @@ #include "streaksremove.h" + +char *streaksProgram = "streakscompare"; static pmConfig * parseArguments(int, char **); Index: /branches/eam_branches/20091201/magic/remove/src/streaksio.c =================================================================== --- /branches/eam_branches/20091201/magic/remove/src/streaksio.c (revision 26869) +++ /branches/eam_branches/20091201/magic/remove/src/streaksio.c (revision 26870) @@ -121,4 +121,14 @@ sf->transparentStreaks = psMetadataLookupF64(&status, config->arguments, "TRANSPARENT_STREAKS"); + sf->stats = psMetadataAlloc(); + psString statsFileName= psMetadataLookupStr(&status, config->arguments, "STATS"); + if (statsFileName) { + sf->statsFile = fopen(statsFileName, "w"); + if (!sf->statsFile) { + psError(PS_ERR_IO, true, "failed to open stats file %s", statsFileName); + streaksExit("", PS_EXIT_CONFIG_ERROR); + } + } + return sf; } @@ -128,5 +138,5 @@ { freeImages(sf); - psFree(sf->warpedPixels); + psFree(sf->diffedPixels); psFree(sf->tiles); psFree(sf->view); @@ -800,4 +810,5 @@ return; } + streaksVersionHeaderFull(sfile->header); if (!psFitsWriteImage(sfile->fits, sfile->header, sfile->image, 0, extname)) { psError(PS_ERR_IO, false, "failed to write image to %s extnum: %d", @@ -817,4 +828,5 @@ return; } + streaksVersionHeaderFull(sfile->header); if (!psFitsWriteImageCube(sfile->fits, sfile->header, imagecube, extname)) { psError(PS_ERR_IO, false, "failed to write image to %s extnum: %d", @@ -851,7 +863,7 @@ sfile->header = NULL; psFree(sfile->image); - sfile->header = NULL; + sfile->image = NULL; psFree(sfile->imagecube); - sfile->header = NULL; + sfile->imagecube = NULL; } @@ -1181,7 +1193,14 @@ } -void -strkGetMaskValues(streakFiles *sfiles, psU32 *maskStreak, psU32 *maskMask) -{ +// Get the mask values that we need. +// If an input mask file is provided get them from there. +// Otherwise get them from the recipe MASKS +void +strkGetMaskValues(streakFiles *sfiles) +{ + if (sfiles->maskStreak) { + // already done + return; + } if (sfiles->inMask && sfiles->inMask->header) { if (!pmConfigMaskReadHeader(sfiles->config, sfiles->inMask->header)) { @@ -1196,5 +1215,5 @@ streaksExit("", PS_EXIT_CONFIG_ERROR); } - *maskStreak = psMetadataLookupU32(&status, masks, "STREAK"); + sfiles->maskStreak = psMetadataLookupU32(&status, masks, "STREAK"); if (!status) { psError(PM_ERR_CONFIG, false, "failed to lookup mask value for STREAK in recipes\n"); @@ -1212,5 +1231,5 @@ } } - *maskMask = ~convPoor; + sfiles->maskMask = ~convPoor; } Index: /branches/eam_branches/20091201/magic/remove/src/streaksio.h =================================================================== --- /branches/eam_branches/20091201/magic/remove/src/streaksio.h (revision 26869) +++ /branches/eam_branches/20091201/magic/remove/src/streaksio.h (revision 26870) @@ -16,5 +16,5 @@ void copyFitsOptions(sFile *out, sFile *rec, sFile *in, psVector *tiles); void setupImageRefs(sFile *out, sFile *recoveryOut, sFile *in, int extnum, bool exciseAll); -void strkGetMaskValues(streakFiles *sfiles, psU32 *maskStreak, psU32 *maskMask); +void strkGetMaskValues(streakFiles *sfiles); void copyExtraExtensions(streakFiles *sfiles); Index: /branches/eam_branches/20091201/magic/remove/src/streaksrelease.c =================================================================== --- /branches/eam_branches/20091201/magic/remove/src/streaksrelease.c (revision 26869) +++ /branches/eam_branches/20091201/magic/remove/src/streaksrelease.c (revision 26870) @@ -4,4 +4,6 @@ static bool readAndCopyToOutput(streakFiles *sf, bool exciseAll); static void writeImages(streakFiles *sf, bool exciseImageCube); + +char *streaksProgram = "streaksrelease"; int @@ -16,8 +18,4 @@ return PS_EXIT_CONFIG_ERROR; } - - // Values to set for masked pixels - psU32 maskStreak = 0; // for the image and weight (usually NAN, MAXINT for integer images) - psU32 maskMask = 0; // value looked up for MASK.STREAK // Does true work here? @@ -50,9 +48,7 @@ // now that we've read the input files, lookup the mask values that we read - if (maskStreak == 0) { - strkGetMaskValues(sfiles, &maskStreak, &maskMask); - } - - setMaskedToNAN(sfiles, maskMask, true); + strkGetMaskValues(sfiles); + + setMaskedToNAN(sfiles, sfiles->maskMask, true); // write out the destreaked temporary images and the recovery images Index: /branches/eam_branches/20091201/magic/remove/src/streaksremove.c =================================================================== --- /branches/eam_branches/20091201/magic/remove/src/streaksremove.c (revision 26869) +++ /branches/eam_branches/20091201/magic/remove/src/streaksremove.c (revision 26870) @@ -12,12 +12,18 @@ 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); + +char *streaksProgram = "streaksremove"; + +// 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 +32,5 @@ psLibInit(NULL); - psTimerStart("STREAKSREMOVE"); + psTimerStart("TOTAL_TIME"); pmConfig *config = parseArguments(argc, argv); @@ -33,8 +39,4 @@ streaksExit("", PS_EXIT_CONFIG_ERROR); } - - // Values to set for masked pixels - psU32 maskStreak = 0; // for the image and weight (usually NAN, MAXINT for integer images) - psU32 maskMask = 0; // value looked up for MASK.STREAK psString streaksFileName = psMetadataLookupStr(NULL, config->arguments, "STREAKS"); @@ -50,4 +52,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 @@ -61,12 +64,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 @@ -76,11 +79,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"); + psMetadataAddF32(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); } @@ -101,4 +105,11 @@ 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; + psF64 ua_t = 0; // Iterate through each component of the input (except for raw images there is only one) do { @@ -118,7 +129,5 @@ // now that we've read the input files, lookup the mask values - if (maskStreak == 0) { - strkGetMaskValues(sfiles, &maskStreak, &maskMask); - } + strkGetMaskValues(sfiles); totalPixels += sfiles->inImage->numRows * sfiles->inImage->numCols; @@ -132,4 +141,6 @@ StreakPixels *pixels = streak_on_component(streaks, sfiles->astrom, sfiles->inImage->numCols, sfiles->inImage->numRows); + gsp_t += psTimerClear("GET_STREAK_PIXELS"); + psMetadataAddF32(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")); @@ -137,23 +148,27 @@ // 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, 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"); + psMetadataAddF32(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, 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"); + psMetadataAddF32(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) { // set any pixels that were masked, to NAN (unless they are already NAN) - setMaskedToNAN(sfiles, maskMask, true); + setMaskedToNAN(sfiles, sfiles->maskMask, true); } @@ -170,14 +185,24 @@ // chip processed files with the data calcuated by psastro at the camera stage // (actually we use a linear approximation) + psTimerStart("UPDATE_ASTROMETRY"); updateAstrometry(sfiles); - } - - censorSources(sfiles, maskStreak); + ua_t += psTimerClear("UPDATE_ASTROMETRY"); + psMetadataAddF32(sfiles->stats, PS_LIST_TAIL, "UPDATE_ASTROMETRY", PS_META_REPLACE, "", ua_t); + psLogMsg("streaksremove", PS_LOG_INFO, "time to update astrometry: %f\n", ua_t); + } + + psTimerStart("CENSOR_SOURCES"); + censorSources(sfiles, sfiles->maskStreak); + cs_t += psTimerClear("CENSOR_SOURCES"); + psMetadataAddF32(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"); + psMetadataAddF32(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")); + } while (streakFilesNextExtension(sfiles)); @@ -185,5 +210,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 streakFraction = (double) totalStreakPixels / totalPixels; + psLogMsg("streaksremove", PS_LOG_INFO, "pixels: %ld streak pixels: %ld %4.2f%%\n", totalPixels, totalStreakPixels, streakFraction * 100); + psMetadataAddF32(sfiles->stats, PS_LIST_TAIL, "STREAK_FRACTION", PS_META_REPLACE, "", streakFraction); // check the weight and mask files for extra extensions that might be in files @@ -194,10 +225,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"); + psMetadataAddF32(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"); @@ -205,4 +237,6 @@ exit(PS_EXIT_UNKNOWN_ERROR); } + psF64 ro_t = psTimerClear("REPLICATE_OUTPUTS"); + psMetadataAddF32(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. @@ -214,9 +248,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"); @@ -224,4 +258,24 @@ // XXX: pick a specific error code for this failure exit(PS_EXIT_UNKNOWN_ERROR); + } + psF64 si_t = psTimerClear("SWAP_INSTANCES"); + psMetadataAddF32(sfiles->stats, PS_LIST_TAIL, "SWAP_INSTANCES", PS_META_REPLACE, "", si_t); + } + + psF64 total_time = psTimerClear("TOTAL_TIME"); + psMetadataAddF32(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(statsMDC); + fclose(sfiles->statsFile); + sfiles->statsFile = NULL; + psFree(sfiles->stats); + sfiles->stats = NULL; } } @@ -234,5 +288,4 @@ streaksNebulousCleanup(); pmConfigDone(); - psLogMsg("streaksremove", PS_LOG_INFO, "time to run streaksremove: %f\n", psTimerClear("STREAKSREMOVE")); psLibFinalize(); @@ -243,5 +296,5 @@ static long -censorPixels(streakFiles *sfiles, psImage *pixels, bool checkNonWarpedPixels, psU16 maskStreak) +censorPixels(streakFiles *sfiles, psImage *pixels, bool checkNonDiffedPixels, psU16 maskStreak) { long streakPixels = 0; @@ -251,5 +304,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); @@ -257,5 +310,5 @@ } else { // This pixel was not included in any warp and has thus already excised - // by exciseNonWarpedPixels + // by exciseNonDiffedPixels } } @@ -274,5 +327,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"); @@ -365,8 +418,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); } @@ -378,10 +431,10 @@ psArgumentRemove(argnum, &argc, argv); double transparentStreaks = atof(argv[argnum]); - psMetadataAddF64(config->arguments, PS_LIST_TAIL, "TRANSPARENT_STREAKS", 0, + psMetadataAddF32(config->arguments, PS_LIST_TAIL, "TRANSPARENT_STREAKS", 0, "value to adjust excised pixels", transparentStreaks); psArgumentRemove(argnum, &argc, argv); } - // 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"); @@ -449,4 +502,11 @@ } + 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); + } + if ((argnum = psArgumentGet(argc, argv, "-tmproot"))) { psArgumentRemove(argnum, &argc, argv); @@ -491,21 +551,28 @@ updateAstrometry(streakFiles *sf) { - // XXX: why do I check this here? Shouldn't it be just around the call to linearizeTransforms? if (sf->bilevelAstrometry) { - if (!linearizeTransforms(sf->astrom)) { - // fit failed, leave the astrometry unchanged + // fit failed, leave the transform in the file unchanged return; } - - if (!pmAstromWriteWCS(sf->outImage->header, sf->inAstrom->fpa, sf->chip, 0.001)) { - psError(PS_ERR_UNKNOWN, false, "failed to update astrometry for extension %d", sf->extnum); - streaksExit("", PS_EXIT_UNKNOWN_ERROR); - } - if (sf->outMask) { - pmAstromWriteWCS(sf->outMask->header, sf->inAstrom->fpa, sf->chip, 0.001); - } - if (sf->outWeight) { - pmAstromWriteWCS(sf->outWeight->header, sf->inAstrom->fpa, sf->chip, 0.001); + } + if (!pmAstromWriteWCS(sf->outImage->header, sf->inAstrom->fpa, sf->chip, 0.001)) { + psError(PS_ERR_UNKNOWN, false, "failed to update astrometry for extension %d", sf->extnum); + streaksExit("", PS_EXIT_UNKNOWN_ERROR); + } + if (sf->outMask) { + pmAstromWriteWCS(sf->outMask->header, sf->inAstrom->fpa, sf->chip, 0.001); + } + if (sf->outWeight) { + pmAstromWriteWCS(sf->outWeight->header, sf->inAstrom->fpa, sf->chip, 0.001); + } +} + +static void +setStreakBits(psImage *maskImage, psU32 maskStreak) +{ + for (int y=0 ; y < maskImage->numRows; y++) { + for (int x=0 ; x < maskImage->numCols; x++) { + maskImage->data.PS_TYPE_IMAGE_MASK_DATA[y][x] |= maskStreak; } } @@ -516,5 +583,4 @@ readAndCopyToOutput(streakFiles *sf, bool exciseAll) { - bool updateAstrometry = false; if (sf->inImage->pmfile) { // image data from pmFPAfile (diff or warp file) @@ -534,15 +600,4 @@ streaksExit("", PS_EXIT_UNKNOWN_ERROR); } - // For the chip level files, copy the WCS from the astrometry file to the header - // XXX: do we want to do this for raw images as well? - if (sf->stage == IPP_STAGE_CHIP) { - if (!sf->bilevelAstrometry) { - updateAstrometry = true; - if (!pmAstromWriteWCS(sf->inImage->header, sf->inAstrom->fpa, sf->chip, 0.001)) { - psError(PS_ERR_UNKNOWN, false, "failed to update astrometry for extension %d", sf->extnum); - streaksExit("", PS_EXIT_UNKNOWN_ERROR); - } - } - } sf->outImage->header = psMemIncrRefCounter(sf->inImage->header); if (sf->recImage) { @@ -577,8 +632,12 @@ } addDestreakKeyword(sf->outMask->header); - if (updateAstrometry) { - pmAstromWriteWCS(sf->outMask->header, sf->inAstrom->fpa, sf->chip, 0.001); - } - setupImageRefs(sf->outMask, sf->recMask, sf->inMask, sf->extnum, exciseAll); + // Note: we don't excise the mask pixels even if exciseAll is true. + setupImageRefs(sf->outMask, sf->recMask, sf->inMask, sf->extnum, false); + if (exciseAll) { + strkGetMaskValues(sf); + + // add the STREAK bit to the mask image pixels + setStreakBits(sf->inMask->image, sf->maskStreak); + } if (sf->outChMask) { sf->outChMask->header = (psMetadata *) psMemIncrRefCounter(sf->outMask->header); @@ -605,7 +664,4 @@ } addDestreakKeyword(sf->outWeight->header); - if (updateAstrometry) { - pmAstromWriteWCS(sf->inWeight->header, sf->inAstrom->fpa, sf->chip, 0.001); - } setupImageRefs(sf->outWeight, sf->recWeight, sf->inWeight, sf->extnum, exciseAll); @@ -630,27 +686,31 @@ } else { // we have an image cube - double initValue; if (exciseImageCube) { // copy the entire input image to the recovery image writeImageCube(sf->recImage, sf->inImage->imagecube, extname, sf->extnum); - initValue = NAN; } else { // otherwise write it to the output writeImageCube(sf->outImage, sf->inImage->imagecube, extname, sf->extnum); - initValue = 0; - } - - // borrow one of the images from the imagecube and set it to init value - psImage *image = psArrayGet (sf->inImage->imagecube, 0); - psMemIncrRefCounter(image); - psImageInit(image, initValue); + } + + // Now deal with the other output image if (exciseImageCube) { - sf->outImage->image = image; - writeImage(sf->outImage, extname, sf->extnum); + // Set the values in the imagecube images to NAN and write them to the output image + for (int i = 0; i < psArrayLength(sf->inImage->imagecube); i++) { + psImage *image = psArrayGet (sf->inImage->imagecube, i); + // XXX: NAN isn't right. It should be the integer equivalent. Should use exciseValue + // but it isn't set with this code path. Fix that. + psImageInit(image, 65535); + } + writeImageCube(sf->outImage, sf->inImage->imagecube, extname, sf->extnum); } else { - // write zero valued image to reccovery if (sf->recImage) { - sf->recImage->image = image; - writeImage(sf->recImage, extname, sf->extnum); + // Set the values in the imagecube images to zero + for (int i = 0; i < psArrayLength(sf->inImage->imagecube); i++) { + psImage *image = psArrayGet (sf->inImage->imagecube, i); + psImageInit(image, 0); + } + // copy the entire zeroed image to the recovery image + writeImageCube(sf->recImage, sf->inImage->imagecube, extname, sf->extnum); } } @@ -683,17 +743,35 @@ } - float imageValue = sfiles->inImage->image->data.F32[y][x]; - if (sfiles->recImage && !isExciseValue(imageValue, sfiles->inImage->exciseValue) ) { - sfiles->recImage->image->data.F32[y][x] = imageValue; - } - - if (sfiles->transparentStreaks == 0) { - sfiles->outImage->image->data.F32[y][x] = exciseValue; + if (sfiles->inImage->image->type.type == PS_TYPE_U16) { + psU16 imageValue = sfiles->inImage->image->data.U16[y][x]; + if (sfiles->recImage && !isExciseValue(imageValue, sfiles->inImage->exciseValue) ) { + sfiles->recImage->image->data.U16[y][x] = imageValue; + } + + if (sfiles->transparentStreaks == 0) { + sfiles->outImage->image->data.U16[y][x] = exciseValue; + } else { + if (streak) { + // as a visualization aid don't mask the pixel, just change the intensity + sfiles->outImage->image->data.U16[y][x] = imageValue + sfiles->transparentStreaks; + } else { + sfiles->outImage->image->data.U16[y][x] = exciseValue; + } + } } else { - if (streak) { - // as a visualization aid don't mask the pixel, just change the intensity - sfiles->outImage->image->data.F32[y][x] = imageValue + sfiles->transparentStreaks; + float imageValue = sfiles->inImage->image->data.F32[y][x]; + if (sfiles->recImage && !isExciseValue(imageValue, sfiles->inImage->exciseValue) ) { + sfiles->recImage->image->data.F32[y][x] = imageValue; + } + + if (sfiles->transparentStreaks == 0) { + sfiles->outImage->image->data.F32[y][x] = exciseValue; } else { - sfiles->outImage->image->data.F32[y][x] = exciseValue; + if (streak) { + // as a visualization aid don't mask the pixel, just change the intensity + sfiles->outImage->image->data.F32[y][x] = imageValue + sfiles->transparentStreaks; + } else { + sfiles->outImage->image->data.F32[y][x] = exciseValue; + } } } @@ -711,11 +789,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; @@ -739,5 +816,5 @@ } - psU8 *pixels = sfiles->warpedPixels->data.U8[yChip]; + psU8 *pixels = sfiles->diffedPixels->data.U8[yChip]; if (xParity == 1) { @@ -762,10 +839,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; } @@ -780,12 +857,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; } Index: /branches/eam_branches/20091201/magic/remove/src/streaksremove.h =================================================================== --- /branches/eam_branches/20091201/magic/remove/src/streaksremove.h (revision 26869) +++ /branches/eam_branches/20091201/magic/remove/src/streaksremove.h (revision 26870) @@ -33,6 +33,9 @@ int yParity; double exciseValue; + psString program; } sFile; +// used for error messages +extern char * streaksProgram; typedef enum { @@ -64,4 +67,6 @@ sFile *inSources; sFile *outSources; + psMetadata *stats; + FILE *statsFile; psString class_id; pmFPAfile *inAstrom; @@ -71,8 +76,10 @@ pmChip *chip; // current chip pmCell *cell; // current cell - psImage *warpedPixels; + psImage *diffedPixels; psVector *tiles; double transparentStreaks; psString program_name; + psU32 maskStreak; + psU32 maskMask; } streakFiles; @@ -86,9 +93,12 @@ extern void cellToChipInt(unsigned int *xChip, unsigned int *yChip, strkAstrom *astrom, int xCell, int yCell); -extern bool computeWarpedPixels(streakFiles *sf); +extern bool computeDiffedPixels(streakFiles *sf); extern void streaksExit(psString, int); extern ippStage parseStage(psString); extern psString pathToDirectory(char *path); extern void setStreakFiles( streakFiles *); + +extern bool streaksVersionHeaderFull(psMetadata *header); +extern psString streaksVersionLong(void); #define CHIP_LEVEL_INPUT(_stage) ((_stage == IPP_STAGE_RAW) || (_stage == IPP_STAGE_CHIP)) Index: /branches/eam_branches/20091201/magic/remove/src/streaksreplace.c =================================================================== --- /branches/eam_branches/20091201/magic/remove/src/streaksreplace.c (revision 26869) +++ /branches/eam_branches/20091201/magic/remove/src/streaksreplace.c (revision 26870) @@ -5,4 +5,6 @@ static void writeImages(streakFiles *sf, bool exciseImageCube); static bool readAndCopyToOutput(streakFiles *sf, bool restoreImageCube); + +char *streaksProgram = "streaksreplace"; int main(int argc, char *argv[]) Index: anches/eam_branches/20091201/magic/remove/src/warpedpixels.c =================================================================== --- /branches/eam_branches/20091201/magic/remove/src/warpedpixels.c (revision 26869) +++ (revision ) @@ -1,330 +1,0 @@ -#include "streaksremove.h" -#include "pmAstrometryWCS.h" - -// #define DEBUG_PRINT 1 - -static void addTouchedPixels(streakFiles *sf, psString fileName); -static void nameCorners(psPlane pt[4]); -static int xLeft(psPlane *pt, int y); -static int xRight(psPlane *pt, int y); - -bool -computeWarpedPixels(streakFiles *sf) -{ - bool status; - pmConfig *config = sf->config; - - psArray *skycells = psMetadataLookupPtr(&status, config->arguments, "SKYCELLS"); - if (skycells == NULL) { - return false; - } - - int n = psArrayLength(skycells); - - psRegion *bounds = pmChipPixels(sf->chip); - - int width = bounds->x1 - bounds->x0; - int height = bounds->y1 - bounds->y0; - - psFree(bounds); - - sf->warpedPixels = psImageAlloc(width, height, PS_TYPE_U8); - psImageInit(sf->warpedPixels, 0); - - for (int i=0; iinAstrom->fpa->toSky->type != PS_PROJ_DIS) { - pmAstromWriteWCS(header, sf->inAstrom->fpa, sf->chip, 0.001); - } - psFits *fits = psFitsOpen("warpedpixels.fits", "w"); - - psFitsWriteImage(fits, header, sf->warpedPixels, 0, NULL); - psFitsClose(fits); - psFree(header); - } - return true; -} - -static void -addTouchedPixels(streakFiles *sf, psString fileName) -{ - // read fits header - psFits *fits = psFitsOpen(fileName, "r"); - if (!fits) { - psError(PS_ERR_IO, false, "failed to open skycell file: %s", fileName); - streaksExit("", PS_EXIT_DATA_ERROR); - } - psMetadata *header = psFitsReadHeader(NULL, fits); - psFitsClose(fits); - if (!header) { - psError(PS_ERR_IO, false, "failed to read fixts header from skycell file: %s", fileName); - streaksExit("", PS_EXIT_DATA_ERROR); - } - // set up astrometry for conversion from skycell image to sky - pmAstromWCS *wcs = pmAstromWCSfromHeader(header); - if (!wcs) { - psError(PS_ERR_IO, false, "failed to read astrometry from skycell file: %s", fileName); - streaksExit("", PS_EXIT_DATA_ERROR); - } - int naxis1 = psMetadataLookupS32(NULL, header, "NAXIS1"); - int naxis2 = psMetadataLookupS32(NULL, header, "NAXIS2"); - - psFree(header); - - /* now set up our wrapper to the chip astrometry to apply to the whole chip */ - sf->astrom = streakSetAstrometry(sf->astrom, sf->stage, sf->inAstrom->fpa, sf->chip, false, NULL, - sf->warpedPixels->numCols, sf->warpedPixels->numRows); - - // convert corners of skycell to sky coordinates - // compute overlap with the chip using astrometry - psPlane pt[4]; - pt[0].x = 0; - pt[0].y = 0; - - pt[1].x = naxis1 - 1; - pt[1].y = 0; - - pt[2].x = naxis1 - 1; - pt[2].y = naxis2 - 1; - - pt[3].x = 0; - pt[3].y = naxis2 - 1; - - for (int i = 0; i < 4 ; i++ ) { - psSphere sky; - - // convert corner of skycell to sky coordinates - if (!pmAstromWCStoSky(&sky, wcs, &pt[i])) { - psErrorStackPrint(stderr, "failed to convert pt %d of %s to sky coords", i, fileName); - streaksExit("", PS_EXIT_DATA_ERROR); - } - strkPt p; - // convert to chip coordinates - if (!skyToCell(&p, sf->astrom, sky.r, sky.d)) { - psErrorStackPrint(stderr, "failed to convert pt %d of %s to sky coords", i, fileName); - streaksExit("", PS_EXIT_DATA_ERROR); - } - pt[i].x = p.x; - pt[i].y = p.y; - } - - psFree(wcs); - - // put the corners in the desired order (see comments below) - - nameCorners(pt); - psString type; - if (pt[0].y < pt[2].y ) { - type = "1"; - } else { - type = "2"; - } -#ifdef DEBUG_PRINT - printf("\nSKYCELL %s Type: %s\n", fileName, type); - for (int i=0; i<4; i++) { - printf("%f %f\n", pt[i].x, pt[i].y); - } -#endif - // Now set the touched pixels - int ymin = fmax(0, pt[1].y ); - int ymax = fmin(pt[3].y + 0.5, sf->warpedPixels->numRows - 1); -#if (DEBUG_PRINT > 1) - printf("\nymin: %d ymax: %d\n", ymin, ymax); -#endif - for (int y = ymin ; y <= ymax; y++) { - int xleft = xLeft(pt, y); - int xright = xRight(pt, y); - - if (xleft < 0) { - xleft = 0; - } - if (xleft > sf->warpedPixels->numCols) { - continue; - } - if (xright < 0) { - continue; - } - if (xright >= sf->warpedPixels->numCols) { - xright = sf->warpedPixels->numCols - 1; - } -#if (DEBUG_PRINT > 1) - printf(" y: %d xleft: %d xright: %d\n", y, xleft, xright); -#endif - - psU8 *rowPixels = sf->warpedPixels->data.U8[y]; - - int n = xright - xleft + 1; - if (n > 0) { - memset(&rowPixels[xleft], 255, n); - } - } -} - -// x as a function of y for the line between two points -// Note: the caller guarentees that the y's of the two points are different -static double xOfY(psPlane *pI, psPlane *pJ, int y) -{ - double m_inv = (double) (pJ->x - pI->x) / (double) (pJ->y - pI->y); - - double x = ((double) y - pI->y) * m_inv + (double) pI->x; - - return x; -} - -static int xLeft(psPlane *pt, int y) -{ - double x_d; - if (y < pt[0].y) { - // left boundary is the line from pt 0 to pt 1 - x_d = xOfY(&pt[0], &pt[1], y); - } else { - // left boundary is the line from pt 0 to pt 3 - x_d = xOfY(&pt[0], &pt[3], y); - } - return (int) (x_d + 0.5); -} - -static int xRight(psPlane *pt, int y) -{ - double x_d; - if (y < pt[2].y) { - // right boundary is the line from pt 1 to pt 2 - x_d = xOfY(&pt[1], &pt[2], y); - - } else { - // right boundary is the line from pt 2 to pt 3 - x_d = xOfY(&pt[2], &pt[3], y); - } - return (int) x_d; -} - -/* - * To compute the overlap of a quadrilateral with a set of - * points we transform the 4 corners to image coordinates - * and name the 4 points of the quad - pt 0 is left most (bottom most corner) - pt 1 is bottom most (right most corner) - pt 2 is right most (top most corner) - pt 3 is top most (left most corner) - -* Then we can divide the quad into 3 regions A B and C based -* on the y coordinate - -* (in the degenerate case of a rectangle aligned to the axes -* we only have 1 region) -* In each of the three regions A, B, and C the test for whether a -* point is contained in the quad on scanline y is simply - - x >= xleft(y) x < xright(y) - -* where xleft(y) and xright(y) are the bounding lines at the given -* y coordinate. - -The following diagrams show some examples. The 4 points are labeled 0 to 3 -The three regions are labeld A B and C, the boundaries between the regions -ocur at pt0.y and pt2.y - -Example 1 - - 3 - C - ---------------2 left boundary: line 0_1 y < pt0.y - B line 0_3 y >= pt0.y - 0-------------- - A right boundary: line 1_2 y < pt2.y - 1 line 2_3 y >= pt.y -************************** -Example 2 - 3 - C - 0---------------- - B - ----------------2 - A - 1 - -************************** -Example 3 - - 3 2 - - B - - 0 1 - -(region A and C are empty in the case where the points form a rectangle) -*/ - -/* - Name the corners - The following algorithm works for points that form a quadrilateral - I think it also works for the situation where 3 points are co-linear - and we have a triangle, but that isn't important for our purposes - -*/ - -static void -nameCorners(psPlane pt[4]) -{ - psPlane pt0, pt1, pt2, pt3; - - /* find pt0 the left most (bottom most) corner */ - int imin = 0; - int min = (int) pt[0].x; - int i; - for (i=1; i<4; i++) { - int x = pt[i].x; - if ((x < min) || - ((x == min) && (pt[i].y < pt[imin].y))) { - imin = i; - min = x; - } - } - // remve pt0 from the list - pt0 = pt[imin]; - for (i=imin; i<3; i++) { - pt[i] = pt[i+1]; - } - - // find the bottom most (right most) corner from the remaining 3 - imin = 0; - min = pt[0].y; - for (i=1; i<3; i++) { - int y = pt[i].y; - if ((y < min) || - ((y == min) && (pt[i].x > pt[imin].x)) ) { - imin = i; - min = y; - } - } - // remve pt1 from the list - pt1 = pt[imin]; - for (i=imin; i<3; i++) { - pt[i] = pt[i+1]; - } - - // now find the right most (top most) of the remaining 2 points - if ((pt[0].x > pt[1].x) || - ((pt[0].x == pt[1].x) && (pt[0].y > pt[1].y))) { - pt2 = pt[0]; - pt3 = pt[1]; - } else { - pt2 = pt[1]; - pt3 = pt[0]; - } - /* write the outputs */ - pt[0] = pt0; - pt[1] = pt1; - pt[2] = pt2; - pt[3] = pt3; -}