Index: /trunk/magic/remove/src/.cvsignore =================================================================== --- /trunk/magic/remove/src/.cvsignore (revision 20815) +++ /trunk/magic/remove/src/.cvsignore (revision 20816) @@ -1,1 +1,3 @@ streaksremove +streaksreplace +streakscompare Index: /trunk/magic/remove/src/Makefile.simple =================================================================== --- /trunk/magic/remove/src/Makefile.simple (revision 20815) +++ /trunk/magic/remove/src/Makefile.simple (revision 20816) @@ -1,5 +1,13 @@ # skeleton Makefile for streaksremove +# +# The intent is that these programs may be built against an +# IPP installation outside the IPP build enviornment. -OBJECTS= \ +COMMON_OBJECTS = \ + streaksio.o \ + streaksutil.o + +REMOVE_OBJECTS= \ + ${COMMON_OBJECTS} \ streaksremove.o \ streaksastrom.o \ @@ -8,15 +16,33 @@ Line.o -CFLAGS=`psmodules-config --cflags` -g -std=gnu99 +REPLACE_OBJECTS= \ + ${COMMON_OBJECTS} \ + streaksreplace.o + +COMPARE_OBJECTS= \ + ${COMMON_OBJECTS} \ + streakscompare.o + +OPTFLAGS= -g +CFLAGS=`psmodules-config --cflags` -std=gnu99 ${OPTFLAGS} LDFLAGS=`psmodules-config --libs` -streaksremove: ${OBJECTS} +PROGRAMS= streaksremove streaksreplace streakscompare -install: streaksremove - cp streaksremove $(PSCONFDIR)/$(PSCONFIG)/bin - chmod 755 $(PSCONFDIR)/$(PSCONFIG)/bin/streaksremove +all: ${PROGRAMS} + +streaksremove: ${REMOVE_OBJECTS} + +streaksreplace: ${REPLACE_OBJECTS} + +streakscompare: ${COMPARE_OBJECTS} + +install: ${PROGRAMS} + install -t $(PSCONFDIR)/$(PSCONFIG)/bin streaksremove + install -t $(PSCONFDIR)/$(PSCONFIG)/bin streaksreplace + install -t $(PSCONFDIR)/$(PSCONFIG)/bin streakscompare clean: - rm -f *.o streaksremove + rm -f *.o ${PROGRAMS} Index: /trunk/magic/remove/src/streaksastrom.c =================================================================== --- /trunk/magic/remove/src/streaksastrom.c (revision 20815) +++ /trunk/magic/remove/src/streaksastrom.c (revision 20816) @@ -150,2 +150,80 @@ } +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; +} + +void +setupAstrometry(streakFiles *sf) +{ + bool status; + // load astrometry file + if (USE_SUPPLIED_ASTROM(sf->stage)) { + sf->inAstrom = pmFPAfileDefineFromArgs(&status, sf->config, "PSWARP.ASTROM", "ASTROM"); + } else { + // otherwise get astrometry from pmfile + if (!sf->inImage->pmfile) { + streaksExit("unexpected null pmFPAfile", PS_EXIT_CONFIG_ERROR); + } + sf->inAstrom = sf->inImage->pmfile; + } + sf->view = pmFPAviewAlloc(0); + + if (!pmFPAfileIOChecks(sf->config, sf->view, PM_FPA_BEFORE)) { + psError(PS_ERR_UNKNOWN, false, "Failed to load input."); + streaksExit("", 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."); + streaksExit("", PS_EXIT_DATA_ERROR); + } + if (!pmFPAfileIOChecks(sf->config, sf->view, PM_FPA_BEFORE)) { + psError(PS_ERR_UNKNOWN, false, "failed to load chip"); + streaksExit("", PS_EXIT_DATA_ERROR); + } + if (!readAstrometry(sf)) { + psError(PS_ERR_UNKNOWN, false, "failed to read astrometry"); + streaksExit("", PS_EXIT_DATA_ERROR); + } +} Index: /trunk/magic/remove/src/streakscompare.c =================================================================== --- /trunk/magic/remove/src/streakscompare.c (revision 20816) +++ /trunk/magic/remove/src/streakscompare.c (revision 20816) @@ -0,0 +1,153 @@ +#include "streaksremove.h" + +static pmConfig * parseArguments(int, char **); + +int main(int argc, char *argv[]) +{ + long i; + bool status; + + psLibInit(NULL); + + pmConfig *config = parseArguments(argc, argv); + if (!config) { + streaksExit("failed to parse arguments", PS_EXIT_CONFIG_ERROR); + } + + ippStage stage = psMetadataLookupS32(&status, config->arguments, "STAGE"); + + sFile *file1 = sFileOpen(config, stage, "INPUT1", NULL, true); + sFile *file2 = sFileOpen(config, stage, "INPUT2", NULL, true); + + readImage(file1, 0, stage); + readImage(file2, 0, stage); + + psImage *image1 = file1->image; + psImage *image2 = file2->image; + + if ((image1->numRows != image2->numRows) || (image1->numCols != image2->numCols)) { + streaksExit("image sizes do not match\n", PS_EXIT_DATA_ERROR); + } + + int width = image1->numCols; + int height = image1->numRows; + int numErrors = 0; + for (int y = 0; y < height; y++) { + for (int x = 0; x < width; x++) { + bool error = false; + + double value1 = psImageGet(image1, x, y); + double value2 = psImageGet(image1, x, y); + + if (!isnan(value1) && !isnan(value2)) { + if (value1 != value2) { + error = true; + } + } else { + // if one's a nan they had both better be + error = ! isnan(value1) && isnan(value2); + } + + if (error) { + numErrors++; + fprintf(stderr, "pixel at %5d %5d do not match %f %f\n", x, y, value1, value2); + } + } + } + + fprintf(stderr, "%d errors\n", numErrors); + + closeImage(file1); + closeImage(file2); + + psFree(config); + psLibFinalize(); + + fprintf(stderr, "Found %d leaks at %s\n", psMemCheckLeaks (0, NULL, NULL, false), "streakscompare"); + + // PAU + return 0; +} +static void usage(void) +{ + fprintf(stderr, "USAGE: streakscompare file1 file2\n"); + + exit(2); +} + +static pmConfig *parseArguments(int argc, char **argv) +{ + int argnum; + + pmConfig *config = pmConfigRead(&argc, argv, NULL); + if (!config) { + return NULL; + } + + ippStage stage = IPP_STAGE_NONE; + if ((argnum = psArgumentGet(argc, argv, "-stage"))) { + psArgumentRemove(argnum, &argc, argv); + stage = parseStage(argv[argnum]); + psMetadataAddS32(config->arguments, PS_LIST_TAIL, "STAGE", 0, + "pipeline stage for streak removal", stage); + psArgumentRemove(argnum, &argc, argv); + } else { + psError(PS_ERR_UNKNOWN, true, "-stage is required\n"); + usage(); + } + + // If we are asked to replace inputs with the outputs and input image is in nebulous + // we require mask, weight, and tmproot to all be + bool nebulousImage = false; + if ((argnum = psArgumentGet(argc, argv, "-image1"))) { + // for raw and chip level images we use psFits for I/O and ... + nebulousImage = IN_NEBULOUS(argv[argnum+1]); + if (CHIP_LEVEL_INPUT(stage)) { + psArgumentRemove(argnum, &argc, argv); + psMetadataAddStr(config->arguments, PS_LIST_TAIL, "INPUT1", 0, + "name of input image", argv[argnum]); + psArgumentRemove(argnum, &argc, argv); + } else { + // ... for warped images we use pmFPAfiles + if (!pmConfigFileSetsMD(config->arguments, &argc, argv, "INPUT1", "-image1", NULL)) { ; + psError(PS_ERR_UNKNOWN, false, "failed to process -image"); + usage(); + } + } + } else { + psError(PS_ERR_UNKNOWN, true, "-image1 is required\n"); + usage(); + } + if ((argnum = psArgumentGet(argc, argv, "-image2"))) { + // for raw and chip level images we use psFits for I/O and ... + nebulousImage = IN_NEBULOUS(argv[argnum+1]); + if (CHIP_LEVEL_INPUT(stage)) { + psArgumentRemove(argnum, &argc, argv); + psMetadataAddStr(config->arguments, PS_LIST_TAIL, "INPUT2", 0, + "name of input image", argv[argnum]); + psArgumentRemove(argnum, &argc, argv); + } else { + // ... for warped images we use pmFPAfiles + if (!pmConfigFileSetsMD(config->arguments, &argc, argv, "INPUT2", "-image2", NULL)) { ; + psError(PS_ERR_UNKNOWN, false, "failed to process -image"); + usage(); + } + } + } else { + psError(PS_ERR_UNKNOWN, true, "-image2 is required\n"); + usage(); + } + + if (argc != 1) { + psString unexpectedArguments = NULL; + for (int i=1; iconfig = 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); + + if (remove) { + sf->recImage = sFileOpen(config, stage, "RECOVERY", inputBasename, false); + } else { + sf->recImage = sFileOpen(config, stage, "RECOVERY.IMAGE", NULL, true); + } + + 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); + if (remove) { + sf->recMask = sFileOpen(config, stage, "RECOVERY", inputBasename, false); + } else { + sf->recMask = sFileOpen(config, stage, "RECOVERY.MASK", NULL, 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); + if (remove) { + sf->recWeight = sFileOpen(config, stage, "RECOVERY", inputBasename, false); + } else { + sf->recWeight = sFileOpen(config, stage, "RECOVERY.WEIGHT", NULL, true); + } + } + + 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; +} + +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; + } + + 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."); + streaksExit("", PS_EXIT_CONFIG_ERROR); + } + } + + if (!serverURI) { + psError(PM_ERR_CONFIG, true, "Could not determine nebulous server URI."); + + streaksExit("", PS_EXIT_CONFIG_ERROR); + } + + ourNebServer = nebServerAlloc(serverURI); + if (!ourNebServer) { + psError(PM_ERR_SYS, true, "failed to create a nebServer object from %s.", serverURI); + streaksExit("", 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); + } + } + } + + 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"); + streaksExit("", 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"); + streaksExit("", 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); + streaksExit("", PS_EXIT_DATA_ERROR); + } + + pmFPAview *view = pmFPAviewAlloc(0); + if (!pmFPAfileIOChecks(config, view, PM_FPA_BEFORE)) { + psError(PS_ERR_UNKNOWN, false, "Failed to load input."); + streaksExit("", 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); + streaksExit("", 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); + streaksExit("", 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); + streaksExit("", 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); + streaksExit("", 1); + } + + return sfile; +} + + +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); + streaksExit("", PS_EXIT_DATA_ERROR); + } + return true; +} + +bool +streakFilesNextExtension(streakFiles *sf) +{ + // unless stage is raw or chip when we get here we're done + if (sf->stage != IPP_STAGE_RAW) { + if (sf->view) { + 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; + } +} + +void +copyPHU(streakFiles *sfiles, bool remove) +{ + 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); + streaksExit("", 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); + streaksExit("", PS_EXIT_DATA_ERROR); + } + // TODO: add keyword indicating that this is the recovery image + if (remove && sfiles->recImage && !psFitsWriteBlank(sfiles->recImage->fits, imageHeader, NULL)) { + psError(PS_ERR_IO, false, "failed to write primary header to %s", + sfiles->recImage->resolved_name); + streaksExit("", 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); + streaksExit("", 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); + streaksExit("", PS_EXIT_DATA_ERROR); + } + // TODO: add keyword indicating that this is the recovery image + if (remove && sfiles->recMask && !psFitsWriteBlank(sfiles->recMask->fits, maskHeader, NULL)) { + psError(PS_ERR_IO, false, "failed to write primary header to %s", + sfiles->recMask->resolved_name); + streaksExit("", 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); + streaksExit("", 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); + streaksExit("", PS_EXIT_DATA_ERROR); + } + // TODO: add keyword indicating that this is a recovery image + if (remove && sfiles->recWeight && !psFitsWriteBlank(sfiles->recWeight->fits, weightHeader, NULL)) { + psError(PS_ERR_IO, false, "failed to write primary header to %s", + sfiles->recWeight->resolved_name); + streaksExit("", 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; +} + +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) { + streaksExit("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); + streaksExit("", PS_EXIT_PROG_ERROR); + } + view->readout = 0; + pmReadout *readout = pmFPAviewThisReadout(view, sf->inImage->pmfile->fpa); + if (!readout) { + streaksExit("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"); + streaksExit("", 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; + } +} + +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); + streaksExit("", 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); + streaksExit("", PS_EXIT_DATA_ERROR); + } + } else { + if (stage != IPP_STAGE_RAW) { + psError(PS_ERR_PROGRAMMING, true, "unexpected image cube found in non-raw image"); + streaksExit("", 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); + streaksExit("", 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; +} + +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 STREAKS_COMPRESS_OUTPUT + setFitsOptions(out, bitpix, bscale, bzero); + setFitsOptions(rec, bitpix, bscale, bzero); +#endif +} + + +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); + streaksExit("", 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); + streaksExit("", 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); + streaksExit("", PS_EXIT_UNKNOWN_ERROR); + } + psImageInit(out->image, initValue); +} + +void +setupImageRefs(sFile *out, sFile *recoveryOut, sFile *in, int extnum, bool exciseAll) +{ + if (!exciseAll) { + // set output image to be the input image + out->image = (psImage*) psMemIncrRefCounter(in->image); + if (recoveryOut) { + // create a recovery image initialized to NAN + createNewImage(recoveryOut , in, extnum, NAN); + } + } else { + // Excising all pixels in the input + // set the output to an image all NANs + createNewImage(out, in, extnum, NAN); + if (recoveryOut) { + // save the whole input as the recovery image + recoveryOut->image = (psImage*) psMemIncrRefCounter(in->image); + } + } +} + +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); + streaksExit("", PS_EXIT_DATA_ERROR); + } + psFree(sfile->image); + sfile->image = NULL; + psFree(sfile->header); + sfile->header = NULL; +} + +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); + streaksExit("", PS_EXIT_DATA_ERROR); + } + psFree(sfile->header); + sfile->header = NULL; +} + +#ifdef notdef +void +writeImages(streakFiles *sf, bool exciseImageCube) +{ + psString extname = NULL; + if (sf->nHDU > 1) { + bool mdok; + extname = psMetadataLookupStr(&mdok, sf->inImage->header, "EXTNAME"); + } + 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 (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); + if (exciseImageCube) { + sf->outImage->image = image; + writeImage(sf->outImage, extname, sf->extnum); + } else { + // 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; + } + if (sf->outMask) { + writeImage(sf->outMask, extname, sf->extnum); + writeImage(sf->recMask, extname, sf->extnum); + } + if (sf->outWeight) { + writeImage(sf->outWeight, extname, sf->extnum); + writeImage(sf->recWeight, extname, sf->extnum); + } +} +#endif + +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); + streaksExit("", PS_EXIT_DATA_ERROR); + } + psFree(sfile->header); + psFree(sfile->image); + psFree(sfile->imagecube); + psFree(sfile->name); + psFree(sfile->resolved_name); + psFree(sfile); +} + +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); + } +} + +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 + +static bool +replicateOutputs(streakFiles *sfiles) +{ + bool status = false; + + // XXX: TODO: need a nebGetXatrr function, but there isn't one + // another option would be to take the number of copies to be + // created as an option. That way the system could decide + // whether to replicate anything other than raw Image files + void *xattr = NULL; + + if (!replicate(sfiles->outImage, xattr)) { + psError(PM_ERR_SYS, false, "failed to replicate outImage."); + return false; + } + +#ifdef notyet + // XXX: don't replicate mask and weight images until we can look up + // the input's xattr. There may be a perl program that can getXattr + if (sfiles->outMask) { + // get xattr from input to see if we need to replicate + if (!replicate(sfiles->outMask, xattr)) { + psError(PM_ERR_SYS, false, "failed to replicate outImage."); + return false; + } + } + if (sfiles->outWeight) { + // get xattr from input to see if we need to replicate + if (!replicate(sfiles->outWeight, xattr)) { + psError(PM_ERR_SYS, false, "failed to replicate outImage."); + return false; + } + } +#endif + +// replicate the recovery images (if in nebulous) +// perhaps whether we do that or not should be configurable. +// Sounds like we need a recipe + + return true; +} + +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; +} + Index: /trunk/magic/remove/src/streaksio.h =================================================================== --- /trunk/magic/remove/src/streaksio.h (revision 20816) +++ /trunk/magic/remove/src/streaksio.h (revision 20816) @@ -0,0 +1,27 @@ +#ifndef STREAKS_IO_H +#define STREAKS_IO_H 1 + +streakFiles *openFiles(pmConfig *config, bool remove); + +sFile *sFileOpen(pmConfig *config, ippStage stage, psString fileSelect, + psString outputFilename, bool required); + +void sFileFree(sFile *); + +void closeImages(streakFiles *sfiles); +void closeImage(sFile *sfile); + +void readImage(sFile *sfile, int extnum, ippStage stage); +void copyPHU(streakFiles *sfiles, bool remove); +void copyTable(sFile *out, sFile *in, int extnum); +void copyFitsOptions(sFile *out, sFile *rec, sFile *in); +void setupImageRefs(sFile *out, sFile *recoveryOut, sFile *in, int extnum, bool exciseAll); + +void writeImage(sFile *sfile, psString extname, int extnum); +void writeImageCube(sFile *sfile, psArray *imagecube, psString extname, int extnum); +bool replicate(sFile *sfile, void *xattr); +void readImageFrom_pmFile(streakFiles *sf); + +bool streakFilesNextExtension(streakFiles *sf); + +#endif // STREAKS_IO_H Index: /trunk/magic/remove/src/streaksremove.c =================================================================== --- /trunk/magic/remove/src/streaksremove.c (revision 20815) +++ /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; -} + Index: /trunk/magic/remove/src/streaksremove.h =================================================================== --- /trunk/magic/remove/src/streaksremove.h (revision 20815) +++ /trunk/magic/remove/src/streaksremove.h (revision 20816) @@ -10,8 +10,4 @@ #include "psmodules.h" #include "nebclient.h" - -typedef struct { - int dummy; -} Warps; #include "streaksextern.h" @@ -70,4 +66,7 @@ } streakFiles; +#include "streaksio.h" + +extern void setupAstrometry(streakFiles *); extern strkAstrom *streakSetAstrometry(strkAstrom *, ippStage stage, pmFPA *, pmChip *, bool fromCell, psMetadata *md, int numCols, int numRows); @@ -75,7 +74,8 @@ extern void cellToChip(PixelPos *chip, strkAstrom *astrom, PixelPos *cell); - extern bool computeWarpedPixels(streakFiles *sf); -extern void streaksremoveExit(psString, int); +extern void streaksExit(psString, int); +extern ippStage parseStage(psString); +extern psString pathToDirectory(char *path); #define CHIP_LEVEL_INPUT(_stage) ((_stage == IPP_STAGE_RAW) || (_stage == IPP_STAGE_CHIP)) @@ -85,4 +85,3 @@ #define IN_NEBULOUS(_filename) (!strncasecmp(_filename, "neb://", strlen("neb://"))) - #endif // STREAKS_H Index: /trunk/magic/remove/src/streaksreplace.c =================================================================== --- /trunk/magic/remove/src/streaksreplace.c (revision 20816) +++ /trunk/magic/remove/src/streaksreplace.c (revision 20816) @@ -0,0 +1,472 @@ +#include "streaksremove.h" + +static pmConfig *parseArguments(int, char **); +static bool readAndCopyToOutput(streakFiles *sf, bool restoreImageCube); +static void writeImages(streakFiles *sf, bool exciseImageCube); +static bool replicateOutputs(streakFiles *sfiles); +static bool readAndCopyToOutput(streakFiles *sf, bool restoreImageCube); + +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; + } + // TODO: perhaps we should save this in the image header and get it from there? + 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; + } + + streakFiles *sfiles = openFiles(config, false); + + if (sfiles->stage == IPP_STAGE_RAW) { + // copy PHU to output files + copyPHU(sfiles, false); + + // 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 restoreImageCube = 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, restoreImageCube)) { + // 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 pixel in recovery image is not NAN, then it contains an excised pixel + // that has been stored in the recovery image + for (int y=0; y < sfiles->inImage->numRows; y++) { + for (int x=0; x < sfiles->inImage->numCols; x++) { + double recValue = psImageGet(sfiles->recImage->image, x, y); + if (!isnan(recValue)) { + psImageSet(sfiles->inImage->image, x, y, recValue); + } + if (sfiles->inMask) { + double recMaskValue = psImageGet(sfiles->recMask->image, x, y); + if (recValue != maskStreak) { + psImageSet(sfiles->inMask->image, x, y, recMaskValue); + } + } + if (sfiles->inWeight) { + double recWeightValue = psImageGet(sfiles->recWeight->image, x, y); + if (recWeightValue != NAN) { + psImageSet(sfiles->inWeight->image, x, y, recWeightValue); + } + } + } + } + + // write out the re-streaked images + writeImages(sfiles, restoreImageCube); + + } 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), "streaksreplace"); + + return 0; +} +static void usage(void) +{ + fprintf(stderr, "USAGE: streaksreplace -stage raw|chip|warp|diff -image IMAGE.fits -tmproot directory\n\n"); + fprintf(stderr, "Optional arguments:\n"); + fprintf(stderr, "\t-mask MASK.fits mask file to re-streak\n"); + fprintf(stderr, "\t-recmask RECOVERYMASK.fits: recovery mask file (required with -mask)\n"); + fprintf(stderr, "\t-weight WEIGHT.fits: weight file to re-streak\n"); + fprintf(stderr, "\t-recweight RECOVERYWEIGHT.fits: recovery weight file (required with -weight)\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"); + + exit(2); +} + +static pmConfig *parseArguments(int argc, char **argv) +{ + pmConfig *config = pmConfigRead(&argc, argv, NULL); + if (!config) { + streaksExit("failed to parse arguments", PS_EXIT_CONFIG_ERROR); + } + + int argnum; + ippStage stage = IPP_STAGE_NONE; + if ((argnum = psArgumentGet(argc, argv, "-stage"))) { + psArgumentRemove(argnum, &argc, argv); + stage = parseStage(argv[argnum]); + psMetadataAddS32(config->arguments, PS_LIST_TAIL, "STAGE", 0, + "pipeline stage for streak removal", stage); + psArgumentRemove(argnum, &argc, argv); + } else { + psError(PS_ERR_UNKNOWN, true, "-stage is required\n"); + usage(); + } + + // If we are asked to replace inputs with the outputs and input image is in nebulous + // we require mask, weight, and tmproot to all be + bool nebulousImage = false; + if ((argnum = psArgumentGet(argc, argv, "-image"))) { + // for raw and chip level images we use psFits for I/O and ... + nebulousImage = IN_NEBULOUS(argv[argnum+1]); + if (CHIP_LEVEL_INPUT(stage)) { + psArgumentRemove(argnum, &argc, argv); + psMetadataAddStr(config->arguments, PS_LIST_TAIL, "INPUT", 0, + "name of input image", argv[argnum]); + psArgumentRemove(argnum, &argc, argv); + } else { + // ... for warped images we use pmFPAfiles + if (!pmConfigFileSetsMD(config->arguments, &argc, argv, "INPUT", "-image", NULL)) { ; + psError(PS_ERR_UNKNOWN, false, "failed to process -image"); + usage(); + } + } + } else { + psError(PS_ERR_UNKNOWN, true, "-image is required\n"); + usage(); + } + + bool gotReplace = false; + if ((argnum = psArgumentGet(argc, argv, "-replace"))) { + gotReplace = true; + psArgumentRemove(argnum, &argc, argv); + psMetadataAddBool(config->arguments, PS_LIST_TAIL, "REPLACE", 0, "replace input files", + true); + } + + bool gotMask = false; + if ((argnum = psArgumentGet(argc, argv, "-mask"))) { + psArgumentRemove(argnum, &argc, argv); + bool nebulousMask = IN_NEBULOUS(argv[argnum]); + if (nebulousMask != nebulousImage) { + psError(PS_ERR_UNKNOWN, true, "mask image must be %snebulous path with %s image path\n", + nebulousImage ? "" : "non ", nebulousImage ? "nebulous" : "regular"); + usage(); + } + gotMask = true; + psMetadataAddStr(config->arguments, PS_LIST_TAIL, "INPUT.MASK", 0, + "name of input mask image", argv[argnum]); + psArgumentRemove(argnum, &argc, argv); + } + + bool gotWeight = false; + if ((argnum = psArgumentGet(argc, argv, "-weight"))) { + psArgumentRemove(argnum, &argc, argv); + bool nebulousWeight = IN_NEBULOUS(argv[argnum]); + if (nebulousWeight != nebulousImage) { + psError(PS_ERR_UNKNOWN, true, "weight image must have %snebulous path with %s image path\n", + nebulousImage ? "" : "non ", nebulousImage ? "nebulous" : "regular"); + usage(); + } + gotWeight = true; + psMetadataAddStr(config->arguments, PS_LIST_TAIL, "INPUT.WEIGHT", 0, + "name of input weight image", argv[argnum]); + psArgumentRemove(argnum, &argc, argv); + } + + if ((argnum = psArgumentGet(argc, argv, "-tmproot"))) { + psArgumentRemove(argnum, &argc, argv); + bool nebulousTmp = IN_NEBULOUS(argv[argnum]); + if (gotReplace && (nebulousTmp != nebulousImage)) { + psError(PS_ERR_UNKNOWN, true, "tmproot must have %snebulous path with %s image path with -replace\n", + nebulousImage ? "" : "non ", nebulousImage ? "nebulous" : "regular"); + usage(); + } + psString dir = pathToDirectory(argv[argnum]); + psMetadataAddStr(config->arguments, PS_LIST_TAIL, "OUTPUT", 0, "directory for temporary files", + dir); + psArgumentRemove(argnum, &argc, argv); + } else { + psError(PS_ERR_UNKNOWN, true, "-tmproot is required\n"); + usage(); + } + + if ((argnum = psArgumentGet(argc, argv, "-recimage"))) { + psArgumentRemove(argnum, &argc, argv); + psMetadataAddStr(config->arguments, PS_LIST_TAIL, "RECOVERY.IMAGE", 0, "recovery image file name", + argv[argnum]); + psArgumentRemove(argnum, &argc, argv); + } else { + psError(PS_ERR_UNKNOWN, true, "-recimage is required\n"); + usage(); + } + if ((argnum = psArgumentGet(argc, argv, "-recmask"))) { + psArgumentRemove(argnum, &argc, argv); + bool nebulousMask = IN_NEBULOUS(argv[argnum]); + if (nebulousMask != nebulousImage) { + psError(PS_ERR_UNKNOWN, true, "mask image must be %snebulous path with %s image path\n", + nebulousImage ? "" : "non ", nebulousImage ? "nebulous" : "regular"); + usage(); + } + psMetadataAddStr(config->arguments, PS_LIST_TAIL, "INPUT.MASK", 0, + "name of input mask image", argv[argnum]); + psArgumentRemove(argnum, &argc, argv); + } else if (gotMask) { + psError(PS_ERR_UNKNOWN, true, "-recmask is required with -mask\n"); + usage(); + } + + if ((argnum = psArgumentGet(argc, argv, "-recweight"))) { + psArgumentRemove(argnum, &argc, argv); + bool nebulousWeight = IN_NEBULOUS(argv[argnum]); + if (nebulousWeight != nebulousImage) { + psError(PS_ERR_UNKNOWN, true, "weight image must have %snebulous path with %s image path\n", + nebulousImage ? "" : "non ", nebulousImage ? "nebulous" : "regular"); + usage(); + } + psMetadataAddStr(config->arguments, PS_LIST_TAIL, "INPUT.WEIGHT", 0, + "name of input weight image", argv[argnum]); + psArgumentRemove(argnum, &argc, argv); + } else if (gotWeight) { + psError(PS_ERR_UNKNOWN, true, "-recweight is required with -weight\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); + } + + if (argc != 1) { + psString unexpectedArguments = NULL; + for (int i=1; iinImage->pmfile) { + // image data from pmFPAfile (diff or warp file) + readImageFrom_pmFile(sf); + } else { + // image data directly from psFits + readImage(sf->inImage, sf->extnum, sf->stage); + } + // read the recovery pixels + readImage(sf->recImage, sf->extnum, sf->stage); + + sf->outImage->header = (psMetadata*) psMemIncrRefCounter(sf->inImage->header); + + if (!SFILE_IS_IMAGE(sf->inImage)) { + // fprintf(stderr, "skipping extnum: %d not an image\n", sf->extnum); + // XXX: This extension is probably a video table. + copyTable(sf->outImage, sf->inImage, sf->extnum); + return false; + } + + if (sf->inImage->image) { + setupImageRefs(sf->outImage, NULL, sf->inImage, sf->extnum, false); + } else { + // Image cubes are handeled specially + } + + // set up the compression parameters +#ifdef STREAKS_COMPRESS_OUTPUT + // compression of the image pixels is disabled for now. Some consortium members + // have problems reading compressed images. + 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); +#endif + if (sf->inMask) { + readImage(sf->inMask, sf->extnum, sf->stage); + readImage(sf->recMask, sf->extnum, sf->stage); + sf->outMask->header = (psMetadata*) psMemIncrRefCounter(sf->inMask->header); + + // XXX: TODO + setupImageRefs(sf->outMask, NULL, sf->inMask, sf->extnum, false); + + // XXX: see note above + copyFitsOptions(sf->outMask, NULL, sf->inMask); + psFitsSetCompression(sf->outMask->fits, PS_FITS_COMPRESS_PLIO, sf->tiles, 8, 0, 0); + } + + if (sf->inWeight) { + readImage(sf->inWeight, sf->extnum, sf->stage); + readImage(sf->recWeight, sf->extnum, sf->stage); + + sf->outWeight->header = (psMetadata*) psMemIncrRefCounter(sf->inWeight->header); + setupImageRefs(sf->outMask, NULL, sf->inMask, sf->extnum, false); + + copyFitsOptions(sf->outWeight, NULL, sf->inWeight); + // XXX: see note above + psFitsSetCompression(sf->outWeight->fits, PS_FITS_COMPRESS_RICE, sf->tiles, 8, 0, 0); + } + + return true; +} + +static void writeImages(streakFiles *sf, bool exciseImageCube) +{ + psString extname = NULL; + if (sf->nHDU > 1) { + bool mdok; + extname = psMetadataLookupStr(&mdok, sf->inImage->header, "EXTNAME"); + } + if (sf->inImage->numZPlanes == 0) { + // note exciseing complete images is handled in readAndCopyToOutput + writeImage(sf->outImage, extname, sf->extnum); + } else { +#ifdef TODO + // we have an image cube + double initValue; + // otherwise write it to the output + writeImageCube(sf->outImage, sf->recImage->imagecube, extname, sf->extnum); + + // 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 (exciseImageCube) { + sf->outImage->image = image; + writeImage(sf->outImage, extname, sf->extnum); + } else { + // write zero valued image to reccovery + if (sf->recImage) { + sf->recImage->image = image; + writeImage(sf->recImage, extname, sf->extnum); + } + } +#endif + // Assumption: there are no mask and weight images for video cells because they only appear + // at the raw stage. TODO: insure this + return; + } + if (sf->outMask) { + writeImage(sf->outMask, extname, sf->extnum); + } + if (sf->outWeight) { + writeImage(sf->outWeight, extname, sf->extnum); + } +} + +static bool replicateOutputs(streakFiles *sfiles) +{ + bool status = false; + + // XXX: TODO: need a nebGetXatrr function, but there isn't one + // another option would be to take the number of copies to be + // created as an option. That way the system could decide + // whether to replicate anything other than raw Image files + void *xattr = NULL; + + if (!replicate(sfiles->outImage, xattr)) { + psError(PM_ERR_SYS, false, "failed to replicate outImage."); + return false; + } + +#ifdef notyet + // XXX: don't replicate mask and weight images until we can look up + // the input's xattr. There may be a perl program that can getXattr + if (sfiles->outMask) { + // get xattr from input to see if we need to replicate + if (!replicate(sfiles->outMask, xattr)) { + psError(PM_ERR_SYS, false, "failed to replicate outImage."); + return false; + } + } + if (sfiles->outWeight) { + // get xattr from input to see if we need to replicate + if (!replicate(sfiles->outWeight, xattr)) { + psError(PM_ERR_SYS, false, "failed to replicate outImage."); + return false; + } + } +#endif + +// replicate the recovery images (if in nebulous) +// perhaps whether we do that or not should be configurable. +// Sounds like we need a recipe + + return true; +} + Index: /trunk/magic/remove/src/streaksutil.c =================================================================== --- /trunk/magic/remove/src/streaksutil.c (revision 20816) +++ /trunk/magic/remove/src/streaksutil.c (revision 20816) @@ -0,0 +1,41 @@ +#include "streaksremove.h" + +psString pathToDirectory(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); + streaksExit("", PS_EXIT_DATA_ERROR); + // notreached + return IPP_STAGE_NONE; + } +} + +// to enhance clarity in these programs we don't propagate errors up the stack +// we just bail out +void streaksExit(psString str, int exitCode) { + psErrorStackPrint(stderr, str); + exit(exitCode); +} + Index: /trunk/magic/remove/src/warpedpixels.c =================================================================== --- /trunk/magic/remove/src/warpedpixels.c (revision 20815) +++ /trunk/magic/remove/src/warpedpixels.c (revision 20816) @@ -57,10 +57,10 @@ if (!fits) { psError(PS_ERR_IO, false, "failed to open skycell file: %s", fileName); - streaksremoveExit("", PS_EXIT_DATA_ERROR); + streaksExit("", PS_EXIT_DATA_ERROR); } psMetadata *header = psFitsReadHeader(NULL, fits); if (!header) { psError(PS_ERR_IO, false, "failed to read fixts header from skycell file: %s", fileName); - streaksremoveExit("", PS_EXIT_DATA_ERROR); + streaksExit("", PS_EXIT_DATA_ERROR); } // set up astrometry for conversion from skycell image to sky @@ -68,5 +68,5 @@ if (!wcs) { psError(PS_ERR_IO, false, "failed to read astrometry from skycell file: %s", fileName); - streaksremoveExit("", PS_EXIT_DATA_ERROR); + streaksExit("", PS_EXIT_DATA_ERROR); } int naxis1 = psMetadataLookupS32(NULL, header, "NAXIS1"); @@ -99,5 +99,5 @@ if (!pmAstromWCStoSky(&sky, wcs, &pt[i])) { psError(PS_ERR_IO, false, "failed to convert pt %d of %s to sky coords: %s", fileName); - streaksremoveExit("", PS_EXIT_DATA_ERROR); + streaksExit("", PS_EXIT_DATA_ERROR); } strkPt p; @@ -105,5 +105,5 @@ if (!skyToCell(&p, sf->astrom, sky.r, sky.d)) { psError(PS_ERR_IO, false, "failed to convert pt %d of %s to sky coords: %s", fileName); - streaksremoveExit("", PS_EXIT_DATA_ERROR); + streaksExit("", PS_EXIT_DATA_ERROR); } pt[i].x = p.x;