Index: trunk/magic/remove/src/streaksremove.c =================================================================== --- trunk/magic/remove/src/streaksremove.c (revision 20810) +++ trunk/magic/remove/src/streaksremove.c (revision 20816) @@ -1,424 +1,196 @@ #include "streaksremove.h" -#include "libgen.h" -#include "unistd.h" - -static nebServer *ourNebServer = NULL; - -// to enhance clarity we don't propagate errors up the stack -// we just bail out -void streaksremoveExit(psString str, int exitCode) { - psErrorStackPrint(stderr, str); - exit(exitCode); -} - -static void sFileFree(sFile *sfile) -{ - psFree(sfile->resolved_name); - psFree(sfile->name); - psFree(sfile); -} - - -// getNebServer() -// -// it gets created the first time a nebulous file is accessed. -// if config is passed in we are to create it and exit if not able. -// if config is null we return NULL indicating that there are no nebulous files in use -// The purpose of this is to avoid passing extra arguments between all of the functions. -// but it is probably a bad idea. -// -static nebServer *getNebServer(pmConfig *config) -{ - if (ourNebServer) { - return ourNebServer; - } - + +static pmConfig *parseArguments(int argc, char **argv); +static bool readAndCopyToOutput(streakFiles *sf, bool exciseAll); +static StreakPixels * getStreakPixels(streakFiles *sfiles, Streaks *streaks); +static void exciseNonWarpedPixels(streakFiles *sfiles, double newMaskValue); +static bool warpedPixel(streakFiles *sfiles, PixelPos *cellCoord); +static void excisePixel(streakFiles *sfiles, unsigned int x, unsigned int y, bool streak, double newMaskValue); +static void writeImages(streakFiles *sf, bool exciseImageCube); +static bool replicateOutputs(streakFiles *sfiles); + +int +main(int argc, char *argv[]) +{ + bool status; + + psLibInit(NULL); + + pmConfig *config = parseArguments(argc, argv); if (!config) { - return NULL; - } - - psString serverURI = getenv("NEB_SERVER"); - - // XXX: Note: all of the flags to psError should be true, but I don't think nebclient - // uses psError - if (!serverURI) { - bool status; - serverURI = psMetadataLookupStr(&status, config->site, "NEB_SERVER"); - if (!status) { - psError(PM_ERR_CONFIG, true, "failed to lookup config value for NEB_SERVER."); - streaksremoveExit("", PS_EXIT_CONFIG_ERROR); - } - } - - if (!serverURI) { - psError(PM_ERR_CONFIG, true, "Could not determine nebulous server URI."); - - streaksremoveExit("", PS_EXIT_CONFIG_ERROR); - } - - ourNebServer = nebServerAlloc(serverURI); - if (!ourNebServer) { - psError(PM_ERR_SYS, true, "failed to create a nebServer object from %s.", serverURI); - streaksremoveExit("", PS_EXIT_CONFIG_ERROR); - } - - return ourNebServer; -} - -static psString -resolveFilename(pmConfig *config, sFile *sfile, bool create) -{ - sfile->inNebulous = IN_NEBULOUS(sfile->name); - - if (sfile->inNebulous) { - // make sure we have our neb server connection - nebServer *server = getNebServer(config); - if (create) { - // delete the existing file, since there may be more than one - // instance. It will get created below in pmConfigConvertFilename - if (nebFind(server, sfile->name)) { - nebDelete(server, sfile->name); + psError(PS_ERR_UNKNOWN, false, "failed to parse arguments\n"); + return PS_EXIT_CONFIG_ERROR; + } + + psMetadata *masks = psMetadataLookupMetadata(&status, config->recipes, "MASKS"); + if (!status) { + psError(PM_ERR_CONFIG, false, "failed to lookup MASKS in recipes\n"); + return PS_EXIT_CONFIG_ERROR; + } + double maskStreak = (double) psMetadataLookupU8(&status, masks, "STREAK"); + if (!status) { + psError(PM_ERR_CONFIG, false, "failed to lookup mask value for STREAK in recipes\n"); + return PS_EXIT_CONFIG_ERROR; + } + + psString streaksFileName = psMetadataLookupStr(NULL, config->arguments, "STREAKS"); + + Streaks *streaks = readStreaksFile(streaksFileName); + if (!streaks) { + psErrorStackPrint(stderr, "failed to read streaks file: %s", streaksFileName); + exit (PS_EXIT_PROG_ERROR); + } + + streakFiles *sfiles = openFiles(config, true); + setupAstrometry(sfiles); + + 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 ) { + // From ICD: + // In the raw and detrended images, the pixels which were not + // included in any of the streak-processed warps must also be masked. + // 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. + // if no skycells are provided sfiles->exciseAll is set to true + + if (! computeWarpedPixels(sfiles) ) { + // we have no choice to excise all pixels + exciseAll = true; + } + } + + if (sfiles->stage == IPP_STAGE_RAW) { + // copy PHU to output files + copyPHU(sfiles, true); + + // advance to the first image extension + if (!streakFilesNextExtension(sfiles)) { + psErrorStackPrint(stderr, "failed to advance to first extension of input images"); + exit (PS_EXIT_PROG_ERROR); + } + } + + // Iterate through each component of the input (only one except for raw stage) + do { + 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; + } + + if (!exciseAll) { + // Identify pixels to mask because of streaks + + StreakPixels *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, maskStreak); } - } - } - - return pmConfigConvertFilename(sfile->name, config, create, create); -} - -sFile *sFileOpen(pmConfig *config, ippStage stage, psString fileSelect, - psString outputFilename, bool required) -{ - bool status; - sFile *sfile = psAlloc(sizeof(sFile)); - memset(sfile, 0, sizeof(sFile)); - - // We use psFits directly to read the image unless the stage is warp - // or diff. In those cases we use the pmFPAfile functions to read the image file - // to make managing the astrometry easy. - // The reason we don't use pmFPAfile in all cases is that I was having trouble getting - // all of the keywords in the raw image files written to the output destreaked files - - if (!CHIP_LEVEL_INPUT(stage) && !strcmp(fileSelect, "INPUT")) { - // stage is warp or diff AND fileSelect eq "INPUT" - // get data from pmFPAfile. - - // we need to know what the nebulous and real filenames are so we steal - // some code from pmFPAfileDefineFromArgs - // XXX: create a pmFPAfile function that does this and use it there - psArray *infiles = psMetadataLookupPtr(&status, config->arguments, "INPUT"); - if (!status) { - psError(PS_ERR_PROGRAMMING, false, "INPUT not found in config->arguments"); - streaksremoveExit("", PS_EXIT_PROG_ERROR); - } - if (infiles->n != 1) { - psError(PS_ERR_IO, false, "Found n == %ld files in %s in arguments expencted 1\n", - infiles->n, "INPUT"); - streaksremoveExit("", PS_EXIT_DATA_ERROR); - } - // end of file name lookup code adapted from pmFPAfileDefineFromArgs - // - sfile->name = psStringCopy(infiles->data[0]); - sfile->inNebulous = IN_NEBULOUS(sfile->name); - - // XXX: I should probably be using a different file rule for diffs, but I don't - // have an input rule that takes a diff image. - // Since they're skycells, this should be compatible - sfile->pmfile = pmFPAfileDefineFromArgs(&status, config, "PPSUB.INPUT", "INPUT"); - if (!sfile->pmfile) { - psError(PS_ERR_IO, false, "Failed to define file for %s", fileSelect); - streaksremoveExit("", PS_EXIT_DATA_ERROR); - } - - pmFPAview *view = pmFPAviewAlloc(0); - if (!pmFPAfileIOChecks(config, view, PM_FPA_BEFORE)) { - psError(PS_ERR_UNKNOWN, false, "Failed to load input."); - streaksremoveExit("", PS_EXIT_DATA_ERROR); - } - psFree(view); - sfile->resolved_name = psStringCopy(sfile->pmfile->filename); - - // copy header from fpu - sfile->header = (psMetadata*) psMemIncrRefCounter(sfile->pmfile->fpa->hdu->header); - - return sfile; - } - - // For all other files we use using psFits for i/o - - psString name = psMetadataLookupStr(&status, config->arguments, fileSelect); - if (!status || !name) { - if (required) { - psError(PS_ERR_IO, false, "Failed to lookup name for %s", fileSelect); - sFileFree(sfile); - streaksremoveExit("", 1); - } - return NULL; - } - - // if outputFilename is not null name it contains the "directory" - // and outputFilename is the basename name of the file (or nebulous key) - // and the file is to be opened for writing - if (outputFilename) { - psStringAppend(&sfile->name, "%s%s", name, outputFilename); - sfile->resolved_name = resolveFilename(config, sfile, true); - if (!sfile->resolved_name) { - psError(PS_ERR_IO, false, "Failed to resolve filename for %s", sfile->name); - sFileFree(sfile); - streaksremoveExit("", 1); - } - sfile->fits = psFitsOpen(sfile->resolved_name, "w"); - if (sfile->fits) { - sfile->fits->options = psFitsOptionsAlloc(); - } - } else { - sfile->name = psStringCopy(name); - sfile->resolved_name = resolveFilename(config, sfile, false); - if (!sfile->resolved_name) { - psError(PS_ERR_IO, false, "Failed to resolve name for %s", sfile->name); - sFileFree(sfile); - streaksremoveExit("", 1); - } - sfile->fits = psFitsOpen(sfile->resolved_name, "r"); - if (sfile->fits) { - sfile->nHDU = psFitsGetSize(sfile->fits); - } - } - - if (!sfile->fits) { - psError(PS_ERR_IO, false, "failed to open fits file %s for %s", - sfile->resolved_name, outputFilename ? "writing" : "reading"); - sFileFree(sfile); - streaksremoveExit("", 1); - } - - return sfile; -} - - -static bool -readAstrometry(streakFiles *sf) -{ - pmHDU *phu = pmFPAviewThisPHU(sf->view, sf->inAstrom->fpa); - if (phu) { - bool status; - char *ctype = psMetadataLookupStr(&status, phu->header, "CTYPE1"); - if (ctype) { - sf->bilevelAstrometry = !strcmp (&ctype[4], "-DIS"); - } - } - - if (sf->bilevelAstrometry) { - // Do we get here for GPC1 ? I don't necessarily have an fpa for the input image - if (!pmAstromReadBilevelMosaic(sf->inAstrom->fpa, phu->header)) { - psError(PS_ERR_UNKNOWN, false, "Unable to read bilevel mosaic astrometry for input FPA."); - return false; - } - } else { - pmHDU *hdu = pmFPAviewThisHDU(sf->view, sf->inAstrom->fpa); - PS_ASSERT_PTR_NON_NULL(hdu, 1) - PS_ASSERT_PTR_NON_NULL(hdu->header, 1) - - // we use a default FPA pixel scale of 1.0 - if (!pmAstromReadWCS (sf->inAstrom->fpa, sf->chip, hdu->header, 1.0)) { - psError(PS_ERR_UNKNOWN, false, "Unable to read WCS astrometry for input FPA."); - return false; - } - } - - return true; -} - -// load the fpa containing the astrometry, find the appropriate chip and process the data -static void -setupAstromFromFPA(streakFiles *sf) -{ - sf->view = pmFPAviewAlloc(0); - - if (!pmFPAfileIOChecks(sf->config, sf->view, PM_FPA_BEFORE)) { - psError(PS_ERR_UNKNOWN, false, "Failed to load input."); - streaksremoveExit("", PS_EXIT_DATA_ERROR); - } - - while ((sf->chip = pmFPAviewNextChip(sf->view, sf->inAstrom->fpa, 1))) { - if (sf->inAstrom->fpa->chips->n == 1) { - // There's only one chip in this FPA and we've got it so break - break; - } - bool status; - psString chip_name = psMetadataLookupStr(&status, sf->chip->concepts, "CHIP.NAME"); - if (!strcmp(chip_name, sf->class_id)) { - break; - } - } - if (!sf->chip) { - psError(PS_ERR_UNKNOWN, true, "Failed to find chip with data."); - streaksremoveExit("", PS_EXIT_DATA_ERROR); - } - if (!pmFPAfileIOChecks(sf->config, sf->view, PM_FPA_BEFORE)) { - psError(PS_ERR_UNKNOWN, false, "failed to load chip"); - streaksremoveExit("", PS_EXIT_DATA_ERROR); - } - - readAstrometry(sf); -} - -streakFiles *openFiles(pmConfig *config) -{ - bool status; - streakFiles *sf = psAlloc(sizeof(streakFiles)); - memset(sf, 0, sizeof(*sf)); - - sf->config = config; - - // error checking is done by sFileOpen. If a file can't be opened we just exit - ippStage stage = psMetadataLookupS32(&status, config->arguments, "STAGE"); - - sf->stage = stage; - sf->extnum = 0; - - sf->class_id = psMetadataLookupStr(&status, config->arguments, "CLASS_ID"); - - sf->inImage = sFileOpen(config, stage, "INPUT", NULL, true); - sf->nHDU = sf->inImage->nHDU; - - // don't need to free inputBasename see basename(3) - // The names of the temporary and recovery files are taken from the input - char *inputBasename = basename(sf->inImage->name); - sf->outImage = sFileOpen(config, stage, "OUTPUT", inputBasename, true); - sf->recImage = sFileOpen(config, stage, "RECOVERY", inputBasename, false); - - sf->inMask = sFileOpen(config, stage, "INPUT.MASK", NULL, false); - if (sf->inMask) { - inputBasename = basename(sf->inMask->name); - sf->outMask = sFileOpen(config, stage, "OUTPUT", inputBasename, true); - sf->recMask = sFileOpen(config, stage, "RECOVERY", inputBasename, true); - } - sf->inWeight = sFileOpen(config, stage, "INPUT.WEIGHT", NULL, false); - if (sf->inWeight) { - inputBasename = basename(sf->inWeight->name); - sf->outWeight = sFileOpen(config, stage, "OUTPUT", inputBasename, true); - sf->recWeight = sFileOpen(config, stage, "RECOVERY", inputBasename, true); - } - - // load astrometry file - if (USE_SUPPLIED_ASTROM(sf->stage)) { - sf->inAstrom = pmFPAfileDefineFromArgs(&status, config, "PSWARP.ASTROM", "ASTROM"); - } else { - // otherwise get astrometry from pmfile - if (!sf->inImage->pmfile) { - streaksremoveExit("unexpected null pmFPAfile", PS_EXIT_CONFIG_ERROR); - } - sf->inAstrom = sf->inImage->pmfile; - } - setupAstromFromFPA(sf); - - psElemType tileType; // Type corresponding to "long" - if (sizeof(long) == sizeof(psS64)) { - tileType = PS_TYPE_S64; - } else if (sizeof(long) == sizeof(psS32)) { - tileType = PS_TYPE_S32; - } else { - psAbort("can't map (long) type to a psLib type"); - } - - sf->tiles = psVectorAlloc(3, tileType); // Tile sizes - if (tileType == PS_TYPE_S64) { - sf->tiles->data.S64[0] = 0; - sf->tiles->data.S64[1] = 1; - sf->tiles->data.S64[2] = 1; - } else { - sf->tiles->data.S32[0] = 0; - sf->tiles->data.S32[1] = 1; - sf->tiles->data.S32[2] = 1; - } - - sf->transparentStreaks = psMetadataLookupF64(&status, config->arguments, "TRANSPARENT_STREAKS"); - - return sf; -} - -static bool -moveExt(sFile *sfile, int extnum) -{ - bool status = psFitsMoveExtNum(sfile->fits, extnum, false); - if (!status) { - psError(PS_ERR_IO, false, - "failed to move to extension %d for %s", extnum, sfile->resolved_name); - streaksremoveExit("", PS_EXIT_DATA_ERROR); - } - return true; -} - -static bool -streakFilesNextExtension(streakFiles *sf) -{ - // unless stage is raw or chip when we get here we're done - if (sf->stage != IPP_STAGE_RAW) { - sf->view->readout = -1; - pmFPAfileIOChecks(sf->config, sf->view, PM_FPA_AFTER); - sf->view->cell = -1; - pmFPAfileIOChecks(sf->config, sf->view, PM_FPA_AFTER); - sf->view->chip = -1; - pmFPAfileIOChecks(sf->config, sf->view, PM_FPA_AFTER); - return false; - } - - ++sf->extnum; - if (sf->nHDU == 1) { - // return true the first time through - return (sf->extnum == 1); - } - if (sf->extnum < sf->nHDU) { - moveExt(sf->inImage, sf->extnum); - if (sf->inMask) { - moveExt(sf->inMask, sf->extnum); - } - if (sf->inWeight) { - moveExt(sf->inWeight, sf->extnum); - } - return true; - } else { - return false; - } -} - -psString checkdir(char *path) -{ - // insure that path ends in / - int len = strlen(path); - psString dir = NULL; - if (strrchr(path, '/') != (path + len - 1)) { - psStringAppend(&dir, "%s/", path); - } else { - dir = psStringCopy(path); - } - return dir; -} - -ippStage -parseStage(psString stageStr) -{ - if (!strcmp("raw", stageStr)) { - return IPP_STAGE_RAW; - } else if (!strcmp("chip", stageStr)) { - return IPP_STAGE_CHIP; - } else if (!strcmp("warp", stageStr)) { - return IPP_STAGE_WARP; - } else if (!strcmp("diff", stageStr)) { - return IPP_STAGE_DIFF; - } else { - psError(PS_ERR_UNKNOWN, true, "unknown stage string: %s\n", stageStr); - streaksremoveExit("", PS_EXIT_DATA_ERROR); - // notreached - return IPP_STAGE_NONE; - } -} - - -pmConfig *parseArguments(int argc, char **argv) { + + if (sfiles->inImage->image) { + 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 + } + } + } else { + // this component contains an image cube, excise it completely + exciseImageCube = true; + } + psArrayElementsFree (pixels); + } + + // write out the destreaked temporary images and the recovery images + writeImages(sfiles, exciseImageCube); + + } while (streakFilesNextExtension(sfiles)); + + // close all files + closeImages(sfiles); + + // 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); + } + +#ifdef NOTYET + 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 + 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, ""); + 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); + } + } + } +#endif // REPLACE, REMOVE + // nebServerFree(ourNebServer); + psFree(config); + pmConceptsDone(); + pmConfigDone(); + psLibFinalize(); + + // PAU + fprintf(stderr, "Found %d leaks at %s\n", psMemCheckLeaks (0, NULL, NULL, false), "streaksremove"); + + return 0; +} + +static pmConfig *parseArguments(int argc, char **argv) +{ pmConfig *config = pmConfigRead(&argc, argv, NULL); if (!config) { - streaksremoveExit("failed to parse arguments", PS_EXIT_CONFIG_ERROR); + streaksExit("failed to parse arguments", PS_EXIT_CONFIG_ERROR); } @@ -550,5 +322,5 @@ return NULL; } - psString dir = checkdir(argv[argnum]); + psString dir = pathToDirectory(argv[argnum]); psMetadataAddStr(config->arguments, PS_LIST_TAIL, "OUTPUT", 0, "directory for temporary files", dir); @@ -561,5 +333,5 @@ if ((argnum = psArgumentGet(argc, argv, "-recovery"))) { psArgumentRemove(argnum, &argc, argv); - psString dir = checkdir(argv[argnum]); + psString dir = pathToDirectory(argv[argnum]); psMetadataAddStr(config->arguments, PS_LIST_TAIL, "RECOVERY", 0, "directory for recovery files", dir); @@ -591,320 +363,4 @@ return config; -} - -static void -copyPHU(streakFiles *sfiles) -{ - psAssert(sfiles->stage == IPP_STAGE_RAW, "copyPHU should only be used for raw stage"); - - psMetadata *imageHeader = psFitsReadHeader(NULL, sfiles->inImage->fits); - if (!imageHeader) { - psError(PS_ERR_IO, false, "failed to read primary header from %s\n", - sfiles->inImage->resolved_name); - streaksremoveExit("", PS_EXIT_DATA_ERROR); - } - - // TODO: add keyword indicating that streaks have been removed - if (!psFitsWriteBlank(sfiles->outImage->fits, imageHeader, NULL)) { - psError(PS_ERR_IO, false, "failed to write primary header to %s", - sfiles->outImage->resolved_name); - streaksremoveExit("", PS_EXIT_DATA_ERROR); - } - // TODO: add keyword indicating that this is the recovery image - 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); - streaksremoveExit("", PS_EXIT_DATA_ERROR); - } - psFree(imageHeader); - - sFile *inMask = sfiles->inMask; - if (inMask) { - psMetadata *maskHeader = psFitsReadHeader(NULL, inMask->fits); - if (!maskHeader) { - psError(PS_ERR_IO, false, "failed to read primary header from %s\n", - sfiles->inMask->resolved_name); - streaksremoveExit("", 1); - } - // TODO: add keyword indicating that streaks have been removed - if (!psFitsWriteBlank(sfiles->outMask->fits, maskHeader, NULL)) { - psError(PS_ERR_IO, false, "failed to write primary header to %s", - sfiles->outMask->resolved_name); - streaksremoveExit("", PS_EXIT_DATA_ERROR); - } - // TODO: add keyword indicating that this is the recovery image - 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); - streaksremoveExit("", PS_EXIT_DATA_ERROR); - } - psFree(maskHeader); - } - sFile *inWeight = sfiles->inWeight; - if (inWeight) { - psMetadata *weightHeader = psFitsReadHeader(NULL, inWeight->fits); - if (!weightHeader) { - psError(PS_ERR_IO, false, "failed to read primary header from %s\n", - sfiles->inWeight->resolved_name); - streaksremoveExit("", 1); - } - // TODO: add keyword indicating that streaks have been removed - if (!psFitsWriteBlank(sfiles->outWeight->fits, weightHeader, NULL)) { - psError(PS_ERR_IO, false, "failed to write primary header to %s", - sfiles->outWeight->resolved_name); - streaksremoveExit("", PS_EXIT_DATA_ERROR); - } - // TODO: add keyword indicating that this is a recovery image - 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); - streaksremoveExit("", PS_EXIT_DATA_ERROR); - } - psFree(weightHeader); - } -} - -// Determine whether the current header for this file is an image or not -// Find a cleaner way to do this -static bool -isImage(sFile *in) -{ - bool status; - psString exttype = psMetadataLookupStr(&status, in->header, "EXTTYPE"); - if (exttype && !strcmp(exttype, "IMAGE")) { - return true; - } - // examine current header and determine if it is an image - psString xtension = psMetadataLookupStr(&status, in->header, "XTENSION"); - if (xtension) { - if (!strcmp(xtension, "IMAGE")) { - return true; - } else if (!strcmp(xtension, "BINTABLE")) { - psString ztension = psMetadataLookupStr(&status, in->header, "ZTENSION"); - if (ztension && !strcmp(ztension, "IMAGE")) { - return true; - } - } - } else if (psMetadataLookupBool(&status, in->header, "ZIMAGE")) { - return true; - } else if (in->nHDU == 1) { - // no extensions in the file, can just return true? For now require - // at least one dimension - int naxis = psMetadataLookupS32(&status, in->header, "NAXIS"); - return naxis > 0; - } - return false; -} - -static void -readImageFrom_pmFile(streakFiles *sf) -{ - // XXX: This function assumes that it is only used for a single - // chip single cell FPA (i.e. a skycell) - pmFPAview *view = sf->view; - assert(view->chip == 0); - view->cell = 0; - sf->cell = pmFPAviewThisCell(view, sf->inImage->pmfile->fpa); - - if (!sf->cell) { - streaksremoveExit("no cells found in chip", PS_EXIT_PROG_ERROR); - } - if (sf->cell->readouts->n != 1) { - psError(PS_ERR_PROGRAMMING, true, "unexpected number of readouts found: %ld", sf->cell->readouts->n); - streaksremoveExit("", PS_EXIT_PROG_ERROR); - } - view->readout = 0; - pmReadout *readout = pmFPAviewThisReadout(view, sf->inImage->pmfile->fpa); - if (!readout) { - streaksremoveExit("readout 0 not found", PS_EXIT_PROG_ERROR); - } - - sf->inImage->image = (psImage*) psMemIncrRefCounter(readout->image); - sf->inImage->numCols = readout->image->numCols; - sf->inImage->numRows = readout->image->numRows; -} - -void -setDataExtent(ippStage stage, sFile *in) -{ - if (stage == IPP_STAGE_RAW) { - psString datasec = psMetadataLookupStr(NULL, in->header, "DATASEC"); - if (!datasec) { - psError(PS_ERR_IO, false, "failed to find DATASEC in header"); - streaksremoveExit("", PS_EXIT_DATA_ERROR); - } - int xmin, xmax, ymin, ymax; - sscanf(datasec, "[%d:%d,%d:%d]", &xmin, &xmax, &ymin, &ymax); - in->numCols = xmax - xmin + 1; - in->numRows = ymax - ymin + 1; - } else { - in->numCols = in->image->numCols; - in->numRows = in->image->numRows; - } -} - -static void -readImage(sFile *in, int extnum, ippStage stage) -{ - psRegion region = {0, 0, 0, 0}; - - if (in->header) psFree(in->header); - if (in->image) { - psFree(in->image); - in->image = NULL; - } - if (in->imagecube) { - psFree(in->imagecube); - in->imagecube = NULL; - } - - in->header = psFitsReadHeader(NULL, in->fits); - if (!in->header) { - psError(PS_ERR_IO, false, "failed to read header from %s extnum: %d", - in->resolved_name, extnum); - streaksremoveExit("", PS_EXIT_DATA_ERROR); - } - - if (!isImage(in)) { - return; - } - - in->numZPlanes = psMetadataLookupS32(NULL, in->header, "NAXIS3"); - - if (in->numZPlanes == 0) { - in->image = psFitsReadImage(in->fits, region, 0); - if (!in->image) { - psError(PS_ERR_IO, false, "failed to read image from %s extnum: %d", - in->resolved_name, extnum); - streaksremoveExit("", PS_EXIT_DATA_ERROR); - } - } else { - if (stage != IPP_STAGE_RAW) { - psError(PS_ERR_PROGRAMMING, true, "unexpected image cube found in non-raw image"); - streaksremoveExit("", PS_EXIT_PROG_ERROR); - } - in->imagecube = psFitsReadImageCube(in->fits, region); - if (!in->imagecube) { - psError(PS_ERR_IO, false, "failed to read image cube from %s extnum: %d", - in->resolved_name, extnum); - streaksremoveExit("", PS_EXIT_DATA_ERROR); - } - psImage *image = (psImage *) (in->imagecube->data[0]); - } - setDataExtent(stage, in); -} - -static void -setFitsOptions(sFile *sfile, int bitpix, float bscale, float bzero) -{ - if (!sfile) { - return; - } - if (sfile->fits->options) { - psFree(sfile->fits->options); - } - sfile->fits->options = psFitsOptionsAlloc(); - sfile->fits->options->scaling = PS_FITS_SCALE_MANUAL; - sfile->fits->options->bitpix = bitpix; - sfile->fits->options->bscale = bscale; - sfile->fits->options->bzero = bzero; -} - -static void -copyFitsOptions(sFile *out, sFile *rec, sFile *in) -{ - // Get current BITPIX, BSCALE, BZERO, EXTNAME - // Probably not necessary to look the numerical values up in this - // way, but guards against changes to psLib and cfitsio FITS - // handling. - psMetadataItem *bitpixItem = psMetadataLookup(in->header, "BITPIX"); - psAssert(bitpixItem, "Every FITS image should have BITPIX"); - int bitpix = psMetadataItemParseS32(bitpixItem); - psMetadataItem *bscaleItem = psMetadataLookup(in->header, "BSCALE"); - - float bscale; - if (!bscaleItem) { - psWarning("BSCALE isn't set; defaulting to unity"); - bscale = 1.0; - } else { - bscale = psMetadataItemParseF32(bscaleItem); - } - psMetadataItem *bzeroItem = psMetadataLookup(in->header, "BZERO"); - float bzero; - if (!bzeroItem) { - psWarning("BZERO isn't set; defaulting to zero"); - bzero = 0.0; - } else { - bzero = psMetadataItemParseF32(bzeroItem); - } - -#ifdef COMPRESS_OUTPUT - setFitsOptions(out, bitpix, bscale, bzero); - setFitsOptions(rec, bitpix, bscale, bzero); -#endif -} - - -static void -copyTable(sFile *out, sFile *in, int extnum) -{ - bool mdok; - psString xtension = psMetadataLookupStr(&mdok, in->header, "XTENSION"); - psString extname = psMetadataLookupStr(&mdok, in->header, "EXTNAME"); - - if (!xtension) { - psWarning("extnum %d has no image and xtension is not defined in %s", - extnum, in->resolved_name); - // XXX abort? - return; - } - if (!extname) { - psWarning("extnum %d has no image and extname not defined in %s", - extnum, in->resolved_name); - // XXX abort? - return; - } - psArray *table = psFitsReadTable(in->fits); - if (!table) { - psError(PS_ERR_UNKNOWN, false, "failed to read table in extension %d from in->resolved name", extnum); - streaksremoveExit("", PS_EXIT_DATA_ERROR); - } - - if (!psFitsWriteTable(out->fits, out->header, table, extname)) { - psError(PS_ERR_UNKNOWN, false, "failed to copy table in extension %d", extnum); - streaksremoveExit("", PS_EXIT_DATA_ERROR); - } -} - -static void -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(out->image, initValue); -} - -static void -setupImageRefs(sFile *out, sFile *recovery, sFile *in, int extnum, bool exciseAll) -{ - if (!exciseAll) { - // set output image to be the input image - out->image = (psImage*) psMemIncrRefCounter(in->image); - if (recovery) { - // create a recovery image initialized to NAN - createNewImage(recovery, in, extnum, NAN); - } - } else { - // Excising all pixels in the input - // 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); - } - } } @@ -929,5 +385,5 @@ if (!sf->astrom) { psError(PS_ERR_UNKNOWN, false, "failed to set up astrometry for extension %d", sf->extnum); - streaksremoveExit("", PS_EXIT_UNKNOWN_ERROR); + streaksExit("", PS_EXIT_UNKNOWN_ERROR); } if (sf->stage == IPP_STAGE_CHIP) { @@ -936,5 +392,5 @@ 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); - streaksremoveExit("", PS_EXIT_UNKNOWN_ERROR); + streaksExit("", PS_EXIT_UNKNOWN_ERROR); } } @@ -958,5 +414,5 @@ // set up the compression parameters -#ifdef COMPRESS_OUTPUT +#ifdef STREAKS_COMPRESS_OUTPUT // compression of the image pixels is disabled for now. Some consortium members // have problems reading them @@ -967,8 +423,8 @@ // 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); -#endif if (sf->recImage) { psFitsSetCompression(sf->recImage->fits, PS_FITS_COMPRESS_RICE, sf->tiles, 8, 0, 0); } +#endif if (sf->inMask) { @@ -1010,39 +466,8 @@ // area (no overscan) - return true; } -static void -writeImage(sFile *sfile, psString extname, int extnum) -{ - if (!sfile || !sfile->image) { - return; - } - if (!psFitsWriteImage(sfile->fits, sfile->header, sfile->image, 0, extname)) { - psError(PS_ERR_IO, false, "failed to write image to %s extnum: %d", - sfile->resolved_name, extnum); - streaksremoveExit("", PS_EXIT_DATA_ERROR); - } - psFree(sfile->image); - sfile->image = NULL; - psFree(sfile->header); - sfile->header = NULL; -} - -static void -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", - sfile->resolved_name, extnum); - streaksremoveExit("", PS_EXIT_DATA_ERROR); - } - psFree(sfile->header); - sfile->header = NULL; -} + static void @@ -1098,61 +523,4 @@ } -static void -closeImage(sFile *sfile) -{ - if (!sfile) { - return; - } - if (sfile->fits && !psFitsClose(sfile->fits)) { - psError(PS_ERR_IO, false, "failed to close image to %s", - sfile->resolved_name); - streaksremoveExit("", PS_EXIT_DATA_ERROR); - } - psFree(sfile->header); - psFree(sfile->image); - psFree(sfile->imagecube); - psFree(sfile->name); - psFree(sfile->resolved_name); - psFree(sfile); -} - -static void -closeImages(streakFiles *sf) -{ - closeImage(sf->inImage); - closeImage(sf->outImage); - closeImage(sf->recImage); - if (sf->inMask) { - closeImage(sf->inMask); - closeImage(sf->outMask); - closeImage(sf->recMask); - } - if (sf->inWeight) { - closeImage(sf->inWeight); - closeImage(sf->outWeight); - closeImage(sf->recWeight); - } -} - -static bool -replicate(sFile *sfile, void *xattr) -{ - if (!sfile->inNebulous) { - return true; - } - nebServer *server = getNebServer(NULL); - - // for now just set "user.copies" to 2 - if (!nebSetXattr(server, sfile->name, "user.copies", "2", NEB_REPLACE)) { - psError(PM_ERR_UNKNOWN, true, "nebSetXattr failed for %s\n%s", sfile->name, nebErr(server)); - return false; - } - if (!nebReplicate(server, sfile->name, NULL, NULL)) { - psError(PM_ERR_UNKNOWN, true, "nebSetXattr failed for %s\n%s", sfile->name, nebErr(server)); - return false; - } - return true; -} - // for any of the outputs that are stored in Nebulous set their extended attributes // and replicate the files @@ -1200,106 +568,4 @@ } -static bool -swapOutputToInput(sFile *in, sFile *out) -{ - - if (in->inNebulous) { - nebServer *server = getNebServer(NULL); - - if (!nebSwap(server, in->name, out->name)) { - psError(PM_ERR_SYS, true, "failed to swap files for: %s.", - in->name, nebErr(server)); - return false; - } - } else { - psString command = NULL; - - // regular files. use mv creating a backup - // NOTE: -b is a gnu option and may not always be available - psStringAppend(&command, "mv -b --suffix=.bak %s %s", out->resolved_name, in->resolved_name); - int status = system(command); - if (status != 0) { - psError(PM_ERR_SYS, true, "%s failed with status %d.", - command, status); - psFree(command); - return false; - } - psFree(command); - // XXX: shouldn't I be doing mv in->resolved_name.bak out->resolved_name - // to complete the swap if !remove? - } - - return true; -} -static bool -swapOutputsToInputs(streakFiles *sfiles) -{ - if (sfiles->inMask) { - if (!swapOutputToInput(sfiles->inMask, sfiles->outMask)) { - psError(PM_ERR_SYS, false, "failed to swap instances for Mask."); - return false; - } - } - - if (sfiles->inWeight) { - if (!swapOutputToInput(sfiles->inWeight, sfiles->outWeight)) { - psError(PM_ERR_SYS, false, "failed to swap instances for Weight."); - return false; - } - } - - if (!swapOutputToInput(sfiles->inImage, sfiles->outImage)) { - psError(PM_ERR_SYS, false, "failed to swap instances for Image."); - return false; - } - - return true; -} - -static bool -deleteFile(sFile *sfile) -{ - - if (sfile->inNebulous) { - nebServer *server = getNebServer(NULL); - if (!nebDelete(server, sfile->name)) { - psError(PM_ERR_SYS, false, "failed to delete %s\n%s.", sfile->name, - nebErr(server)); - return false; - } - } else { - if (unlink(sfile->resolved_name)) { - psError(PM_ERR_SYS, false, "failed to delete %s\n%s.", sfile->resolved_name, - strerror(errno)); - return false; - } - } - return true; -} - -static bool -deleteTemps(streakFiles *sfiles) -{ - if (sfiles->outMask) { - if (!deleteFile(sfiles->outMask)) { - psError(PM_ERR_SYS, false, "failed to delete Mask."); - return false; - } - } - - if (sfiles->outWeight) { - if (!deleteFile(sfiles->outWeight)) { - psError(PM_ERR_SYS, false, "failed to delete Weight."); - return false; - } - } - - if (!deleteFile(sfiles->outImage)) { - psError(PM_ERR_SYS, false, "failed to delete Image."); - return false; - } - - return true; -} static void @@ -1312,5 +578,5 @@ double imageValue = psImageGet (sfiles->inImage->image, x, y); - if (sfiles->recImage) { + if (sfiles->recImage && !isnan(imageValue) ) { psImageSet (sfiles->recImage->image, x, y, imageValue); } @@ -1327,8 +593,4 @@ } - // TODO: - // 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 (sfiles->inWeight) { double weightValue = psImageGet (sfiles->inWeight->image, x, y); @@ -1385,5 +647,5 @@ } -bool +static bool warpedPixel(streakFiles *sfiles, PixelPos *cellCoord) { @@ -1419,5 +681,5 @@ } -bool +static bool streakEndOnComponent(streak *streak, double r_min, double r_max, double d_min, double d_max) { @@ -1434,5 +696,5 @@ -Streaks * +static Streaks * restrictStreaksToComponent(streakFiles *sfiles, Streaks *streaks) { @@ -1494,181 +756,3 @@ return pixels; } -int -main(int argc, char *argv[]) -{ - long i; - bool status; - StreakPixels *pixels; - PixelPos *pixelPos; - - psLibInit(NULL); - - pmConfig *config = parseArguments(argc, argv); - if (!config) { - psError(PS_ERR_UNKNOWN, false, "failed to parse arguments\n"); - return PS_EXIT_CONFIG_ERROR; - } - - psMetadata *masks = psMetadataLookupMetadata(&status, config->recipes, "MASKS"); - if (!status) { - psError(PM_ERR_CONFIG, false, "failed to lookup MASKS in recipes\n"); - return PS_EXIT_CONFIG_ERROR; - } - double maskStreak = (double) psMetadataLookupU8(&status, masks, "STREAK"); - if (!status) { - psError(PM_ERR_CONFIG, false, "failed to lookup mask value for STREAK in recipes\n"); - return PS_EXIT_CONFIG_ERROR; - } - - psString streaksFileName = psMetadataLookupStr(NULL, config->arguments, "STREAKS"); - - Streaks *streaks = readStreaksFile(streaksFileName); - if (!streaks) { - psErrorStackPrint(stderr, "failed to read streaks file: %s", streaksFileName); - exit (PS_EXIT_PROG_ERROR); - } - - streakFiles *sfiles = openFiles(config); - - 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 ) { - // From ICD: - // In the raw and detrended images, the pixels which were not - // included in any of the streak-processed warps must also be masked. - // 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. - // if no skycells are provided sfiles->exciseAll is set to true - - if (! computeWarpedPixels(sfiles) ) { - // we have no choice to excise all pixels - exciseAll = true; - } - } - - - if (sfiles->stage == IPP_STAGE_RAW) { - // copy PHU to output files - copyPHU(sfiles); - - // advance to the first image extension - if (!streakFilesNextExtension(sfiles)) { - psErrorStackPrint(stderr, "failed to advance to first extension of input images"); - exit (PS_EXIT_PROG_ERROR); - } - } - - // Iterate through components of input files - do { - 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; - } - - 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, maskStreak); - } - - 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, true, maskStreak); - - } else { - // This pixel was not included in any warp and has thus already excised - // by exciseNonWarpedPixels - } - } - } else { - // this component contains an image cube, excise it completely - exciseImageCube = true; - } - psArrayElementsFree (pixels); - } - - // write out the destreaked temporary images and the recovery images - writeImages(sfiles, exciseImageCube); - - } while (streakFilesNextExtension(sfiles)); - - // close all files - closeImages(sfiles); - - // 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); - } - -#ifdef NOTYET - 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 - 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, ""); - 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); - } - } - } -#endif // REPLACE, REMOVE - nebServer(ourNebServer); - psFree(config); - pmConceptsDone(); - pmConfigDone(); - psLibFinalize(); - - // PAU - fprintf(stderr, "Found %d leaks at %s\n", psMemCheckLeaks (0, NULL, NULL, false), "streaksremove"); - - return 0; -} +