Changeset 20816 for trunk/magic/remove/src

Timestamp:

Nov 21, 2008, 1:33:36 PM (17 years ago)

Author:

bills

Message:

Added a program to replace streaks and to compare two image files
(comparsion not complete yet)
Split streaksremove.c up so to make it less wieldy and to facilitate
reusing the io functions.

Location:

trunk/magic/remove/src

Files:

• . (modified) (1 prop)
• .cvsignore (modified) (1 diff)
• Makefile.simple (modified) (2 diffs)
• streaksastrom.c (modified) (1 diff)
• streakscompare.c (added)
• streaksio.c (added)
• streaksio.h (added)
• streaksremove.c (modified) (17 diffs)
• streaksremove.h (modified) (4 diffs)
• streaksreplace.c (added)
• streaksutil.c (added)
• warpedpixels.c (modified) (4 diffs)

trunk/magic/remove/src

@@ -1 +1 @@
 streaksremove
+streaksreplace
+streakscompare

@@ -1 +1 @@
 streaksremove
+streaksreplace
+streakscompare

trunk/magic/remove/src/.cvsignore

@@ -1 +1 @@
 streaksremove
+streaksreplace
+streakscompare

trunk/magic/remove/src/Makefile.simple

@@ -1 +1 @@
 # 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 +16 @@
     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}
     
 

trunk/magic/remove/src/streaksastrom.c

@@ -150 +150 @@
 }
  
+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);
+    }
+}

trunk/magic/remove/src/streaksremove.c

@@ -1 +1 @@
 #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 +322 @@
             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 +333 @@
     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 +363 @@
 
     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 +385 @@
     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 +392 @@
         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 +414 @@
 
     // 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 +423 @@
     // 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 +466 @@
     // 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 +523 @@
 }
 
-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 +568 @@
 }
 
-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 +578 @@
 
     double imageValue  = psImageGet (sfiles->inImage->image,  x, y);
-    if (sfiles->recImage) {
+    if (sfiles->recImage && !isnan(imageValue) ) {
         psImageSet (sfiles->recImage->image,  x, y, imageValue);
     }
@@ -1327 +593 @@
     }
 
-    // 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 +647 @@
 }
 
-bool
+static bool
 warpedPixel(streakFiles *sfiles, PixelPos *cellCoord)
 {
@@ -1419 +681 @@
 }
 
-bool
+static bool
 streakEndOnComponent(streak *streak, double r_min, double r_max, double d_min, double d_max)
 {
@@ -1434 +696 @@
 
 
-Streaks *
+static Streaks *
 restrictStreaksToComponent(streakFiles *sfiles, Streaks *streaks)
 {
@@ -1494 +756 @@
     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;
-}
+

trunk/magic/remove/src/streaksremove.h

@@ -10 +10 @@
 #include "psmodules.h"
 #include "nebclient.h"
-
-typedef struct {
-    int dummy;
-} Warps;
 
 #include "streaksextern.h"
@@ -70 +66 @@
 } 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 +74 @@
 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 +85 @@
 #define IN_NEBULOUS(_filename) (!strncasecmp(_filename, "neb://", strlen("neb://")))
 
-
 #endif // STREAKS_H

trunk/magic/remove/src/warpedpixels.c

@@ -57 +57 @@
     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 +68 @@
     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 +99 @@
         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 +105 @@
         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;

@@ -1 +1 @@
 streaksremove
+streaksreplace
+streakscompare