Index: branches/pap/magic/remove/src/streaksremove.c =================================================================== --- branches/pap/magic/remove/src/streaksremove.c (revision 23948) +++ branches/pap/magic/remove/src/streaksremove.c (revision 25027) @@ -1,2 +1,11 @@ +/* + * streaksremove + * + * Convert satellite streak detctions into masks and remove the covered pixels from the + * input images. + * Optionally swap the inputs with the outputs so that subsequent references to the original + * images access the destreaked versions. + */ + #include "streaksremove.h" @@ -4,8 +13,9 @@ static bool readAndCopyToOutput(streakFiles *sf, bool exciseAll); static void exciseNonWarpedPixels(streakFiles *sfiles, double newMaskValue); -static bool warpedPixel(streakFiles *sfiles, PixelPos *cellCoord); +static bool warpedPixel(streakFiles *sfiles, int x, int y); static void excisePixel(streakFiles *sfiles, unsigned int x, unsigned int y, bool streak, double newMaskValue); static void writeImages(streakFiles *sf, bool exciseImageCube); static void updateAstrometry(streakFiles *sfiles); +static void censorSources(streakFiles *sfiles, psU32 maskStreak); int @@ -23,18 +33,24 @@ } - psU32 maskStreak = 0; - psU32 maskMask = 0; + // 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"); + // call Paul Sydney's code to parse the streaks file that DetectStreaks produced Streaks *streaks = readStreaksFile(streaksFileName); if (!streaks) { - psErrorStackPrint(stderr, "failed to read streaks file: %s", streaksFileName); + psError(PS_ERR_UNKNOWN, "failed to read streaks file: %s", streaksFileName); streaksExit("", PS_EXIT_PROG_ERROR); } - streakFiles *sfiles = openFiles(config, true); + // open all of the input and output files, save their descriptions in the streakFiles struct + streakFiles *sfiles = openFiles(config, true, argv[0]); setupAstrometry(sfiles); + // Optionally we can set pixels that are masked to NAN since they couldn't have been + // examined for streaks. Usually this is done by the distribution system just prior + // to release bool nanForRelease = psMetadataLookupBool(&status, config->arguments, "NAN_FOR_RELEASE"); if (nanForRelease && (sfiles->inMask == NULL)) { @@ -52,5 +68,5 @@ if (checkNonWarpedPixels ) { - // From ICD: + // From magic ICD: // In the raw and detrended images, the pixels which were not // included in any of the streak-processed warps must also be masked. @@ -65,9 +81,11 @@ } psF64 cwp_t = psTimerClear("COMPUTE_WARPED_PIXELS"); - printf("time to compute warped pixels: %f\n", cwp_t); + psLogMsg("streaksremove", PS_LOG_INFO, "time to compute warped pixels: %f\n", cwp_t); } if (sfiles->stage == IPP_STAGE_RAW) { - // copy PHU to output files + // Except for raw stage, all of our (GPC1) files have one image extension. + // Raw files have a phu and multiple extensions, one per chip + // Since this is a raw file, copy it's PHU to output files copyPHU(sfiles, true); @@ -82,5 +100,5 @@ int totalStreakPixels = 0; - // Iterate through each component of the input (there is only one except for raw images) + // Iterate through each component of the input (except for raw images there is only one) do { bool exciseImageCube = false; @@ -110,10 +128,11 @@ psTimerStart("GET_STREAK_PIXELS"); - StreakPixels *pixels = streak_on_component (streaks, sfiles->astrom, - sfiles->inImage->numCols, sfiles->inImage->numRows); - - printf("time to get streak pixels: %f\n", psTimerClear("GET_STREAK_PIXELS")); - + // call Paul Sydney's code to compute the set of pixels that are covered by the detected streaks + StreakPixels *pixels = streak_on_component(streaks, sfiles->astrom, sfiles->inImage->numCols, + sfiles->inImage->numRows); + psLogMsg("streaksremove", PS_LOG_INFO, "time to get streak pixels: %f\n", psTimerClear("GET_STREAK_PIXELS")); + // if this extension contained an image, excise the streaked pixels. + // otherwise it contained an image cube (video cell) which is handled in the if block if (sfiles->inImage->image) { if (checkNonWarpedPixels) { @@ -124,91 +143,103 @@ exciseNonWarpedPixels(sfiles, maskStreak); - printf("time to excise non warped pixels: %f\n", psTimerClear("EXCISE_NON_WARPED")); + psLogMsg("streaksremove", PS_LOG_INFO, "time to excise non warped pixels: %f\n", psTimerClear("EXCISE_NON_WARPED")); } - totalStreakPixels += psArrayLength(pixels); + + psTimerStart("REMOVE_STREAKS"); - for (int i = 0; i < psArrayLength (pixels); ++i) { - PixelPos *pixelPos = psArrayGet (pixels, i); - - if (!checkNonWarpedPixels || warpedPixel(sfiles, pixelPos)) { - - excisePixel(sfiles, pixelPos->x, pixelPos->y, true, maskStreak); - - } else { - // This pixel was not included in any warp and has thus already excised - // by exciseNonWarpedPixels + + for (int y=0 ; y < sfiles->inImage->numRows; y++) { + for (int x = 0; x < sfiles->inImage->numCols; x++) { + if (psImageGet(pixels, x, y)) { + ++totalStreakPixels; + if (!checkNonWarpedPixels || warpedPixel(sfiles, x, y)) { + + excisePixel(sfiles, x, y, true, maskStreak); + + } else { + // This pixel was not included in any warp and has thus already excised + // by exciseNonWarpedPixels + } + } } } - printf("time to remove streak pixels: %f\n", psTimerClear("REMOVE_STREAKS")); + + psLogMsg("streaksremove", PS_LOG_INFO, "time to remove streak pixels: %f\n", psTimerClear("REMOVE_STREAKS")); if (nanForRelease) { + // set any pixels that were masked, to NAN (unless they are already NAN) setMaskedToNAN(sfiles, maskMask, true); } } else { - // this component contains an image cube, excise it completely + // this component contains an image cube + // For now excise it completely exciseImageCube = true; } - psArrayElementsFree (pixels); psFree(pixels); } if (sfiles->stage == IPP_STAGE_CHIP) { + // as a convience to the user of the output, replace the bogus WCS transform in the + // chip processed files with the data calcuated by psastro at the camera stage + // (actually we use a linear approximation) updateAstrometry(sfiles); } - // write out the destreaked temporary images and the recovery images + censorSources(sfiles, maskStreak); + + // write the destreaked "temporary" images and the recovery images writeImages(sfiles, exciseImageCube); - printf("time to process component %d: %f\n", sfiles->extnum, psTimerClear("PROCESS_COMPONENT")); + + psLogMsg("streaksremove", PS_LOG_INFO, "time to process component %d: %f\n", sfiles->extnum, psTimerClear("PROCESS_COMPONENT")); } while (streakFilesNextExtension(sfiles)); + psFree(streaks); - printf("pixels: %ld streak pixels: %ld %4.2f%%\n", totalPixels, totalStreakPixels, 100. * totalStreakPixels / totalPixels); + psLogMsg("streaksremove", PS_LOG_INFO, "pixels: %ld streak pixels: %ld %4.2f%%\n", totalPixels, totalStreakPixels, 100. * totalStreakPixels / totalPixels); + + // check the weight and mask files for extra extensions that might be in files + // (covariance matrix for example) + copyExtraExtensions(sfiles); + + // all done close the files. This is where the files are written so it can take a long time. psTimerStart("CLOSE_IMAGES"); - // close all files + closeImages(sfiles); - printf("time to close images: %f\n", psTimerClear("CLOSE_IMAGES")); + + psLogMsg("streaksremove", PS_LOG_INFO, "time to close images: %f\n", psTimerClear("CLOSE_IMAGES")); + + + if (!replicateOutputs(sfiles)) { + psError(PS_ERR_UNKNOWN, false, "failed to replicate output files"); + deleteTemps(sfiles); + psErrorStackPrint(stderr, ""); + exit(PS_EXIT_UNKNOWN_ERROR); + } // NOTE: from here on we can't just quit if something goes wrong. // especially if we're working at the raw stage - - if (!replicateOutputs(sfiles)) { - psError(PS_ERR_UNKNOWN, false, "failed to replicate output files"); - psErrorStackPrint(stderr, ""); - exit(PS_EXIT_UNKNOWN_ERROR); - } + // turn off automatic deletion of output files by streaksExit + setStreakFiles(NULL); if (psMetadataLookupBool(&status, config->arguments, "REPLACE")) { // swap the instances for the input and output - // Note this is a database operation. No file I/O is performed + // Note this is a nebulous database operation. No file I/O is performed if (!swapOutputsToInputs(sfiles)) { - psError(PS_ERR_UNKNOWN, false, "failed to swap files"); - - // XXX: Now what? I guess swapOutputsToInputs will need to undo anything that - // it has done and give a detailed report of what happened - - psErrorStackPrint(stderr, ""); + // 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 + + psErrorStackPrint(stderr, "failed to swap files"); + + // XXX: pick a specific error code for this failure exit(PS_EXIT_UNKNOWN_ERROR); } - - if (psMetadataLookupBool(&status, config->arguments, "REMOVE")) { - // delete the temporary storage objects (which now points to the original image(s) - if (!deleteTemps(sfiles)) { - psError(PS_ERR_UNKNOWN, false, "failed to delete temporary files"); - // XXX: Now what? At this point the output files have been swapped, so we can't - // repeat the operation. - - // Returning error status here is problematic. The inputs have been streak removed - // but they're still lying around - // Maybe just print an error message and - // let other system tools clean up - psErrorStackPrint(stderr, ""); - exit(PS_EXIT_UNKNOWN_ERROR); - } - } - } + } + // all done. Clean up to look for memory leaks. psFree(sfiles); @@ -218,5 +249,5 @@ streaksNebulousCleanup(); pmConfigDone(); - printf("time to run streaksremove: %f\n", psTimerClear("STREAKSREMOVE")); + psLogMsg("streaksremove", PS_LOG_INFO, "time to run streaksremove: %f\n", psTimerClear("STREAKSREMOVE")); psLibFinalize(); @@ -235,5 +266,4 @@ 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-remove: remove the original image after processing (requires -replace)\n"); fprintf(stderr, "\t-keepnonwarped: 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"); @@ -311,4 +341,8 @@ bool gotReplace = false; if ((argnum = psArgumentGet(argc, argv, "-replace"))) { + if (!nebulousImage) { + psError(PS_ERR_UNKNOWN, true, "replace is only supported for nebulous files"); + usage(); + } gotReplace = true; psArgumentRemove(argnum, &argc, argv); @@ -394,4 +428,16 @@ psArgumentRemove(argnum, &argc, argv); } + if ((argnum = psArgumentGet(argc, argv, "-sources"))) { + psArgumentRemove(argnum, &argc, argv); + bool nebulousSources = IN_NEBULOUS(argv[argnum]); + if (nebulousSources != nebulousImage) { + psError(PS_ERR_UNKNOWN, true, "sources file must have %snebulous path with %s image path\n", + nebulousImage ? "" : "non ", nebulousImage ? "nebulous" : "regular"); + usage(); + } + psMetadataAddStr(config->arguments, PS_LIST_TAIL, "SOURCES", 0, + "name of input sources file", argv[argnum]); + psArgumentRemove(argnum, &argc, argv); + } if ((argnum = psArgumentGet(argc, argv, "-tmproot"))) { @@ -403,8 +449,6 @@ usage(); } - psString dir = pathToDirectory(argv[argnum]); - psMetadataAddStr(config->arguments, PS_LIST_TAIL, "OUTPUT", 0, "directory for temporary files", - dir); - psFree(dir); + psMetadataAddStr(config->arguments, PS_LIST_TAIL, "OUTPUT", 0, "path for (temporary if replae) output files", + argv[argnum]); psArgumentRemove(argnum, &argc, argv); } else { @@ -415,22 +459,10 @@ if ((argnum = psArgumentGet(argc, argv, "-recovery"))) { psArgumentRemove(argnum, &argc, argv); - psString dir = pathToDirectory(argv[argnum]); - psMetadataAddStr(config->arguments, PS_LIST_TAIL, "RECOVERY", 0, "directory for recovery files", - dir); - psFree(dir); + psMetadataAddStr(config->arguments, PS_LIST_TAIL, "RECOVERY", 0, "path for recovery files", + argv[argnum]); psArgumentRemove(argnum, &argc, argv); } else if ((stage == IPP_STAGE_RAW) && gotReplace) { psError(PS_ERR_UNKNOWN, true, "-recovery is required for -stage raw with -replace\n"); usage(); - } - - if ((argnum = psArgumentGet(argc, argv, "-remove"))) { - if (!gotReplace) { - psError(PS_ERR_UNKNOWN, true, "-replace is required with -remove\n"); - usage(); - } - psArgumentRemove(argnum, &argc, argv); - psMetadataAddBool(config->arguments, PS_LIST_TAIL, "REMOVE", 0, "remove original files", - true); } @@ -451,7 +483,11 @@ updateAstrometry(streakFiles *sf) { + // XXX: why do I check this here? Shouldn't it be just around the call to linearizeTransforms? if (sf->bilevelAstrometry) { - linearizeTransforms(sf->astrom); + if (!linearizeTransforms(sf->astrom)) { + // fit failed, leave the astrometry unchanged + return; + } if (!pmAstromWriteWCS(sf->outImage->header, sf->inAstrom->fpa, sf->chip, 0.001)) { @@ -501,8 +537,10 @@ } } - sf->outImage->header = (psMetadata*) psMemIncrRefCounter(sf->inImage->header); + sf->outImage->header = psMemIncrRefCounter(sf->inImage->header); if (sf->recImage) { - sf->recImage->header = (psMetadata*) psMemIncrRefCounter(sf->inImage->header); - } + sf->recImage->header = psMetadataCopy(NULL, sf->inImage->header); + addRecoveryKeyword(sf->recImage->header); + } + addDestreakKeyword(sf->outImage->header); if (!SFILE_IS_IMAGE(sf->inImage)) { @@ -520,17 +558,5 @@ // set up the compression parameters -#ifdef STREAKS_COMPRESS_OUTPUT - // compression of the image pixels is disabled for now. Some consortium members - // have problems reading them - copyFitsOptions(sf->outImage, sf->recImage, sf->inImage); - - // XXX: TODO: can we derive these values from the input header? - // psFitsCompressionGet(sf->inImage->image) gives compression none - // perhaps we should just use the definition of COMP_IMG in the configuration - psFitsSetCompression(sf->outImage->fits, PS_FITS_COMPRESS_RICE, sf->tiles, 8, 0, 0); - if (sf->recImage) { - psFitsSetCompression(sf->recImage->fits, PS_FITS_COMPRESS_RICE, sf->tiles, 8, 0, 0); - } -#endif + copyFitsOptions(sf->outImage, sf->recImage, sf->inImage, sf->tiles); if (sf->inMask) { @@ -539,6 +565,8 @@ sf->outMask->header = (psMetadata*) psMemIncrRefCounter(sf->inMask->header); if (sf->recMask) { - sf->recMask->header = (psMetadata*) psMemIncrRefCounter(sf->inMask->header); - } + sf->recMask->header = psMetadataCopy(NULL, sf->outMask->header); + addRecoveryKeyword(sf->recMask->header); + } + addDestreakKeyword(sf->outMask->header); if (updateAstrometry) { pmAstromWriteWCS(sf->outMask->header, sf->inAstrom->fpa, sf->chip, 0.001); @@ -554,19 +582,8 @@ } -#ifdef STREAKS_COMPRESS_OUTPUT - // XXX: see note above - copyFitsOptions(sf->outMask, sf->recMask, sf->inMask); - psFitsSetCompression(sf->outMask->fits, PS_FITS_COMPRESS_PLIO, sf->tiles, 8, 0, 0); - if (sf->recMask) { - psFitsSetCompression(sf->recMask->fits, PS_FITS_COMPRESS_PLIO, sf->tiles, 8, 0, 0); - } + copyFitsOptions(sf->outMask, sf->recMask, sf->inMask, sf->tiles); if (sf->outChMask) { - copyFitsOptions(sf->outChMask, sf->recChMask, sf->inMask); - psFitsSetCompression(sf->outChMask->fits, PS_FITS_COMPRESS_PLIO, sf->tiles, 8, 0, 0); - if (sf->recChMask) { - psFitsSetCompression(sf->recChMask->fits, PS_FITS_COMPRESS_PLIO, sf->tiles, 8, 0, 0); - } - } -#endif + copyFitsOptions(sf->outChMask, sf->recChMask, sf->inMask, sf->tiles); + } } } @@ -576,6 +593,8 @@ sf->outWeight->header = (psMetadata*) psMemIncrRefCounter(sf->inWeight->header); if (sf->recWeight) { - sf->recWeight->header = (psMetadata*) psMemIncrRefCounter(sf->inWeight->header); - } + sf->recWeight->header = psMetadataCopy(NULL, sf->outWeight->header); + addRecoveryKeyword(sf->recWeight->header); + } + addDestreakKeyword(sf->outWeight->header); if (updateAstrometry) { pmAstromWriteWCS(sf->inWeight->header, sf->inAstrom->fpa, sf->chip, 0.001); @@ -583,20 +602,9 @@ setupImageRefs(sf->outWeight, sf->recWeight, sf->inWeight, sf->extnum, exciseAll); -#ifdef STREAKS_COMPRESS_OUTPUT - copyFitsOptions(sf->outWeight, sf->recWeight, sf->inWeight); - // XXX: see note above - psFitsSetCompression(sf->outWeight->fits, PS_FITS_COMPRESS_RICE, sf->tiles, 8, 0, 0); - if (sf->recWeight) { - psFitsSetCompression(sf->recWeight->fits, PS_FITS_COMPRESS_RICE, sf->tiles, 8, 0, 0); - } -#endif - } - // and for raw images, create sub images that represent the actual image - // area (no overscan) + copyFitsOptions(sf->outWeight, sf->recWeight, sf->inWeight, sf->tiles); + } return true; } - - static void @@ -746,5 +754,5 @@ static bool -warpedPixel(streakFiles *sfiles, PixelPos *cellCoord) +warpedPixel(streakFiles *sfiles, int x, int y) { PixelPos chipCoord; @@ -757,10 +765,10 @@ // we clip so that the streak calculation code doesn't have to // clipping here insures that we don't touch the overscan regions - if ((cellCoord->x < 0) || (cellCoord->x >= sfiles->inImage->numCols) || - (cellCoord->y < 0) || (cellCoord->y >= sfiles->inImage->numRows)) { + if ((x < 0) || (x >= sfiles->inImage->numCols) || + (y < 0) || (y >= sfiles->inImage->numRows)) { return false; } - cellToChipInt(&chipCoord.x, &chipCoord.y, sfiles->astrom, cellCoord->x, cellCoord->y); + cellToChipInt(&chipCoord.x, &chipCoord.y, sfiles->astrom, x, y); if (chipCoord.x < 0 || chipCoord.x >= sfiles->warpedPixels->numCols) { @@ -774,2 +782,80 @@ } +// read a sources file (.cmf) and remove any rows whose coordinate is convered by a +// streak mask +static void +censorSources(streakFiles *sfiles, psU32 maskStreak) +{ + if ((!sfiles->inSources) || (!sfiles->outMask)) { + return; + } + psImage *maskImage = sfiles->outMask->image; + if (!maskImage) { + psError(PS_ERR_IO, false, "maskImage is null!"); + streaksExit("", PS_EXIT_PROG_ERROR); + } + + sFile *in = sfiles->inSources; + sFile *out = sfiles->outSources; + + in->header = psFitsReadHeader(NULL, in->fits); + if (!in->header) { + psError(PS_ERR_IO, false, "failed to read header from %s", in->resolved_name); + streaksExit("", PS_EXIT_DATA_ERROR); + } + + bool status; + psString extname = psMetadataLookupStr(&status, in->header, "EXTNAME"); + if (!extname) { + psError(PS_ERR_IO, false, "failed to find extname in header of %s", in->resolved_name); + streaksExit("", PS_EXIT_DATA_ERROR); + } + + psArray *inTable = psFitsReadTable(in->fits); + if (!inTable->n) { + psError(PS_ERR_IO, false, "table in %s is empty", in->resolved_name); + streaksExit("", PS_EXIT_DATA_ERROR); + } + + psArray *outTable = psArrayAllocEmpty(inTable->n); + int j = 0; + int numCensored = 0; + for (int i = 0 ; i < inTable->n; i++) { + psMetadata *row = inTable->data[i]; + + psF32 x = psMetadataLookupF32 (&status, row, "X_PSF"); + psF32 y = psMetadataLookupF32 (&status, row, "Y_PSF"); + + psU32 mask = psImageGet(maskImage, x, y); + + // Key the source if the center pixel is not masked with maskStreak + if (! (mask & maskStreak) ) { + psArraySet(outTable, j++, row); + } else { + numCensored++; + } + } + + // get rid of unused elements (don't know if this is necessary) + psArrayRealloc(outTable, j); + + addDestreakKeyword(in->header); + if (psArrayLength(outTable) > 0) { + printf("Censored %d sources\n", numCensored); + if (! psFitsWriteTable(out->fits, in->header, outTable, extname)) { + psError(PS_ERR_IO, false, "failed to write table to %s", out->resolved_name); + streaksExit("", PS_EXIT_DATA_ERROR); + } + } else { + printf("Censored ALL %d sources\n", numCensored); + if (! psFitsWriteTableEmpty(out->fits, in->header, inTable->data[0], extname)) { + psError(PS_ERR_IO, false, "failed to write empty table to %s", out->resolved_name); + streaksExit("", PS_EXIT_DATA_ERROR); + } + } + + if (!psFitsClose(out->fits)) { + psError(PS_ERR_IO, false, "failed to close table %s", out->resolved_name); + streaksExit("", PS_EXIT_DATA_ERROR); + } +}