Index: /trunk/magic/remove/src/streaksremove.c =================================================================== --- /trunk/magic/remove/src/streaksremove.c (revision 20779) +++ /trunk/magic/remove/src/streaksremove.c (revision 20780) @@ -290,5 +290,5 @@ char *inputBasename = basename(sf->inImage->name); sf->outImage = sFileOpen(config, stage, "OUTPUT", inputBasename, true); - sf->recImage = sFileOpen(config, stage, "RECOVERY", inputBasename, true); + sf->recImage = sFileOpen(config, stage, "RECOVERY", inputBasename, false); sf->inMask = sFileOpen(config, stage, "INPUT.MASK", NULL, false); @@ -490,11 +490,12 @@ } - - if (!pmConfigFileSetsMD(config->arguments, &argc, argv, "SKYCELLS", "-skycell", "-skycelllist")) { ; - if (CHIP_LEVEL_INPUT(stage)) { - psError(PS_ERR_UNKNOWN, true, "-skycelllist is required for raw and chip stages\n"); - return NULL; - } - } + 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 skycells are not provided then we have to execise all pixels unless -keepnonwarped + pmConfigFileSetsMD(config->arguments, &argc, argv, "SKYCELLS", "-skycell", "-skycelllist"); if (!pmConfigFileSetsMD(config->arguments, &argc, argv, "ASTROM", "-astrom", NULL)) { ; @@ -504,4 +505,5 @@ } } + if ((argnum = psArgumentGet(argc, argv, "-mask"))) { psArgumentRemove(argnum, &argc, argv); @@ -516,4 +518,5 @@ psArgumentRemove(argnum, &argc, argv); } + if ((argnum = psArgumentGet(argc, argv, "-weight"))) { psArgumentRemove(argnum, &argc, argv); @@ -528,4 +531,5 @@ psArgumentRemove(argnum, &argc, argv); } + if ((argnum = psArgumentGet(argc, argv, "-tmproot"))) { psArgumentRemove(argnum, &argc, argv); @@ -544,4 +548,5 @@ return NULL; } + if ((argnum = psArgumentGet(argc, argv, "-recovery"))) { psArgumentRemove(argnum, &argc, argv); @@ -550,8 +555,9 @@ dir); psArgumentRemove(argnum, &argc, argv); - } else if (stage == IPP_STAGE_RAW) { - psError(PS_ERR_UNKNOWN, true, "-recovery is required for -stage raw\n"); + } else if ((stage == IPP_STAGE_RAW) && gotReplace) { + psError(PS_ERR_UNKNOWN, true, "-recovery is required for -stage raw with -replace\n"); return NULL; } + if ((argnum = psArgumentGet(argc, argv, "-remove"))) { if (!gotReplace) { @@ -589,5 +595,4 @@ } - // TODO: add keyword indicating that streaks have been removed if (!psFitsWriteBlank(sfiles->outImage->fits, imageHeader, NULL)) { @@ -597,5 +602,5 @@ } // TODO: add keyword indicating that this is the recovery image - if (!psFitsWriteBlank(sfiles->recImage->fits, imageHeader, NULL)) { + if (sfiles->recImage && !psFitsWriteBlank(sfiles->recImage->fits, imageHeader, NULL)) { psError(PS_ERR_IO, false, "failed to write primary header to %s", sfiles->recImage->resolved_name); @@ -619,5 +624,5 @@ } // TODO: add keyword indicating that this is the recovery image - if (!psFitsWriteBlank(sfiles->recMask->fits, maskHeader, NULL)) { + if (sfiles->recMask && !psFitsWriteBlank(sfiles->recMask->fits, maskHeader, NULL)) { psError(PS_ERR_IO, false, "failed to write primary header to %s", sfiles->recMask->resolved_name); @@ -641,5 +646,5 @@ } // TODO: add keyword indicating that this is a recovery image - if (!psFitsWriteBlank(sfiles->recWeight->fits, weightHeader, NULL)) { + if (sfiles->recWeight && !psFitsWriteBlank(sfiles->recWeight->fits, weightHeader, NULL)) { psError(PS_ERR_IO, false, "failed to write primary header to %s", sfiles->recWeight->resolved_name); @@ -764,5 +769,4 @@ streaksremoveExit("", PS_EXIT_DATA_ERROR); } - in->numRows = in->image->numRows; } else { if (stage != IPP_STAGE_RAW) { @@ -784,4 +788,7 @@ setFitsOptions(sFile *sfile, int bitpix, float bscale, float bzero) { + if (!sfile) { + return; + } if (sfile->fits->options) { psFree(sfile->fits->options); @@ -822,6 +829,8 @@ } -// setFitsOptions(out, bitpix, bscale, bzero); -// setFitsOptions(rec, bitpix, bscale, bzero); +#ifdef COMPRESS_OUTPUT + setFitsOptions(out, bitpix, bscale, bzero); + setFitsOptions(rec, bitpix, bscale, bzero); +#endif } @@ -859,17 +868,37 @@ static void -createRecoveryImage(sFile *rec, sFile *in, int extnum) -{ - rec->image = psImageRecycle(rec->image, in->image->numCols, in->image->numRows, in->image->type.type); - if (!rec->image) { - psError(PS_ERR_UNKNOWN, false, "failed to allocate recoveyr image for extnsion %d", extnum); +createNewImage(sFile *out, sFile *in, int extnum, double initValue) +{ + out->image = psImageRecycle(out->image, in->image->numCols, in->image->numRows, in->image->type.type); + if (!out->image) { + psError(PS_ERR_UNKNOWN, false, "failed to allocate image for extnsion %d", extnum); streaksremoveExit("", PS_EXIT_UNKNOWN_ERROR); } - psImageInit(rec->image, 0); + psImageInit(out->image, initValue); +} + +static void +setupImageRefs(sFile *out, sFile *recovery, sFile *in, int extnum, bool exciseAll) +{ + if (exciseAll) { + // set the output to an image all NANs + createNewImage(out, in, extnum, NAN); + if (recovery) { + // save the whole input as the recovery image + recovery->image = (psImage*) psMemIncrRefCounter(in->image); + } + } else { + // output is the input image + out->image = (psImage*) psMemIncrRefCounter(in->image); + if (recovery) { + // create a recovery image initialized to zero + createNewImage(recovery, in, extnum, 0); + } + } } // set the image for the output files to the contents of the corresponding input file. static bool -readAndCopyToOutput(streakFiles *sf) +readAndCopyToOutput(streakFiles *sf, bool exciseAll) { bool updateAstrometry = false; @@ -893,5 +922,5 @@ if (sf->stage == IPP_STAGE_CHIP) { // For the chip level files, copy the WCS from the astrometry file to the header - updateAstrometry = true; + // 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); @@ -900,5 +929,7 @@ } sf->outImage->header = (psMetadata*) psMemIncrRefCounter(sf->inImage->header); - sf->recImage->header = (psMetadata*) psMemIncrRefCounter(sf->inImage->header); + if (sf->recImage) { + sf->recImage->header = (psMetadata*) psMemIncrRefCounter(sf->inImage->header); + } if (!SFILE_IS_IMAGE(sf->inImage)) { @@ -910,33 +941,38 @@ if (sf->inImage->image) { - sf->outImage->image = (psImage*) psMemIncrRefCounter(sf->inImage->image); - createRecoveryImage(sf->recImage, sf->inImage, sf->extnum); + setupImageRefs(sf->outImage, sf->recImage, sf->inImage, sf->extnum, exciseAll); } else { - // we don't do a reference if we have an imagecube + // Image cubes are handeled specially } // set up the compression parameters - // copyFitsOptions(sf->outImage, sf->recImage, sf->inImage); +#ifdef 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); - - psFitsSetCompression(sf->recImage->fits, PS_FITS_COMPRESS_RICE, sf->tiles, 8, 0, 0); + psFitsSetCompression(sf->outImage->fits, PS_FITS_COMPRESS_RICE, sf->tiles, 8, 0, 0); +#endif + if (sf->recImage) { + psFitsSetCompression(sf->recImage->fits, PS_FITS_COMPRESS_RICE, sf->tiles, 8, 0, 0); + } if (sf->inMask) { readImage(sf->inMask, sf->extnum, sf->stage); - sf->outMask->header = (psMetadata*) psMemIncrRefCounter(sf->inMask->header); if (updateAstrometry) { pmAstromWriteWCS(sf->outMask->header, sf->inAstrom->fpa, sf->chip, 0.001); } - sf->outMask->image = (psImage*) psMemIncrRefCounter(sf->inMask->image); - sf->recMask->image = (psImage*) psMemIncrRefCounter(sf->inMask->image); - - createRecoveryImage(sf->recMask, sf->inMask, sf->extnum); + sf->outMask->header = (psMetadata*) psMemIncrRefCounter(sf->inMask->header); + if (sf->recMask) { + sf->recMask->header = (psMetadata*) psMemIncrRefCounter(sf->inMask->header); + } + + setupImageRefs(sf->outMask, sf->recMask, sf->inMask, sf->extnum, exciseAll); + + // XXX: see note above copyFitsOptions(sf->outMask, sf->recMask, sf->inMask); - - // XXX: see note above psFitsSetCompression(sf->outMask->fits, PS_FITS_COMPRESS_PLIO, sf->tiles, 8, 0, 0); psFitsSetCompression(sf->recMask->fits, PS_FITS_COMPRESS_PLIO, sf->tiles, 8, 0, 0); @@ -945,12 +981,12 @@ if (sf->inWeight) { readImage(sf->inWeight, sf->extnum, sf->stage); + if (updateAstrometry) { + pmAstromWriteWCS(sf->inWeight->header, sf->inAstrom->fpa, sf->chip, 0.001); + } sf->outWeight->header = (psMetadata*) psMemIncrRefCounter(sf->inWeight->header); - if (updateAstrometry) { - pmAstromWriteWCS(sf->outWeight->header, sf->inAstrom->fpa, sf->chip, 0.001); - } - sf->outWeight->image = (psImage*) psMemIncrRefCounter(sf->inWeight->image); - - sf->recWeight->header = (psMetadata*) psMemIncrRefCounter(sf->inWeight->header); - createRecoveryImage(sf->recWeight, sf->inWeight, sf->extnum); + if (sf->recWeight) { + sf->recWeight->header = (psMetadata*) psMemIncrRefCounter(sf->inWeight->header); + } + setupImageRefs(sf->outMask, sf->recMask, sf->inMask, sf->extnum, exciseAll); copyFitsOptions(sf->outWeight, sf->recWeight, sf->inWeight); @@ -970,5 +1006,5 @@ writeImage(sFile *sfile, psString extname, int extnum) { - if (!sfile->image) { + if (!sfile || !sfile->image) { return; } @@ -987,4 +1023,7 @@ writeImageCube(sFile *sfile, psArray *imagecube, psString extname, int extnum) { + if (!imagecube) { + return; + } if (!psFitsWriteImageCube(sfile->fits, sfile->header, imagecube, extname)) { psError(PS_ERR_IO, false, "failed to write image to %s extnum: %d", @@ -997,5 +1036,5 @@ static void -writeImages(streakFiles *sf, bool exciseAll) +writeImages(streakFiles *sf, bool exciseImageCube) { psString extname = NULL; @@ -1005,28 +1044,35 @@ } if (sf->inImage->numZPlanes == 0) { + // note exciseing complete images is handled in readAndCopyToOutput writeImage(sf->outImage, extname, sf->extnum); writeImage(sf->recImage, extname, sf->extnum); } else { + // we have an image cube double initValue; - if (exciseAll) { + 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 NAN + // 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); - if (exciseAll) { + if (exciseImageCube) { sf->outImage->image = image; writeImage(sf->outImage, extname, sf->extnum); } else { - sf->recImage->image = image; - writeImage(sf->recImage, extname, sf->extnum); - } - // Assumption there are no mask and weight images for video cells + // write zero valued image to reccovery + if (sf->recImage) { + sf->recImage->image = image; + writeImage(sf->recImage, extname, sf->extnum); + } + } + // Assumption: there are no mask and weight images for video cells return; } @@ -1044,8 +1090,8 @@ closeImage(sFile *sfile) { - if (!sfile->fits) { + if (!sfile) { return; } - if (!psFitsClose(sfile->fits)) { + if (sfile->fits && !psFitsClose(sfile->fits)) { psError(PS_ERR_IO, false, "failed to close image to %s", sfile->resolved_name); @@ -1253,5 +1299,7 @@ double imageValue = psImageGet (sfiles->inImage->image, x, y); - psImageSet (sfiles->recImage->image, x, y, imageValue); + if (sfiles->recImage) { + psImageSet (sfiles->recImage->image, x, y, imageValue); + } #define ACTUALLY_MASK_PIXEL @@ -1407,4 +1455,31 @@ } +static StreakPixels * +getStreakPixels(streakFiles *sfiles, Streaks *streaks) +{ + StreakPixels *pixels; +#ifdef notyet + Streaks *componentStreaks = NULL; + Streaks *streaksToProcess; + + bool only_component_streaks = false; // XXX flesh this out and make it an option + if (only_component_streaks) { + componentStreaks = restrictStreaksToComponent(sfiles, streaks); + streaksToProcess = componentStreaks; + } else { + streaksToProcess = streaks; + } + + pixels = streak_on_component (streaksToProcess, sfiles->astrom, + sfiles->inImage->numCols, sfiles->inImage->numRows); + if (componentStreaks) { + psFree(componentStreaks); + } +#else + pixels = streak_on_component (streaks, sfiles->astrom, + sfiles->inImage->numCols, sfiles->inImage->numRows); +#endif + return pixels; +} int main(int argc, char *argv[]) @@ -1433,7 +1508,14 @@ streakFiles *sfiles = openFiles(config); + bool exciseAll = false; + // --keepnonwarped is a test and debug mode bool status; - bool checkWarpedPixels = ! psMetadataLookupBool(&status, config->arguments, "INCLUDE_NON_WARPED"); - if (checkWarpedPixels) { + 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 ) { // From ICD: // In the raw and detrended images, the pixels which were not @@ -1441,6 +1523,10 @@ // Note that the warp and difference skycells are only generated if more // than a small fraction of the pixels are lit by the input image. - - computeWarpedPixels(sfiles); + // if no skycells are provided sfiles->exciseAll is set to true + + if (! computeWarpedPixels(sfiles) ) { + // we have no choice to excise all pixels + exciseAll = true; + } } @@ -1457,71 +1543,54 @@ } - // Iterate through components of image + // Iterate through components of input files do { - bool exciseAll = false; - - // read the images and copy the pixels image from the inputs to the outputs - // This also sets up the astrometry and initializes the recovery images - if (!readAndCopyToOutput(sfiles)) { - // this extension is not an image skip (video descriptor?) - // XXX is there anything else that could be in an input file that we need to handle? + bool exciseImageCube = false; + + // read the images and copy the data from the inputs to the outputs + // This also sets up the astrometry and initializes the recovery image(s) + + if (!readAndCopyToOutput(sfiles, exciseAll)) { + // false return value indicates that this extension is not an image and + // the contents has already been copied to the output file. + // we don't need to process this extesion for streaks. continue; } - // Identify pixels to mask because of streaks - -#ifdef notyet - Streaks *componentStreaks = NULL; - Streaks *streaksToProcess; - - bool only_component_streaks = false; // XXX flesh this out and make it an option - if (only_component_streaks) { - componentStreaks = restrictStreaksToComponent(sfiles, streaks); - streaksToProcess = componentStreaks; - } else { - streaksToProcess = streaks; - } - - pixels = streak_on_component (streaksToProcess, sfiles->astrom, - sfiles->inImage->numCols, sfiles->inImage->numRows); - if (componentStreaks) { - psFree(componentStreaks); - } -#else - pixels = streak_on_component (streaks, sfiles->astrom, - sfiles->inImage->numCols, sfiles->inImage->numRows); -#endif - // XXX: - // - // PFS: Need to get mask weight for PM_MASK_DETECTOR. Is this stored - // in a config somewhere? Not sure how to properly set the maskValue. - // WES just pick some numbers - - if (checkWarpedPixels) { - exciseNonWarpedPixels(sfiles, MASK_BAD_WARP); - } - - if (sfiles->inImage->image) { - for (i = 0; i < psArrayLength (pixels); ++i) { - pixelPos = psArrayGet (pixels, i); - - if (!checkWarpedPixels || warpedPixel(sfiles, pixelPos)) { - - excisePixel(sfiles, pixelPos->x, pixelPos->y, MASK_STREAK); - - } else { - // This pixel was not included in any warp and has thus already excised above + if (!exciseAll) { + // Identify pixels to mask because of streaks + + pixels = getStreakPixels(sfiles, streaks); + + // XXX: + // + // PFS: Need to get mask weight for PM_MASK_DETECTOR. Is this stored + // in a config somewhere? Not sure how to properly set the maskValue. + + + if (checkNonWarpedPixels) { + exciseNonWarpedPixels(sfiles, MASK_BAD_WARP); + } + + if (sfiles->inImage->image) { + for (i = 0; i < psArrayLength (pixels); ++i) { + pixelPos = psArrayGet (pixels, i); + + if (!checkNonWarpedPixels || warpedPixel(sfiles, pixelPos)) { + + excisePixel(sfiles, pixelPos->x, pixelPos->y, MASK_STREAK); + + } else { + // This pixel was not included in any warp and has thus already excised above + } } + } else { + // this component contains an image cube, excise it completely + exciseImageCube = true; } - } else { - // This video cell is touched by a streak excise - // For now we excise all video cells. - // if (psArrayLength(pixels)) - exciseAll = true; - } - psArrayElementsFree (pixels); + psArrayElementsFree (pixels); + } // write out the destreaked temporary images and the recovery images - writeImages(sfiles, exciseAll); + writeImages(sfiles, exciseImageCube); } while (streakFilesNextExtension(sfiles));