Index: branches/tap_branches/magic/remove/src/Makefile.am
===================================================================
--- branches/tap_branches/magic/remove/src/Makefile.am	(revision 27838)
+++ branches/tap_branches/magic/remove/src/Makefile.am	(revision 27838)
@@ -0,0 +1,95 @@
+if HAVE_SVNVERSION
+STREAKS_VERSION=`$(SVNVERSION) ..`
+else
+STREAKS_VERSION="UNKNOWN"
+endif
+
+if HAVE_SVN
+STREAKS_BRANCH=`$(SVN) info .. | $(SED) -n -e '/URL:/ h' -e '/Repository Root:/ { x; H; x; s|Repository Root: \(.*\)\nURL: \1\(.*\)|\2| ; s|^/|| ; s|/[a-zA-Z]*/src.*|| ; p }'`
+STREAKS_SOURCE=`$(SVN) info | $(SED) -n -e 's/Repository UUID: // p'`
+else
+STREAKS_BRANCH="UNKNOWN"
+STREAKS_SOURCE="UNKNOWN"
+endif
+
+# Force recompilation of streaksVersion.c, since it gets the version information
+streaksVersion.c: streaksVersionDefinitions.h
+streaksVersionDefinitions.h: streaksVersionDefinitions.h.in FORCE
+	-$(RM) streaksVersionDefinitions.h
+	$(SED) -e "s|@STREAKS_VERSION@|\"$(STREAKS_VERSION)\"|" -e "s|@STREAKS_BRANCH@|\"$(STREAKS_BRANCH)\"|" -e "s|@STREAKS_SOURCE@|\"$(STREAKS_SOURCE)\"|" streaksVersionDefinitions.h.in > streaksVersionDefinitions.h
+FORCE: ;
+
+bin_PROGRAMS = streaksremove streaksrelease streaksreplace streakscompare isdestreaked
+
+STREAKSREMOVE_CFLAGS=-DSTREAKS_COMPRESS_OUTPUT=1
+
+STREAKSCFLAGS = $(STREAKSREMOVE_CFLAGS) $(PSMODULE_CFLAGS) $(PSLIB_CFLAGS)
+STREAKSLDFLAGS = $(STREAKSREMOVE_LIBS) $(PSMODULE_LIBS) $(PSLIB_LIBS)
+
+streaksremove_CFLAGS = $(STREAKSCFLAGS)
+streaksremove_LDFLAGS = $(STREAKSLDFLAGS)
+streaksreplace_CFLAGS = $(STREAKSCFLAGS)
+streaksreplace_LDFLAGS = $(STREAKSLDFLAGS)
+streaksrelease_CFLAGS = $(STREAKSCFLAGS)
+streaksrelease_LDFLAGS = $(STREAKSLDFLAGS)
+streakscompare_CFLAGS = $(STREAKSCFLAGS)
+streakscompare_LDFLAGS = $(STREAKSLDFLAGS)
+isdestreaked_CFLAGS = $(STREAKSCFLAGS)
+isdestreaked_LDFLAGS = $(STREAKSLDFLAGS)
+
+COMMON_SOURCES = 	\
+	streaksio.c 	\
+	streaksutil.c 	\
+	streaksastrom.c \
+	streaksVersion.c
+
+streaksremove_SOURCES = \
+	streaksremove.c	\
+	streaksextern.c	\
+	diffedpixels.c	\
+	Line.c		\
+	$(COMMON_SOURCES)
+
+streaksrelease_SOURCES = \
+	streaksrelease.c\
+	streaksio.c	\
+	streaksutil.c	\
+	streaksastrom.c \
+	streaksVersion.c
+
+streaksreplace_SOURCES = \
+	streaksreplace.c\
+	streaksio.c	\
+	streaksutil.c	\
+	streaksastrom.c \
+	streaksVersion.c
+
+streakscompare_SOURCES = \
+	streakscompare.c\
+	streaksio.c	\
+	streaksutil.c	\
+	streaksastrom.c \
+	streaksVersion.c
+
+isdestreaked_SOURCES = \
+	isdestreaked.c\
+	streaksio.c	\
+	streaksutil.c	\
+	streaksastrom.c \
+	streaksVersion.c
+
+include_HEADERS = 	\
+	streaksremove.h \
+	streaksastrom.h \
+	streaksextern.h \
+	streaksio.h 	\
+	streaksVersionDefinitions.h
+
+clean-local:
+	-rm -f TAGS
+	-rm -f tags
+
+# Tags for emacs and vi
+tags:
+	etags `find . -name \*.[ch] -print`
+	ctags `find . -name \*.[ch] -print`
Index: branches/tap_branches/magic/remove/src/Makefile.simple
===================================================================
--- branches/tap_branches/magic/remove/src/Makefile.simple	(revision 25900)
+++ branches/tap_branches/magic/remove/src/Makefile.simple	(revision 27838)
@@ -10,24 +10,19 @@
 
 REMOVE_OBJECTS=    \
-    ${COMMON_OBJECTS} \
     streaksremove.o \
     streaksextern.o \
-    warpedpixels.o \
+    diffedpixels.o \
     Line.o
 
 REPLACE_OBJECTS=      \
-    ${COMMON_OBJECTS} \
     streaksreplace.o
 
 COMPARE_OBJECTS=      \
-    ${COMMON_OBJECTS} \
     streakscompare.o
 
 RELEASE_OBJECTS=      \
-    ${COMMON_OBJECTS} \
     streaksrelease.o
 
 ISDESTREAKED_OBJECTS=      \
-    ${COMMON_OBJECTS} \
     isdestreaked.o
 
@@ -45,16 +40,27 @@
 HEADERS=Line.h streaksastrom.h streaksextern.h streaksio.h streaksremove.h
 
-all:	${PROGRAMS}
+all: ${PROGRAMS}
 
-${REMOVE_OBJECTS}:	${HEADERS}
-streaksremove:  ${REMOVE_OBJECTS}
+# programs all depend on the common objects and the common headers
+PROGRAM_OBJECTS = \
+ ${REMOVE_OBJECTS} \
+ ${REPLACE_OBJECTS} \
+ ${COMPARE_OBJECTS} \
+ ${RELEASE_OBJECTS} \
+ ${ISDESTREAKED_OBJECTS}
 
-streaksreplace:  ${REPLACE_OBJECTS}
+${COMMON_OBJECTS}: ${HEADERS}
 
-streakscompare:  ${COMPARE_OBJECTS}
+${PROGRAM_OBJECTS}: ${HEADERS} ${COMMON_OBJECTS}
 
-streaksrelease:  ${RELEASE_OBJECTS}
+streaksremove:  ${REMOVE_OBJECTS} ${COMMON_OBJECTS}
 
-isdestreaked:	${ISDESTREAKED_OBJECTS}
+streaksreplace:  ${REPLACE_OBJECTS} ${COMMON_OBJECTS}
+
+streakscompare:  ${COMPARE_OBJECTS} ${COMMON_OBJECTS}
+
+streaksrelease:  ${RELEASE_OBJECTS} ${COMMON_OBJECTS}
+
+isdestreaked:	${ISDESTREAKED_OBJECTS} ${COMMON_OBJECTS}
 
 
Index: branches/tap_branches/magic/remove/src/diffedpixels.c
===================================================================
--- branches/tap_branches/magic/remove/src/diffedpixels.c	(revision 27838)
+++ branches/tap_branches/magic/remove/src/diffedpixels.c	(revision 27838)
@@ -0,0 +1,277 @@
+#include "streaksremove.h"
+#include "pmAstrometryWCS.h"
+
+#define DEBUG_PRINT 1
+
+static void addTouchedPixels(streakFiles *sf, psString fileName);
+static void nameCorners(psPlane pt[4]);
+static int xLeft(psPlane *pt, int y);
+static int xRight(psPlane *pt, int y);
+
+#define  BL 0
+#define  BR 1
+#define  TL 2
+#define  TR 3
+
+// Examine the set of diff skycells and compute the pixels that were
+// included in difference processing
+bool
+computeDiffedPixels(streakFiles *sf)
+{
+    bool status;
+    pmConfig    *config = sf->config;
+
+    psArray *skycells = psMetadataLookupPtr(&status, config->arguments, "SKYCELLS");
+    if (skycells == NULL) {
+        return false;
+    }
+
+    int n = psArrayLength(skycells);
+
+    psRegion     *bounds = pmChipPixels(sf->chip);
+
+    int width  = bounds->x1 - bounds->x0;
+    int height = bounds->y1 - bounds->y0;
+
+    psFree(bounds);
+
+    sf->diffedPixels = psImageAlloc(width, height, PS_TYPE_U8);
+    psImageInit(sf->diffedPixels, 0);
+
+    for (int i=0; i<n; i++) {
+        psString filename = psArrayGet(skycells, i);
+        psString resolved_name = pmConfigConvertFilename(filename, config, false, false);
+        if (resolved_name == NULL) {
+            psError(PS_ERR_IO, false, "failed to resolve skycell file: %s", filename);
+            streaksExit("", PS_EXIT_DATA_ERROR);
+        }
+        addTouchedPixels(sf, resolved_name);
+        psFree(resolved_name);
+    }
+
+    bool writeTouchedPixels = false; // XXX: make this an argument to the program
+    if (writeTouchedPixels) {
+        // write out the diffed pixels
+        psMetadata * header = psMetadataAlloc();
+        if (sf->inAstrom->fpa->toSky->type != PS_PROJ_DIS) {
+            pmAstromWriteWCS(header, sf->inAstrom->fpa, sf->chip, 0.001);
+        }
+        psFits *fits = psFitsOpen("diffedpixels.fits", "w");
+
+        psFitsWriteImage(fits, header, sf->diffedPixels, 0, NULL);
+        psFitsClose(fits);
+        psFree(header);
+    }
+    return true;
+}
+
+static void
+addTouchedPixels(streakFiles *sf, psString fileName)
+{
+    // read fits header
+    psFits *fits = psFitsOpen(fileName, "r");
+    if (!fits) {
+        psError(PS_ERR_IO, false, "failed to open skycell file: %s", fileName);
+        streaksExit("", PS_EXIT_DATA_ERROR);
+    }
+    psMetadata *header = psFitsReadHeader(NULL, fits);
+    psFitsClose(fits);
+    if (!header) {
+        psError(PS_ERR_IO, false, "failed to read fixts header from skycell file: %s", fileName);
+        streaksExit("", PS_EXIT_DATA_ERROR);
+    }
+    // set up astrometry for conversion from skycell image to sky
+    pmAstromWCS *wcs = pmAstromWCSfromHeader(header);
+    if (!wcs) {
+        psError(PS_ERR_IO, false, "failed to read astrometry  from skycell file: %s", fileName);
+        streaksExit("", PS_EXIT_DATA_ERROR);
+    }
+    int naxis1 = psMetadataLookupS32(NULL, header, "NAXIS1");
+    int naxis2 = psMetadataLookupS32(NULL, header, "NAXIS2");
+
+    psFree(header);
+
+    /* now set up our wrapper to the chip astrometry to apply to the whole chip */
+    sf->astrom = streakSetAstrometry(sf->astrom, sf->stage, sf->inAstrom->fpa, sf->chip, false, NULL,
+        sf->diffedPixels->numCols, sf->diffedPixels->numRows);
+
+    // convert corners of skycell to chip coordinates
+    // compute overlap with the chip using astrometry
+    psPlane pt[4];
+    pt[0].x = 0;
+    pt[0].y = 0;
+
+    pt[1].x = naxis1 - 1;
+    pt[1].y = 0;
+
+    pt[2].x = naxis1 - 1;
+    pt[2].y = naxis2 - 1;
+
+    pt[3].x = 0;
+    pt[3].y = naxis2 - 1;
+
+    for (int i = 0; i < 4 ; i++ ) {
+        psSphere sky;
+
+        // convert corner of skycell to sky coordinates
+        if (!pmAstromWCStoSky(&sky, wcs, &pt[i])) {
+            psErrorStackPrint(stderr, "failed to convert pt %d of %s to sky coords", i, fileName);
+            streaksExit("", PS_EXIT_DATA_ERROR);
+        }
+        strkPt p;
+        // convert to chip coordinates
+        if (!skyToCell(&p, sf->astrom, sky.r, sky.d)) {
+            psErrorStackPrint(stderr, "failed to convert pt %d of %s to sky coords", i, fileName);
+            streaksExit("", PS_EXIT_DATA_ERROR);
+        }
+        pt[i].x = round(p.x);
+        pt[i].y = round(p.y);
+    }
+
+    psFree(wcs);
+
+    // Identify the bottom left, bottom right, top left, and top right corners
+    nameCorners(pt);
+
+#ifdef DEBUG_PRINT
+    fprintf(stderr, "\nSKYCELL %s\n", fileName);
+    fprintf(stderr, "BL %f %f\n", pt[BL].x, pt[BL].y);
+    fprintf(stderr, "BR %f %f\n", pt[BR].x, pt[BR].y);
+    fprintf(stderr, "TL %f %f\n", pt[TL].x, pt[TL].y);
+    fprintf(stderr, "TR %f %f\n", pt[TR].x, pt[TR].y);
+#endif
+
+    int xmax = sf->diffedPixels->numCols - 1;
+    int ymin;
+    int ymax; 
+    
+    if (pt[BL].y <= pt[BR].y) {
+        ymin = fmax(0, pt[BL].y);
+    } else {
+        ymin = fmax(0, pt[BR].y);
+    }
+    if (pt[TR].y >= pt[TL].y) {
+        ymax = fmin(pt[TR].y, sf->diffedPixels->numRows - 1);
+    } else {
+        ymax = fmin(pt[TL].y, sf->diffedPixels->numRows - 1);
+    }
+        
+#if (DEBUG_PRINT > 1)
+    fprintf(stderr, "\nxmin: %d xmax: %d\n", 0, xmax);
+    fprintf(stderr, "\nymin: %d ymax: %d\n", ymin, ymax);
+#endif
+
+    // Now set the touched pixels
+    for (int y = ymin ; y <= ymax; y++) {
+
+        int xleft  = xLeft(pt, y) - 1;
+        int xright = xRight(pt, y) + 1;
+
+        if (xright < 0) {
+            continue;
+        }
+        if (xleft > xmax) {
+            continue;
+        }
+        if (xleft < 0) {
+            xleft = 0;
+        }
+        if (xright > xmax) {
+            xright = xmax;
+        }
+#if (DEBUG_PRINT > 1)
+        fprintf(stderr, "  y: %d xleft: %d xright: %d\n", y, xleft, xright);
+#endif
+
+        psU8 *rowPixels = sf->diffedPixels->data.U8[y];
+
+        int n = xright - xleft + 1;
+        if (n > 0) {
+            memset(&rowPixels[xleft], 255, n);
+        }
+    }
+}
+
+// x as a function of y for the line between two points
+// Note: the caller guarentees that the y's of the two points are different
+static double xOfY(psPlane *pI, psPlane *pJ, int y)
+{
+    double m_inv = (double) (pJ->x - pI->x) / (double) (pJ->y - pI->y);
+
+    double x = ((double) y - pI->y) * m_inv + (double) pI->x;
+
+    return x;
+}
+
+static int xLeft(psPlane *pt, int y)
+{
+    double x_d;
+
+    if ((y <= pt[BL].y ) && ((pt[BL].y >= pt[BR].y))) {
+        // left boundary is the line from bottom left to bottom right
+        x_d = xOfY(&pt[BL], &pt[BR], y);
+    } else if ( y <= pt[TL].y ) {
+        // left boundary is the line from Bottom left to top left
+        x_d = xOfY(&pt[BL], &pt[TL], y);
+    } else {
+        // left boundary is the line from top left to top right
+        x_d= xOfY(&pt[TL], &pt[TR], y);
+    }
+
+    return (int) floor(x_d);
+}
+
+static int xRight(psPlane *pt, int y)
+{
+    double x_d;
+    if ((y <= pt[BR].y) && (pt[BR].y > pt[BL].y)) {
+        // right boundary is the line from Bottom right to Top right
+        x_d = xOfY(&pt[BL], &pt[BR], y);
+    } else if ( y<= pt[TR].y) {
+        // right boundary is the line from Top left to top right
+        x_d = xOfY(&pt[BR], &pt[TR], y);
+    } else {
+        x_d = xOfY(&pt[TL], &pt[TR], y);
+    }
+    return (int) ceil(x_d);
+}
+
+static void
+nameCorners(psPlane pt[4])
+{
+    // sort points top to bottom (in case of ties take favor rightmost)
+    int i;
+    int j;
+    for (i=0; i<4; i++) {
+        for (j = i + 1; j < 4; j++) {
+            if ((pt[j].y > pt[i].y) || 
+                ((pt[j].y == pt[i].y) && (pt[j].x > pt[i].x))) {
+                psPlane tmpPt = pt[j];
+                pt[j] = pt[i];
+                pt[i] = tmpPt;
+            }
+        }
+    }
+    psPlane bl, br;
+    // sort bottom two points in x
+    if (pt[3].x < pt[2].x) {
+        bl = pt[3];
+        br = pt[2];
+    } else {
+        bl = pt[2];
+        br = pt[3];
+    }
+    // sort the top two points in x
+    psPlane tl, tr;
+    if (pt[0].x < pt[1].x) {
+        tl = pt[0];
+        tr = pt[1];
+    } else {
+        tl = pt[1];
+        tr = pt[0];
+    }
+    pt[BL] = bl;
+    pt[BR] = br;
+    pt[TL] = tl;
+    pt[TR] = tr;
+}
Index: branches/tap_branches/magic/remove/src/isdestreaked.c
===================================================================
--- branches/tap_branches/magic/remove/src/isdestreaked.c	(revision 25900)
+++ branches/tap_branches/magic/remove/src/isdestreaked.c	(revision 27838)
@@ -1,3 +1,5 @@
 #include "streaksremove.h"
+
+char *streaksProgram = "isdestreaked";
 
 int
Index: branches/tap_branches/magic/remove/src/streaksVersion.c
===================================================================
--- branches/tap_branches/magic/remove/src/streaksVersion.c	(revision 27838)
+++ branches/tap_branches/magic/remove/src/streaksVersion.c	(revision 27838)
@@ -0,0 +1,113 @@
+/** @file streaksVersion.c
+ *
+ *  @brief
+ *
+ *  @ingroup streaks
+ *
+ *  Copyright 2009 Institute for Astronomy, University of Hawaii
+ */
+
+#include "streaksremove.h"
+#include "streaksVersionDefinitions.h"
+
+#ifndef STREAKS_VERSION
+#error "STREAKS_VERSION is not set"
+#endif
+#ifndef STREAKS_BRANCH
+#error "STREAKS_BRANCH is not set"
+#endif
+#ifndef STREAKS_SOURCE
+#error "STREAKS_SOURCE is not set"
+#endif
+
+psString streaksVersion(void)
+{
+    char *value = NULL;
+    psStringAppend(&value, "%s@%s", STREAKS_BRANCH, STREAKS_VERSION);
+    return value;
+}
+
+psString streaksSource(void)
+{
+  return psStringCopy(STREAKS_SOURCE);
+}
+
+psString streaksVersionLong(void)
+{
+    psString version = streaksVersion();  // Version, to return
+    psString source = streaksSource();    // Source
+
+    psStringPrepend(&version, "%s ", streaksProgram);
+    psStringAppend(&version, " from %s, built %s, %s", source, __DATE__, __TIME__);
+    psFree(source);
+
+#ifdef __OPTIMIZE__
+    psStringAppend(&version, " optimised");
+#else
+    psStringAppend(&version, " unoptimised");
+#endif
+
+    return version;
+};
+
+bool streaksVersionHeader(psMetadata *header)
+{
+    PS_ASSERT_METADATA_NON_NULL(header, false);
+
+    psString version = streaksVersion(); // Software version
+    psString source = streaksSource();   // Software source
+
+    psStringPrepend(&version, "%s version: ", streaksProgram);
+    psStringPrepend(&source, "%s source: ", streaksProgram);
+
+    psMetadataAddStr(header, PS_LIST_TAIL, "HISTORY", PS_META_DUPLICATE_OK, NULL, version);
+    psMetadataAddStr(header, PS_LIST_TAIL, "HISTORY", PS_META_DUPLICATE_OK, NULL, source);
+
+    psFree(version);
+    psFree(source);
+
+    return true;
+}
+
+bool streaksVersionHeaderFull(psMetadata *header)
+{
+    PS_ASSERT_METADATA_NON_NULL(header, false);
+
+    psTime *time = psTimeGetNow(PS_TIME_TAI); // The time now
+    psString timeString = psTimeToISO(time); // The time in an ISO string
+    psFree(time);
+    psString history = NULL;               // History string
+    psStringAppend(&history, "%s at %s", streaksProgram, timeString);
+    psFree(timeString);
+    psMetadataAddStr(header, PS_LIST_TAIL, "HISTORY", PS_META_DUPLICATE_OK, NULL, history);
+    psFree(history);
+
+    psLibVersionHeader(header);
+    psModulesVersionHeader(header);
+    streaksVersionHeader(header);
+
+    return true;
+}
+
+void streaksVersionPrint(void)
+{
+    psTime *time = psTimeGetNow(PS_TIME_TAI); // The time now
+    psString timeString = psTimeToISO(time); // The time in an ISO string
+    psFree(time);
+    psLogMsg("streaks", PS_LOG_INFO, "%s at %s", streaksProgram, timeString);
+    psFree(timeString);
+
+    psString pslib = psLibVersionLong();// psLib version
+    psString psmodules = psModulesVersionLong(); // psModules version
+    psString streaks = streaksVersionLong(); // streaks version
+
+    psLogMsg("streaks", PS_LOG_INFO, "%s", pslib);
+    psLogMsg("streaks", PS_LOG_INFO, "%s", psmodules);
+    psLogMsg("streaks", PS_LOG_INFO, "%s", streaks);
+
+    psFree(pslib);
+    psFree(psmodules);
+    psFree(streaks);
+
+    return;
+}
Index: branches/tap_branches/magic/remove/src/streaksVersionDefinitions.h.in
===================================================================
--- branches/tap_branches/magic/remove/src/streaksVersionDefinitions.h.in	(revision 27838)
+++ branches/tap_branches/magic/remove/src/streaksVersionDefinitions.h.in	(revision 27838)
@@ -0,0 +1,8 @@
+#ifndef STREAKS_VERSION_DEFINITIONS_H
+#define STREAKS_VERSION_DEFINITIONS_H
+
+#define STREAKS_VERSION @STREAKS_VERSION@ // SVN version
+#define STREAKS_BRANCH  @STREAKS_BRANCH@  // SVN branch
+#define STREAKS_SOURCE  @STREAKS_SOURCE@  // SVN source
+
+#endif
Index: branches/tap_branches/magic/remove/src/streaksastrom.c
===================================================================
--- branches/tap_branches/magic/remove/src/streaksastrom.c	(revision 25900)
+++ branches/tap_branches/magic/remove/src/streaksastrom.c	(revision 27838)
@@ -378,7 +378,6 @@
 linearizeTransforms(strkAstrom *astrom)
 {
-    if (!pmAstromLinearizeTransforms((pmFPA *) astrom->fpa, (pmChip *) astrom->chip)) {
-        // streaksExit("linear fit to astrometry failed\n", PS_EXIT_UNKNOWN_ERROR);
-        psError(PS_ERR_UNKNOWN, false, "linear fit to astrometry failed ignoring");
+    if (!pmAstromLinearizeTransforms((pmFPA *) astrom->fpa, (pmChip *) astrom->chip, NULL, NULL, NULL, 0, 0)) {
+        psErrorStackPrint(stderr, "linear fit to astrometry failed. ignoring\n");
         return false;
     }
Index: branches/tap_branches/magic/remove/src/streakscompare.c
===================================================================
--- branches/tap_branches/magic/remove/src/streakscompare.c	(revision 25900)
+++ branches/tap_branches/magic/remove/src/streakscompare.c	(revision 27838)
@@ -1,3 +1,5 @@
 #include "streaksremove.h"
+
+char *streaksProgram = "streakscompare";
 
 static pmConfig * parseArguments(int, char **);
Index: branches/tap_branches/magic/remove/src/streaksextern.c
===================================================================
--- branches/tap_branches/magic/remove/src/streaksextern.c	(revision 25900)
+++ branches/tap_branches/magic/remove/src/streaksextern.c	(revision 27838)
@@ -78,5 +78,5 @@
         }
     }
-    printf("%d streaks on this component\n", streaksOnComponent);
+    fprintf(stderr, "%d streaks on this component\n", streaksOnComponent);
     return pixels;
 }
Index: branches/tap_branches/magic/remove/src/streaksio.c
===================================================================
--- branches/tap_branches/magic/remove/src/streaksio.c	(revision 25900)
+++ branches/tap_branches/magic/remove/src/streaksio.c	(revision 27838)
@@ -121,4 +121,14 @@
     sf->transparentStreaks = psMetadataLookupF64(&status, config->arguments, "TRANSPARENT_STREAKS");
 
+    sf->stats = psMetadataAlloc();
+    psString statsFileName= psMetadataLookupStr(&status, config->arguments, "STATS");
+    if (statsFileName) {
+        sf->statsFile = fopen(statsFileName, "w");
+        if (!sf->statsFile) {
+            psError(PS_ERR_IO, true, "failed to open stats file %s", statsFileName);
+            streaksExit("", PS_EXIT_CONFIG_ERROR);
+        }
+    }
+
     return sf;
 }
@@ -128,5 +138,5 @@
 {
     freeImages(sf);
-    psFree(sf->warpedPixels);
+    psFree(sf->diffedPixels);
     psFree(sf->tiles);
     psFree(sf->view);
@@ -800,4 +810,6 @@
         return;
     }
+    streaksVersionHeaderFull(sfile->header);
+
     if (!psFitsWriteImage(sfile->fits, sfile->header, sfile->image, 0, extname)) {
         psError(PS_ERR_IO, false, "failed to write image to %s extnum: %d",
@@ -817,4 +829,5 @@
         return;
     }
+    streaksVersionHeaderFull(sfile->header);
     if (!psFitsWriteImageCube(sfile->fits, sfile->header, imagecube, extname)) {
         psError(PS_ERR_IO, false, "failed to write image to %s extnum: %d",
@@ -851,7 +864,7 @@
     sfile->header = NULL;
     psFree(sfile->image);
-    sfile->header = NULL;
+    sfile->image = NULL;
     psFree(sfile->imagecube);
-    sfile->header = NULL;
+    sfile->imagecube = NULL;
 }
 
@@ -859,4 +872,14 @@
 closeImages(streakFiles *sf)
 {
+    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);
+        }
+    }
     closeImage(sf->inImage);
     closeImage(sf->outImage);
@@ -913,6 +936,6 @@
     bool free_user_copies = true;
     if (user_copies == NULL) {
-        user_copies = "2";
-        free_user_copies = false;
+        // input does not have replication requested
+        return true;
     }
     if (!nebSetXattr(server, outFile->name, "user.copies", user_copies,  NEB_REPLACE)) {
@@ -920,5 +943,5 @@
         return false;
     }
-    if (!nebReplicate(server, outFile->name, "any", NULL)) {
+    if (!nebReplicate(server, outFile->name, NULL, NULL)) {
         psError(PM_ERR_UNKNOWN, true, "nebReplicate failed for %s\n%s", outFile->name, nebErr(server));
         return false;
@@ -1111,10 +1134,16 @@
 
         psImage *image = sfiles->outImage->image;
+        int imageType = image->type.type;
+        double exciseValue = sfiles->inImage->exciseValue;
+
         psImage *mask = sfiles->inMask->image;
         psImage *weight = NULL;
+        int weightType = 0;
+        double weightExciseValue = NAN;
         if (sfiles->outWeight) {
             weight = sfiles->outWeight->image;
-        }
-        double exciseValue = sfiles->inImage->exciseValue;
+            weightType = weight->type.type;
+            weightExciseValue = sfiles->inWeight->exciseValue;
+        }
 
         if (printCounts) {
@@ -1126,6 +1155,17 @@
                 // these gets are not necessary, we could just set the pixels to nan
                 // but I want to get the counts
-                double imageVal  = image->data.F32[y][x];
-                psU32 maskVal;
+                double imageVal = 0; // avoid compiler warning
+                if (imageType == PS_TYPE_F32) {
+                    imageVal  = image->data.F32[y][x];
+                } else if (imageType == PS_TYPE_S16) {
+                    imageVal = image->data.S16[y][x];
+                } else if (imageType == PS_TYPE_U16) {
+                    imageVal = image->data.U16[y][x];
+                } else {
+                    psError(PS_ERR_PROGRAMMING, true, "unexpected image type found: %d\n",
+                                imageType);
+                    streaksExit("", PS_EXIT_PROG_ERROR);
+                }
+                psU32 maskVal = 0;
                 if (sfiles->stage == IPP_STAGE_RAW) {
                     unsigned int xChip, yChip;
@@ -1137,13 +1177,41 @@
                 if (maskVal & maskMask) {
                     ++maskedPixels;
-                    if (!isExciseValue(imageVal, sfiles->inImage->exciseValue)) {
+                    if (!isExciseValue(imageVal, exciseValue)) {
                         ++nandPixels;
-                        image->data.F32[y][x] = exciseValue;
+                        if (imageType == PS_TYPE_F32) {
+                            image->data.F32[y][x] = exciseValue;
+                        } else if (imageType == PS_TYPE_S16) {
+                            image->data.S16[y][x] = exciseValue;
+                        } else if (imageType == PS_TYPE_U16) {
+                            image->data.U16[y][x] = exciseValue;
+                        } else {
+                            psError(PS_ERR_PROGRAMMING, true, "unexpected image type found: %d\n",
+                                imageType);
+                            streaksExit("", PS_EXIT_PROG_ERROR);
+                        }
                     }
                     if (weight) {
-                        double weightVal = weight ? weight->data.F32[y][x] : 0;
-                        if (!isnan(weightVal)) {
-                            ++nandWeights;
-                            weight->data.F32[y][x] = NAN;
+                        if (weightType == PS_TYPE_F32) {
+                            double weightVal = weight->data.F32[y][x];
+                            if (!isnan(weightVal)) {
+                                ++nandWeights;
+                                weight->data.F32[y][x] = NAN;
+                            }
+                        } else if(weightType == PS_TYPE_S16) {
+                            double weightVal = weight->data.S16[y][x];
+                            if (!isExciseValue(weightVal, weightExciseValue)) {
+                                ++nandWeights;
+                                image->data.S16[y][x] = weightExciseValue;
+                            }
+                        } else if (weightType == PS_TYPE_U16) {
+                            double weightVal = weight->data.U16[y][x];
+                            if (!isExciseValue(weightVal, weightExciseValue)) {
+                                ++nandWeights;
+                                image->data.U16[y][x] = weightExciseValue;
+                            }
+                        } else {
+                            psError(PS_ERR_PROGRAMMING, true, "unexpected image type found: %d\n",
+                                    weightType);
+                                streaksExit("", PS_EXIT_PROG_ERROR);
                         }
                     }
@@ -1176,12 +1244,26 @@
         psMetadataAddU16(in->header, PS_LIST_TAIL, "BLANK", 0, "", 65535);
         psMetadataAddU16(in->header, PS_LIST_TAIL, "ZBLANK", 0, "", 65535);
-    } else {
+    } else if (in->image->type.type == PS_TYPE_S16) {
+        in->exciseValue = 32767;
+        psMetadataAddU16(in->header, PS_LIST_TAIL, "BLANK", 0, "", 32767);
+        psMetadataAddU16(in->header, PS_LIST_TAIL, "ZBLANK", 0, "", 32767);
+    } else if (in->image->type.type == PS_TYPE_F32) {
         in->exciseValue = NAN;
-    }
-}
-
-void
-strkGetMaskValues(streakFiles *sfiles, psU32 *maskStreak, psU32 *maskMask)
-{
+    } else {
+        psError(PS_ERR_PROGRAMMING, true, "unexpected image type found: %d\n", in->image->type.type);
+        streaksExit("", PS_EXIT_PROG_ERROR);
+    }
+}
+
+// Get the mask values that we need.
+// If an input mask file is provided get them from there.
+// Otherwise get them from the recipe MASKS
+void
+strkGetMaskValues(streakFiles *sfiles)
+{
+    if (sfiles->maskStreak) {
+        // already done
+        return;
+    }
     if (sfiles->inMask && sfiles->inMask->header) {
         if (!pmConfigMaskReadHeader(sfiles->config, sfiles->inMask->header)) {
@@ -1196,5 +1278,5 @@
         streaksExit("", PS_EXIT_CONFIG_ERROR);
     }
-    *maskStreak = psMetadataLookupU32(&status, masks, "STREAK");
+    sfiles->maskStreak = psMetadataLookupU32(&status, masks, "STREAK");
     if (!status) {
         psError(PM_ERR_CONFIG, false, "failed to lookup mask value for STREAK in recipes\n");
@@ -1212,5 +1294,5 @@
         }
     }
-    *maskMask = ~convPoor;
+    sfiles->maskMask = ~convPoor;
 }
 
Index: branches/tap_branches/magic/remove/src/streaksio.h
===================================================================
--- branches/tap_branches/magic/remove/src/streaksio.h	(revision 25900)
+++ branches/tap_branches/magic/remove/src/streaksio.h	(revision 27838)
@@ -16,5 +16,5 @@
 void copyFitsOptions(sFile *out, sFile *rec, sFile *in, psVector *tiles);
 void setupImageRefs(sFile *out, sFile *recoveryOut, sFile *in, int extnum, bool exciseAll);
-void strkGetMaskValues(streakFiles *sfiles, psU32 *maskStreak, psU32 *maskMask);
+void strkGetMaskValues(streakFiles *sfiles);
 void copyExtraExtensions(streakFiles *sfiles);
 
Index: branches/tap_branches/magic/remove/src/streaksrelease.c
===================================================================
--- branches/tap_branches/magic/remove/src/streaksrelease.c	(revision 25900)
+++ branches/tap_branches/magic/remove/src/streaksrelease.c	(revision 27838)
@@ -4,4 +4,6 @@
 static bool readAndCopyToOutput(streakFiles *sf, bool exciseAll);
 static void writeImages(streakFiles *sf, bool exciseImageCube);
+
+char *streaksProgram = "streaksrelease";
 
 int
@@ -16,8 +18,4 @@
         return PS_EXIT_CONFIG_ERROR;
     }
-
-    // Values to set for masked pixels
-    psU32 maskStreak = 0;           // for the image and weight (usually NAN, MAXINT for integer images)
-    psU32 maskMask = 0;             // value looked up for MASK.STREAK
 
     // Does true work here?
@@ -50,9 +48,9 @@
 
         // now that we've read the input files, lookup the mask values that we read
-        if (maskStreak == 0) {
-            strkGetMaskValues(sfiles, &maskStreak, &maskMask);
-        }
-
-        setMaskedToNAN(sfiles, maskMask, true);
+        strkGetMaskValues(sfiles);
+
+        if (!sfiles->inImage->imagecube) {
+            setMaskedToNAN(sfiles, sfiles->maskMask, true);
+        }
 
         // write out the destreaked temporary images and the recovery images
@@ -226,9 +224,6 @@
     if (sf->inImage->image) {
         setupImageRefs(sf->outImage, sf->recImage, sf->inImage, sf->extnum, exciseAll);
-    } else if (sf->inImage->imagecube) {
-        // Image cubes should have been excised in the destreaking process
-        streaksExit("unexpected imagecube found", PS_EXIT_CONFIG_ERROR);
     }  else {
-        return false;
+        // image cubes are handled specially
     }
 
@@ -278,4 +273,5 @@
 
 
+// XXXX: why isn't this in streaksio.c ?? See also streakremove.c
 static void
 writeImages(streakFiles *sf, bool exciseImageCube)
@@ -303,18 +299,4 @@
         }
 
-        // 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;
Index: branches/tap_branches/magic/remove/src/streaksremove.c
===================================================================
--- branches/tap_branches/magic/remove/src/streaksremove.c	(revision 25900)
+++ branches/tap_branches/magic/remove/src/streaksremove.c	(revision 27838)
@@ -12,12 +12,18 @@
 static pmConfig *parseArguments(int argc, char **argv);
 static bool readAndCopyToOutput(streakFiles *sf, bool exciseAll);
-static void exciseNonWarpedPixels(streakFiles *sfiles, psImageMaskType newMaskValue);
-static bool warpedPixel(streakFiles *sfiles, int x, int y);
+static long exciseNonDiffedPixels(streakFiles *sfiles, psImageMaskType newMaskValue);
+static bool diffedPixel(streakFiles *sfiles, int x, int y);
 static void excisePixel(streakFiles *sfiles, unsigned int x, unsigned int y, bool streak, psImageMaskType newMaskValue);
 static void writeImages(streakFiles *sf, bool exciseImageCube);
 static void updateAstrometry(streakFiles *sfiles);
 static void censorSources(streakFiles *sfiles, psImageMaskType maskStreak);
-static long censorPixels(streakFiles *sfiles, psImage * pixels, bool checkNonWarpedPixels, psU16 maskStreak);
-
+static long censorPixels(streakFiles *sfiles, psImage * pixels, bool checkNonDiffedPixels, psU16 maskStreak);
+
+char *streaksProgram = "streaksremove";
+
+// Note: For clarity the flow of this program is in main().
+// There is not a lot of error checking is done in main.
+// Until the end, where we might be doing Nebulous operations, called functions exit when an error
+// is encountered.
 int
 main(int argc, char *argv[])
@@ -26,5 +32,5 @@
 
     psLibInit(NULL);
-    psTimerStart("STREAKSREMOVE");
+    psTimerStart("TOTAL_TIME");
 
     pmConfig *config = parseArguments(argc, argv);
@@ -33,8 +39,4 @@
         streaksExit("", PS_EXIT_CONFIG_ERROR);
     }
-
-    // Values to set for masked pixels
-    psU32 maskStreak = 0;           // for the image and weight (usually NAN, MAXINT for integer images)
-    psU32 maskMask = 0;             // value looked up for MASK.STREAK
 
     psString streaksFileName = psMetadataLookupStr(NULL, config->arguments, "STREAKS");
@@ -50,4 +52,5 @@
     streakFiles *sfiles = openFiles(config, true, argv[0]);
     setupAstrometry(sfiles);
+    sfiles->stats = psMetadataAlloc();
 
     // Optionally we can set pixels that are masked to NAN since they couldn't have been
@@ -61,12 +64,12 @@
 
     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 ) {
+    // --keepnondiffed is a test and debug mode
+    bool keepNonDiffedPixels = psMetadataLookupBool(&status, config->arguments, "KEEP_NON_DIFFED");
+
+    // we need to check for non diffed pixels unless we've been asked not to or the stage is diff
+    // (By definition pixels in diff images were included in the difference images)
+    bool checkNonDiffedPixels = ! (keepNonDiffedPixels || (sfiles->stage == IPP_STAGE_DIFF) );
+
+    if (checkNonDiffedPixels ) {
         // From magic ICD:
         // In the raw and detrended images, the pixels which were not
@@ -76,11 +79,12 @@
         // if no skycells are provided sfiles->exciseAll is set to true
 
-        psTimerStart("COMPUTE_WARPED_PIXELS");
-        if (! computeWarpedPixels(sfiles) ) {
+        psTimerStart("COMPUTE_DIFFED_PIXELS");
+        if (! computeDiffedPixels(sfiles) ) {
             // we have no choice to excise all pixels
             exciseAll = true;
         }
-        psF64 cwp_t = psTimerClear("COMPUTE_WARPED_PIXELS");
-        psLogMsg("streaksremove", PS_LOG_INFO, "time to compute warped pixels: %f\n", cwp_t);
+        psF64 cdp_t = psTimerClear("COMPUTE_DIFFED_PIXELS");
+        psMetadataAddF32(sfiles->stats, PS_LIST_TAIL, "COMPUTE_UNDIFFED_PIXELS", PS_META_REPLACE, "time to compute non-diffedpixels", cdp_t);
+        psLogMsg("streaksremove", PS_LOG_INFO, "time to compute diffed pixels: %f\n", cdp_t);
     }
 
@@ -100,5 +104,13 @@
     long totalPixels = 0;
     long totalStreakPixels = 0;
-
+    long nonDiffedPixels = 0;
+
+    // accumulators for the various timers
+    psF64 gsp_t = 0;
+    psF64 enw_t = 0;
+    psF64 rms_t = 0;
+    psF64 cs_t = 0;
+    psF64 wi_t = 0;
+    psF64 ua_t = 0;
     // Iterate through each component of the input (except for raw images there is only one)
     do {
@@ -118,7 +130,5 @@
 
         // now that we've read the input files, lookup the mask values
-        if (maskStreak == 0) {
-            strkGetMaskValues(sfiles, &maskStreak, &maskMask);
-        }
+        strkGetMaskValues(sfiles);
 
         totalPixels += sfiles->inImage->numRows * sfiles->inImage->numCols;
@@ -132,4 +142,6 @@
             StreakPixels *pixels = streak_on_component(streaks, sfiles->astrom, sfiles->inImage->numCols,
                                                         sfiles->inImage->numRows);
+            gsp_t +=  psTimerClear("GET_STREAK_PIXELS");
+            psMetadataAddF32(sfiles->stats, PS_LIST_TAIL, "GET_STREAK_PIXELS", PS_META_REPLACE, "", gsp_t);
             psLogMsg("streaksremove", PS_LOG_INFO, "time to get streak pixels: %f\n", psTimerClear("GET_STREAK_PIXELS"));
 
@@ -137,23 +149,27 @@
             // otherwise it contained an image cube (video cell) which is handled in the if block
             if (sfiles->inImage->image) {
-                if (checkNonWarpedPixels) {
-                    psTimerStart("EXCISE_NON_WARPED");
-
-                    // set non-warped pixels and variance to NAN, mask to maskStreak (since the pixel
+                if (checkNonDiffedPixels) {
+                    psTimerStart("EXCISE_NON_DIFFED");
+
+                    // set non-diffed pixels and variance to NAN, mask to maskStreak (since the pixel
                     // is excised as part of the destreaking process)
-                    exciseNonWarpedPixels(sfiles, maskStreak);
-
-                    psLogMsg("streaksremove", PS_LOG_INFO, "time to excise non warped pixels: %f\n", psTimerClear("EXCISE_NON_WARPED"));
+                    nonDiffedPixels += exciseNonDiffedPixels(sfiles, sfiles->maskStreak);
+
+                    enw_t +=  psTimerClear("EXCISE_NON_DIFFED");
+                    psMetadataAddF32(sfiles->stats, PS_LIST_TAIL, "EXCISE_NON_DIFFED", PS_META_REPLACE, "", enw_t);
+                    psLogMsg("streaksremove", PS_LOG_INFO, "time to excise non diffed pixels: %f\n", enw_t);
                 }
 
                 psTimerStart("REMOVE_STREAKS");
 
-                totalStreakPixels += censorPixels(sfiles, pixels, checkNonWarpedPixels, maskStreak);
-
-                psLogMsg("streaksremove", PS_LOG_INFO, "time to remove streak pixels: %f\n", psTimerClear("REMOVE_STREAKS"));
+                totalStreakPixels += censorPixels(sfiles, pixels, checkNonDiffedPixels, sfiles->maskStreak);
+
+                rms_t += psTimerClear("REMOVE_STREAKS");
+                psMetadataAddF32(sfiles->stats, PS_LIST_TAIL, "REMOVE_STREAKS", PS_META_REPLACE, "", enw_t);
+                psLogMsg("streaksremove", PS_LOG_INFO, "time to remove streak pixels: %f\n", rms_t);
 
                 if (nanForRelease) {
                     // set any pixels that were masked, to NAN (unless they are already NAN)
-                    setMaskedToNAN(sfiles, maskMask, true);
+                    setMaskedToNAN(sfiles, sfiles->maskMask, true);
                 }
 
@@ -170,14 +186,24 @@
             // chip processed files with the data calcuated by psastro at the camera stage
             // (actually we use a linear approximation)
+            psTimerStart("UPDATE_ASTROMETRY");
             updateAstrometry(sfiles);
-        }
-
-        censorSources(sfiles, maskStreak);
+            ua_t += psTimerClear("UPDATE_ASTROMETRY");
+            psMetadataAddF32(sfiles->stats, PS_LIST_TAIL, "UPDATE_ASTROMETRY", PS_META_REPLACE, "", ua_t);
+            psLogMsg("streaksremove", PS_LOG_INFO, "time to update astrometry: %f\n", ua_t);
+        }
+
+        psTimerStart("CENSOR_SOURCES");
+        censorSources(sfiles, sfiles->maskStreak);
+        cs_t += psTimerClear("CENSOR_SOURCES");
+        psMetadataAddF32(sfiles->stats, PS_LIST_TAIL, "CENSOR_SOURCES", PS_META_REPLACE, "", cs_t);
 
         // write the destreaked "temporary" images and the recovery images
+        psTimerStart("WRITE_IMAGES");
         writeImages(sfiles, exciseImageCube);
-
+        wi_t += psTimerClear("WRITE_IMAGES");
+        psMetadataAddF32(sfiles->stats, PS_LIST_TAIL, "WRITE_IMAGES", PS_META_REPLACE, "", wi_t);
 
         psLogMsg("streaksremove", PS_LOG_INFO, "time to process component %d: %f\n", sfiles->extnum, psTimerClear("PROCESS_COMPONENT"));
+
     } while (streakFilesNextExtension(sfiles));
 
@@ -185,5 +211,16 @@
     psFree(streaks);
 
-    psLogMsg("streaksremove", PS_LOG_INFO, "pixels: %ld streak pixels: %ld %4.2f%%\n", totalPixels, totalStreakPixels, 100. * totalStreakPixels / totalPixels);
+    if (exciseAll) {
+        totalStreakPixels = totalPixels;
+    }
+
+    psF64 streakFraction = (double) totalStreakPixels / totalPixels;
+    psLogMsg("streaksremove", PS_LOG_INFO, "   total pixels:  %8ld\n", totalPixels);
+    psLogMsg("streaksremove", PS_LOG_INFO, "  streak pixels:  %8ld %4.2f%%\n", totalStreakPixels, streakFraction * 100);
+    psMetadataAddF32(sfiles->stats, PS_LIST_TAIL, "STREAK_FRACTION", PS_META_REPLACE, "", streakFraction);
+
+    psF64 nonDiffedFraction = (double) nonDiffedPixels / totalPixels;
+    psLogMsg("streaksremove", PS_LOG_INFO, "nondiffed pixels:  %8ld %4.2f%%\n", nonDiffedPixels, nonDiffedFraction * 100);
+    psMetadataAddF32(sfiles->stats, PS_LIST_TAIL, "NONDIFFED_FRACTION", PS_META_REPLACE, "", nonDiffedFraction);
 
     // check the weight and mask files for extra extensions that might be in files
@@ -194,10 +231,12 @@
 
     psTimerStart("CLOSE_IMAGES");
-
     closeImages(sfiles);
-
-    psLogMsg("streaksremove", PS_LOG_INFO, "time to close images: %f\n", psTimerClear("CLOSE_IMAGES"));
-
-
+    psF64 ci_t = psTimerClear("CLOSE_IMAGES");
+    psMetadataAddF32(sfiles->stats, PS_LIST_TAIL, "CLOSE_IMAGES", PS_META_REPLACE, "", ci_t);
+
+    psLogMsg("streaksremove", PS_LOG_INFO, "time to close images: %f\n", ci_t);
+
+#ifdef DO_REPLICATE
+    psTimerStart("REPLICATE_OUTPUTS");
     if (!replicateOutputs(sfiles)) {
         psErrorStackPrint(stderr, "failed to replicate output files");
@@ -205,4 +244,7 @@
         exit(PS_EXIT_UNKNOWN_ERROR);
     }
+    psF64 ro_t = psTimerClear("REPLICATE_OUTPUTS");
+    psMetadataAddF32(sfiles->stats, PS_LIST_TAIL, "REPLICATE_OUTPUTS", PS_META_REPLACE, "", ro_t);
+#endif
 
     // NOTE: from here on we can't just quit if something goes wrong.
@@ -214,9 +256,9 @@
         //     swap the instances for the input and output
         //     Note this is a nebulous database operation. No file I/O is performed
+        psTimerStart("SWAP_INSTANCES");
         if (!swapOutputsToInputs(sfiles)) {
-            // XXX: Now what?
             // It is up to the program that reverts failed destreak runs to insure that
-            // any input files that have been swapped are restored and that the de-streaked
-            // versions are deleted
+            // any original non-destreaked input files that have been swapped are restored and that the de-streaked
+            // versions are deleted.
 
             psErrorStackPrint(stderr, "failed to swap files");
@@ -224,4 +266,24 @@
             // XXX: pick a specific error code for this failure
             exit(PS_EXIT_UNKNOWN_ERROR);
+        }
+        psF64 si_t = psTimerClear("SWAP_INSTANCES");
+        psMetadataAddF32(sfiles->stats, PS_LIST_TAIL, "SWAP_INSTANCES", PS_META_REPLACE, "", si_t);
+    }
+
+    psF64 total_time = psTimerClear("TOTAL_TIME");
+    psMetadataAddF32(sfiles->stats, PS_LIST_TAIL, "TOTAL_TIME", PS_META_REPLACE, "", total_time);
+    psLogMsg("streaksremove", PS_LOG_INFO, "time to run streaksremove: %f\n", total_time);
+
+    if (sfiles->statsFile) {
+        const char *statsMDC = psMetadataConfigFormat(sfiles->stats);
+        if (!statsMDC || strlen(statsMDC) == 0) {
+            psError(PS_ERR_IO, false, "Unable to get statistics MDC file.\n");
+        } else {
+            fprintf(sfiles->statsFile, "%s", statsMDC);
+            psFree(statsMDC);
+            fclose(sfiles->statsFile);
+            sfiles->statsFile = NULL;
+            psFree(sfiles->stats);
+            sfiles->stats = NULL;
         }
     }
@@ -234,5 +296,4 @@
     streaksNebulousCleanup();
     pmConfigDone();
-    psLogMsg("streaksremove", PS_LOG_INFO, "time to run streaksremove: %f\n", psTimerClear("STREAKSREMOVE"));
     psLibFinalize();
 
@@ -243,5 +304,5 @@
 
 static long
-censorPixels(streakFiles *sfiles, psImage *pixels, bool checkNonWarpedPixels, psU16 maskStreak)
+censorPixels(streakFiles *sfiles, psImage *pixels, bool checkNonDiffedPixels, psU16 maskStreak)
 {
     long streakPixels = 0;
@@ -250,6 +311,6 @@
         for (int x = 0; x < sfiles->inImage->numCols; x++) {
             if (psImageGet(pixels, x, y)) {
-                ++streakPixels;
-                if (!checkNonWarpedPixels || warpedPixel(sfiles, x, y)) {
+                if (!checkNonDiffedPixels || diffedPixel(sfiles, x, y)) {
+                    ++streakPixels;
 
                     excisePixel(sfiles, x, y, true, maskStreak);
@@ -257,5 +318,5 @@
                 } else {
                     // This pixel was not included in any warp and has thus already excised
-                    // by exciseNonWarpedPixels
+                    // by exciseNonDiffedPixels
                 }
             }
@@ -274,5 +335,5 @@
     fprintf(stderr, "\t-weight WEIGHT.fits: weight file to de-streak\n");
     fprintf(stderr, "\t-replace: replace the input images with the output\n");
-    fprintf(stderr, "\t-keepnonwarped: do not exise pixels that were not part of difference processing\n");
+    fprintf(stderr, "\t-keepnondiffed: do not exise pixels that were not part of difference processing\n");
     fprintf(stderr, "\t-transparent val: instead of setting excicsed pixel to NAN add val\n");
     fprintf(stderr, "\t-chip_mask MASK.fits: name of mask for chip stage (camera stage mask is passed tih -mask)\n");
@@ -365,8 +426,8 @@
     }
 
-    if ((argnum = psArgumentGet(argc, argv, "-keepnonwarped"))) {
-        psArgumentRemove(argnum, &argc, argv);
-        psMetadataAddBool(config->arguments, PS_LIST_TAIL, "KEEP_NON_WARPED", 0,
-            "skip excising of non warped pixels", true);
+    if ((argnum = psArgumentGet(argc, argv, "-keepnondiffed"))) {
+        psArgumentRemove(argnum, &argc, argv);
+        psMetadataAddBool(config->arguments, PS_LIST_TAIL, "KEEP_NON_DIFFED", 0,
+            "skip excising of non diffed pixels", true);
     }
 
@@ -378,10 +439,10 @@
         psArgumentRemove(argnum, &argc, argv);
         double transparentStreaks = atof(argv[argnum]);
-        psMetadataAddF64(config->arguments, PS_LIST_TAIL, "TRANSPARENT_STREAKS", 0,
+        psMetadataAddF32(config->arguments, PS_LIST_TAIL, "TRANSPARENT_STREAKS", 0,
             "value to adjust excised pixels", transparentStreaks);
         psArgumentRemove(argnum, &argc, argv);
     }
 
-    // if skycells are not provided then we have to execise all pixels  unless -keepnonwarped
+    // if skycells are not provided then we have to execise all pixels  unless -keepnondiffed
     pmConfigFileSetsMD(config->arguments, &argc, argv, "SKYCELLS", "-skycell", "-skycelllist");
 
@@ -449,4 +510,11 @@
     }
 
+    if ((argnum = psArgumentGet(argc, argv, "-stats"))) {
+        psArgumentRemove(argnum, &argc, argv);
+        psMetadataAddStr(config->arguments, PS_LIST_TAIL, "STATS", 0,
+                "name of input stats file", argv[argnum]);
+        psArgumentRemove(argnum, &argc, argv);
+    }
+
     if ((argnum = psArgumentGet(argc, argv, "-tmproot"))) {
         psArgumentRemove(argnum, &argc, argv);
@@ -491,21 +559,28 @@
 updateAstrometry(streakFiles *sf)
 {
-    // XXX: why do I check this here? Shouldn't it be just around the call to linearizeTransforms?
     if (sf->bilevelAstrometry) {
-
         if (!linearizeTransforms(sf->astrom)) {
-            // fit failed, leave the astrometry unchanged
+            // fit failed, leave the transform in the file unchanged
             return;
         }
-
-        if (!pmAstromWriteWCS(sf->outImage->header, sf->inAstrom->fpa, sf->chip, 0.001)) {
-            psError(PS_ERR_UNKNOWN, false, "failed to update astrometry for extension %d", sf->extnum);
-            streaksExit("", PS_EXIT_UNKNOWN_ERROR);
-        }
-        if (sf->outMask) {
-            pmAstromWriteWCS(sf->outMask->header, sf->inAstrom->fpa, sf->chip, 0.001);
-        }
-        if (sf->outWeight) {
-            pmAstromWriteWCS(sf->outWeight->header, sf->inAstrom->fpa, sf->chip, 0.001);
+    }
+    if (!pmAstromWriteWCS(sf->outImage->header, sf->inAstrom->fpa, sf->chip, 0.001)) {
+        psError(PS_ERR_UNKNOWN, false, "failed to update astrometry for extension %d", sf->extnum);
+        streaksExit("", PS_EXIT_UNKNOWN_ERROR);
+    }
+    if (sf->outMask) {
+        pmAstromWriteWCS(sf->outMask->header, sf->inAstrom->fpa, sf->chip, 0.001);
+    }
+    if (sf->outWeight) {
+        pmAstromWriteWCS(sf->outWeight->header, sf->inAstrom->fpa, sf->chip, 0.001);
+    }
+}
+
+static void 
+setStreakBits(psImage *maskImage, psU32 maskStreak)
+{
+    for (int y=0 ; y < maskImage->numRows; y++) {
+        for (int x=0 ; x < maskImage->numCols; x++) {
+            maskImage->data.PS_TYPE_IMAGE_MASK_DATA[y][x] |= maskStreak;
         }
     }
@@ -516,5 +591,4 @@
 readAndCopyToOutput(streakFiles *sf, bool exciseAll)
 {
-    bool    updateAstrometry = false;
     if (sf->inImage->pmfile) {
         // image data from pmFPAfile (diff or warp file)
@@ -534,15 +608,4 @@
         streaksExit("", PS_EXIT_UNKNOWN_ERROR);
     }
-    // For the chip level files, copy the WCS from the astrometry file to the header
-    // XXX: do we want to do this for raw images as well?
-    if (sf->stage == IPP_STAGE_CHIP) {
-        if (!sf->bilevelAstrometry) {
-            updateAstrometry = true;
-            if (!pmAstromWriteWCS(sf->inImage->header, sf->inAstrom->fpa, sf->chip, 0.001)) {
-                psError(PS_ERR_UNKNOWN, false, "failed to update astrometry for extension %d", sf->extnum);
-                streaksExit("", PS_EXIT_UNKNOWN_ERROR);
-            }
-        }
-    }
     sf->outImage->header =  psMemIncrRefCounter(sf->inImage->header);
     if (sf->recImage) {
@@ -577,8 +640,12 @@
             }
             addDestreakKeyword(sf->outMask->header);
-            if (updateAstrometry) {
-                pmAstromWriteWCS(sf->outMask->header, sf->inAstrom->fpa, sf->chip, 0.001);
-            }
-            setupImageRefs(sf->outMask, sf->recMask, sf->inMask, sf->extnum, exciseAll);
+            // Note: we don't excise the mask pixels even if exciseAll is true.
+            setupImageRefs(sf->outMask, sf->recMask, sf->inMask, sf->extnum, false);
+            if (exciseAll) {
+                strkGetMaskValues(sf);
+                
+                // add the STREAK bit to the mask image pixels
+                setStreakBits(sf->inMask->image, sf->maskStreak);
+            }
             if (sf->outChMask) {
                 sf->outChMask->header = (psMetadata *) psMemIncrRefCounter(sf->outMask->header);
@@ -605,7 +672,4 @@
         }
         addDestreakKeyword(sf->outWeight->header);
-        if (updateAstrometry) {
-            pmAstromWriteWCS(sf->inWeight->header, sf->inAstrom->fpa, sf->chip, 0.001);
-        }
         setupImageRefs(sf->outWeight, sf->recWeight, sf->inWeight, sf->extnum, exciseAll);
 
@@ -630,27 +694,31 @@
     } 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);
+        }
+
+        // Now deal with the other output image
         if (exciseImageCube) {
-            sf->outImage->image = image;
-            writeImage(sf->outImage, extname, sf->extnum);
+            // Set the values in the imagecube images to NAN and write them to the output image
+            for (int i = 0; i < psArrayLength(sf->inImage->imagecube); i++) {
+                psImage *image = psArrayGet (sf->inImage->imagecube, i);
+                // XXX: NAN isn't right. It should be the integer equivalent. Should use exciseValue
+                // but it isn't set with this code path. Fix that.
+                psImageInit(image, 65535);
+            }
+            writeImageCube(sf->outImage, sf->inImage->imagecube, extname, sf->extnum);
         } else {
-            // write zero valued image to reccovery
             if (sf->recImage) {
-                sf->recImage->image = image;
-                writeImage(sf->recImage, extname, sf->extnum);
+                // Set the values in the imagecube images to zero
+                for (int i = 0; i < psArrayLength(sf->inImage->imagecube); i++) {
+                    psImage *image = psArrayGet (sf->inImage->imagecube, i);
+                    psImageInit(image, 0);
+                }
+                // copy the entire zeroed image to the recovery image
+                writeImageCube(sf->recImage, sf->inImage->imagecube, extname, sf->extnum);
             }
         }
@@ -683,17 +751,35 @@
     }
 
-    float imageValue  = sfiles->inImage->image->data.F32[y][x];
-    if (sfiles->recImage && !isExciseValue(imageValue, sfiles->inImage->exciseValue) ) {
-        sfiles->recImage->image->data.F32[y][x] = imageValue;
-    }
-
-    if (sfiles->transparentStreaks == 0) {
-        sfiles->outImage->image->data.F32[y][x] = exciseValue;
+    if (sfiles->inImage->image->type.type == PS_TYPE_U16) {
+        psU16 imageValue  = sfiles->inImage->image->data.U16[y][x];
+        if (sfiles->recImage && !isExciseValue(imageValue, sfiles->inImage->exciseValue) ) {
+            sfiles->recImage->image->data.U16[y][x] = imageValue;
+        }
+
+        if (sfiles->transparentStreaks == 0) {
+            sfiles->outImage->image->data.U16[y][x] = exciseValue;
+        } else {
+            if (streak) {
+                // as a visualization aid don't mask the pixel, just change the intensity
+                sfiles->outImage->image->data.U16[y][x] = imageValue + sfiles->transparentStreaks;
+            } else {
+                sfiles->outImage->image->data.U16[y][x] = exciseValue;
+            }
+        }
     } else {
-        if (streak) {
-            // as a visualization aid don't mask the pixel, just change the intensity
-            sfiles->outImage->image->data.F32[y][x] = imageValue + sfiles->transparentStreaks;
+        float imageValue  = sfiles->inImage->image->data.F32[y][x];
+        if (sfiles->recImage && !isExciseValue(imageValue, sfiles->inImage->exciseValue) ) {
+            sfiles->recImage->image->data.F32[y][x] = imageValue;
+        }
+
+        if (sfiles->transparentStreaks == 0) {
+            sfiles->outImage->image->data.F32[y][x] = exciseValue;
         } else {
-            sfiles->outImage->image->data.F32[y][x] = exciseValue;
+            if (streak) {
+                // as a visualization aid don't mask the pixel, just change the intensity
+                sfiles->outImage->image->data.F32[y][x] = imageValue + sfiles->transparentStreaks;
+            } else {
+                sfiles->outImage->image->data.F32[y][x] = exciseValue;
+            }
         }
     }
@@ -711,11 +797,10 @@
                 sfiles->inMask->image->data.PS_TYPE_IMAGE_MASK_DATA[y][x];
         }
-        sfiles->outMask->image->data.PS_TYPE_IMAGE_MASK_DATA[y][x] =
-            sfiles->inMask->image->data.PS_TYPE_IMAGE_MASK_DATA[y][x] | newMaskValue;
-    }
-}
-
-static void
-exciseNonWarpedPixels(streakFiles *sfiles, psImageMaskType newMaskValue)
+        sfiles->outMask->image->data.PS_TYPE_IMAGE_MASK_DATA[y][x] |= newMaskValue;
+    }
+}
+
+static long
+exciseNonDiffedPixels(streakFiles *sfiles, psImageMaskType newMaskValue)
 {
     int cell_x0 = sfiles->astrom->cell_x0;
@@ -725,4 +810,6 @@
     int numCols = sfiles->inImage->numCols; // for raw images this was calculated from the width of datasec
     int numRows = sfiles->inImage->numRows; // for raw images this was calculated from the height of datasec
+
+    long excisedPixels = 0;
 
 //    printf("%2d x0: %4d y0: %4d xpar: %d ypar: %d\n", sfiles->extnum, cell_x0, cell_y0, xParity, yParity);
@@ -739,5 +826,5 @@
         }
 
-        psU8 *pixels = sfiles->warpedPixels->data.U8[yChip];
+        psU8 *pixels = sfiles->diffedPixels->data.U8[yChip];
 
         if (xParity == 1) {
@@ -746,4 +833,5 @@
                 if (! *pixels ) {
                     excisePixel(sfiles, xCell, yCell, false, newMaskValue);
+                    excisedPixels++;
                 }
             }
@@ -755,17 +843,19 @@
                 if (!*pixels) {
                     excisePixel(sfiles, xCell, yCell, false, newMaskValue);
+                    excisedPixels++;
                 }
             }
         }
     }
+    return excisedPixels;
 }
 
 static bool
-warpedPixel(streakFiles *sfiles, int x, int y)
+diffedPixel(streakFiles *sfiles, int x, int y)
 {
     PixelPos chipCoord;
 
     if (!CHIP_LEVEL_INPUT(sfiles->stage)) {
-        // if we're here on a skycell image by definition this pixel was warped
+        // if we're here on a skycell image by definition this pixel was diffed
         return true;
     }
@@ -780,12 +870,12 @@
     cellToChipInt(&chipCoord.x, &chipCoord.y, sfiles->astrom, x, y);
 
-    if (chipCoord.x < 0 || chipCoord.x >= sfiles->warpedPixels->numCols) {
+    if (chipCoord.x < 0 || chipCoord.x >= sfiles->diffedPixels->numCols) {
         return false;
     }
-    if (chipCoord.y < 0 || chipCoord.y >= sfiles->warpedPixels->numRows) {
+    if (chipCoord.y < 0 || chipCoord.y >= sfiles->diffedPixels->numRows) {
         return false;
     }
 
-    return psImageGet(sfiles->warpedPixels, chipCoord.x, chipCoord.y) ? true : false;
+    return psImageGet(sfiles->diffedPixels, chipCoord.x, chipCoord.y) ? true : false;
 }
 
Index: branches/tap_branches/magic/remove/src/streaksremove.h
===================================================================
--- branches/tap_branches/magic/remove/src/streaksremove.h	(revision 25900)
+++ branches/tap_branches/magic/remove/src/streaksremove.h	(revision 27838)
@@ -33,15 +33,9 @@
     int         yParity;
     double      exciseValue;
+    psString    program;
 } sFile;
 
-
-typedef enum {
-    IPP_STAGE_NONE = 0,
-    IPP_STAGE_RAW,
-    IPP_STAGE_CHIP,
-    IPP_STAGE_WARP,
-    IPP_STAGE_DIFF
-} ippStage;
-
+// used for error messages
+extern char * streaksProgram;
 
 typedef struct {
@@ -64,4 +58,6 @@
     sFile *inSources;
     sFile *outSources;
+    psMetadata *stats;
+    FILE *statsFile;
     psString class_id;
     pmFPAfile  *inAstrom;
@@ -71,8 +67,10 @@
     pmChip  *chip;  // current chip
     pmCell  *cell;  // current cell
-    psImage *warpedPixels;
+    psImage *diffedPixels;
     psVector    *tiles;
     double  transparentStreaks;
     psString    program_name;
+    psU32   maskStreak;
+    psU32   maskMask;
 } streakFiles;
 
@@ -86,9 +84,12 @@
 extern void cellToChipInt(unsigned int *xChip, unsigned int *yChip, strkAstrom *astrom, int xCell, int yCell);
 
-extern bool computeWarpedPixels(streakFiles *sf);
+extern bool computeDiffedPixels(streakFiles *sf);
 extern void streaksExit(psString, int);
 extern ippStage parseStage(psString);
 extern psString pathToDirectory(char *path);
 extern void setStreakFiles( streakFiles *);
+
+extern bool streaksVersionHeaderFull(psMetadata *header);
+extern psString streaksVersionLong(void);
 
 #define CHIP_LEVEL_INPUT(_stage) ((_stage == IPP_STAGE_RAW) || (_stage == IPP_STAGE_CHIP))
Index: branches/tap_branches/magic/remove/src/streaksreplace.c
===================================================================
--- branches/tap_branches/magic/remove/src/streaksreplace.c	(revision 25900)
+++ branches/tap_branches/magic/remove/src/streaksreplace.c	(revision 27838)
@@ -5,4 +5,6 @@
 static void writeImages(streakFiles *sf, bool exciseImageCube);
 static bool readAndCopyToOutput(streakFiles *sf, bool restoreImageCube);
+
+char *streaksProgram = "streaksreplace";
 
 int main(int argc, char *argv[])
Index: branches/tap_branches/magic/remove/src/warpedpixels.c
===================================================================
--- branches/tap_branches/magic/remove/src/warpedpixels.c	(revision 25900)
+++ 	(revision )
@@ -1,330 +1,0 @@
-#include "streaksremove.h"
-#include "pmAstrometryWCS.h"
-
-// #define DEBUG_PRINT 1
-
-static void addTouchedPixels(streakFiles *sf, psString fileName);
-static void nameCorners(psPlane pt[4]);
-static int xLeft(psPlane *pt, int y);
-static int xRight(psPlane *pt, int y);
-
-bool
-computeWarpedPixels(streakFiles *sf)
-{
-    bool status;
-    pmConfig    *config = sf->config;
-
-    psArray *skycells = psMetadataLookupPtr(&status, config->arguments, "SKYCELLS");
-    if (skycells == NULL) {
-        return false;
-    }
-
-    int n = psArrayLength(skycells);
-
-    psRegion     *bounds = pmChipPixels(sf->chip);
-
-    int width  = bounds->x1 - bounds->x0;
-    int height = bounds->y1 - bounds->y0;
-
-    psFree(bounds);
-
-    sf->warpedPixels = psImageAlloc(width, height, PS_TYPE_U8);
-    psImageInit(sf->warpedPixels, 0);
-
-    for (int i=0; i<n; i++) {
-        psString filename = psArrayGet(skycells, i);
-        psString resolved_name = pmConfigConvertFilename(filename, config, false, false);
-
-        addTouchedPixels(sf, resolved_name);
-        psFree(resolved_name);
-    }
-
-    bool writeTouchedPixels = false; // XXX: make this an argument to the program
-    if (writeTouchedPixels) {
-        // write out the warped pixels
-        psMetadata * header = psMetadataAlloc();
-        if (sf->inAstrom->fpa->toSky->type != PS_PROJ_DIS) {
-            pmAstromWriteWCS(header, sf->inAstrom->fpa, sf->chip, 0.001);
-        }
-        psFits *fits = psFitsOpen("warpedpixels.fits", "w");
-
-        psFitsWriteImage(fits, header, sf->warpedPixels, 0, NULL);
-        psFitsClose(fits);
-        psFree(header);
-    }
-    return true;
-}
-
-static void
-addTouchedPixels(streakFiles *sf, psString fileName)
-{
-    // read fits header
-    psFits *fits = psFitsOpen(fileName, "r");
-    if (!fits) {
-        psError(PS_ERR_IO, false, "failed to open skycell file: %s", fileName);
-        streaksExit("", PS_EXIT_DATA_ERROR);
-    }
-    psMetadata *header = psFitsReadHeader(NULL, fits);
-    psFitsClose(fits);
-    if (!header) {
-        psError(PS_ERR_IO, false, "failed to read fixts header from skycell file: %s", fileName);
-        streaksExit("", PS_EXIT_DATA_ERROR);
-    }
-    // set up astrometry for conversion from skycell image to sky
-    pmAstromWCS *wcs = pmAstromWCSfromHeader(header);
-    if (!wcs) {
-        psError(PS_ERR_IO, false, "failed to read astrometry  from skycell file: %s", fileName);
-        streaksExit("", PS_EXIT_DATA_ERROR);
-    }
-    int naxis1 = psMetadataLookupS32(NULL, header, "NAXIS1");
-    int naxis2 = psMetadataLookupS32(NULL, header, "NAXIS2");
-
-    psFree(header);
-
-    /* now set up our wrapper to the chip astrometry to apply to the whole chip */
-    sf->astrom = streakSetAstrometry(sf->astrom, sf->stage, sf->inAstrom->fpa, sf->chip, false, NULL,
-        sf->warpedPixels->numCols, sf->warpedPixels->numRows);
-
-    // convert corners of skycell to sky coordinates
-    // compute overlap with the chip using astrometry
-    psPlane pt[4];
-    pt[0].x = 0;
-    pt[0].y = 0;
-
-    pt[1].x = naxis1 - 1;
-    pt[1].y = 0;
-
-    pt[2].x = naxis1 - 1;
-    pt[2].y = naxis2 - 1;
-
-    pt[3].x = 0;
-    pt[3].y = naxis2 - 1;
-
-    for (int i = 0; i < 4 ; i++ ) {
-        psSphere sky;
-
-        // convert corner of skycell to sky coordinates
-        if (!pmAstromWCStoSky(&sky, wcs, &pt[i])) {
-            psErrorStackPrint(stderr, "failed to convert pt %d of %s to sky coords", i, fileName);
-            streaksExit("", PS_EXIT_DATA_ERROR);
-        }
-        strkPt p;
-        // convert to chip coordinates
-        if (!skyToCell(&p, sf->astrom, sky.r, sky.d)) {
-            psErrorStackPrint(stderr, "failed to convert pt %d of %s to sky coords", i, fileName);
-            streaksExit("", PS_EXIT_DATA_ERROR);
-        }
-        pt[i].x = p.x;
-        pt[i].y = p.y;
-    }
-
-    psFree(wcs);
-
-    // put the corners in the desired order (see comments below)
-
-    nameCorners(pt);
-    psString type;
-    if (pt[0].y < pt[2].y ) {
-        type = "1";
-    } else {
-        type = "2";
-    }
-#ifdef DEBUG_PRINT
-    printf("\nSKYCELL %s Type: %s\n", fileName, type);
-    for (int i=0; i<4; i++) {
-        printf("%f %f\n", pt[i].x, pt[i].y);
-    }
-#endif
-    // Now set the touched pixels
-    int ymin = fmax(0, pt[1].y );
-    int ymax = fmin(pt[3].y + 0.5, sf->warpedPixels->numRows - 1);
-#if (DEBUG_PRINT > 1)
-    printf("\nymin: %d ymax: %d\n", ymin, ymax);
-#endif
-    for (int y = ymin ; y <= ymax; y++) {
-        int xleft  = xLeft(pt, y);
-        int xright = xRight(pt, y);
-
-        if (xleft < 0) {
-            xleft = 0;
-        }
-        if (xleft > sf->warpedPixels->numCols) {
-            continue;
-        }
-        if (xright < 0) {
-            continue;
-        }
-        if (xright >= sf->warpedPixels->numCols) {
-            xright = sf->warpedPixels->numCols - 1;
-        }
-#if (DEBUG_PRINT > 1)
-        printf("  y: %d xleft: %d xright: %d\n", y, xleft, xright);
-#endif
-
-        psU8 *rowPixels = sf->warpedPixels->data.U8[y];
-
-        int n = xright - xleft + 1;
-        if (n > 0) {
-            memset(&rowPixels[xleft], 255, n);
-        }
-    }
-}
-
-// x as a function of y for the line between two points
-// Note: the caller guarentees that the y's of the two points are different
-static double xOfY(psPlane *pI, psPlane *pJ, int y)
-{
-    double m_inv = (double) (pJ->x - pI->x) / (double) (pJ->y - pI->y);
-
-    double x = ((double) y - pI->y) * m_inv + (double) pI->x;
-
-    return x;
-}
-
-static int xLeft(psPlane *pt, int y)
-{
-    double x_d;
-    if (y < pt[0].y) {
-        // left boundary is the line from pt 0 to pt 1
-        x_d = xOfY(&pt[0], &pt[1], y);
-    } else {
-        // left boundary is the line from pt 0 to pt 3
-        x_d = xOfY(&pt[0], &pt[3], y);
-    }
-    return (int) (x_d + 0.5);
-}
-
-static int xRight(psPlane *pt, int y)
-{
-    double x_d;
-    if (y < pt[2].y) {
-        // right boundary is the line from pt 1 to pt 2
-        x_d = xOfY(&pt[1], &pt[2], y);
-
-    } else {
-        // right boundary is the line from pt 2 to pt 3
-        x_d = xOfY(&pt[2], &pt[3], y);
-    }
-    return (int) x_d;
-}
-
-/*
- * To compute the overlap of a quadrilateral with a set of
- * points we transform the 4 corners to image coordinates
- * and name the 4 points of the quad
-            pt 0 is left most (bottom most corner)
-            pt 1 is bottom most (right most corner)
-            pt 2 is right most (top most corner)
-            pt 3 is top most (left most corner)
-
-* Then we can divide the quad into 3 regions A B and C based
-* on the y coordinate
-
-* (in the degenerate case of a rectangle aligned to the axes
-* we only have 1 region)
-* In each of the three regions A, B, and C the test for whether a
-* point is contained in the quad on scanline y is simply
-
-        x >= xleft(y) x < xright(y)
-
-* where xleft(y) and xright(y) are the bounding lines at the given
-* y coordinate.
-
-The following diagrams show some examples. The 4 points are labeled 0 to 3
-The three regions are labeld A B and C, the boundaries between the regions
-ocur at pt0.y and pt2.y
-
-Example 1
-
-                3
-                 C
-         ---------------2           left boundary:  line 0_1 y < pt0.y
-                 B                                  line 0_3 y >= pt0.y
-          0--------------
-                 A                  right boundary: line 1_2 y < pt2.y
-                   1                                line 2_3  y >= pt.y
-**************************
-Example 2
-                3
-              C
-        0----------------
-              B
-      ----------------2
-              A
-            1
-
-**************************
-Example 3
-
-        3       2
-
-            B
-
-        0       1
-
-(region A and C are empty in the case where the points form a rectangle)
-*/
-
-/*
-    Name the corners
-    The following algorithm works for points that form a quadrilateral
-    I think it also works for the situation where 3 points are co-linear
-    and we have a triangle, but that isn't important for our purposes
-
-*/
-
-static void
-nameCorners(psPlane pt[4])
-{
-    psPlane   pt0, pt1, pt2, pt3;
-
-    /* find pt0 the left most (bottom most) corner */
-    int imin = 0;
-    int  min = (int) pt[0].x;
-    int i;
-    for (i=1; i<4; i++) {
-        int x = pt[i].x;
-        if ((x < min) ||
-            ((x == min) && (pt[i].y < pt[imin].y))) {
-            imin = i;
-            min = x;
-        }
-    }
-    // remve pt0 from the list
-    pt0 = pt[imin];
-    for (i=imin; i<3; i++) {
-        pt[i] = pt[i+1];
-    }
-
-    // find the bottom most (right most) corner from the remaining 3
-    imin = 0;
-    min = pt[0].y;
-    for (i=1; i<3; i++) {
-        int y = pt[i].y;
-        if ((y < min) ||
-            ((y == min) && (pt[i].x > pt[imin].x)) ) {
-            imin = i;
-            min = y;
-        }
-    }
-    // remve pt1 from the list
-    pt1 = pt[imin];
-    for (i=imin; i<3; i++) {
-        pt[i] = pt[i+1];
-    }
-
-    // now find the right most (top most) of the remaining 2 points
-    if ((pt[0].x > pt[1].x) ||
-        ((pt[0].x == pt[1].x) && (pt[0].y > pt[1].y))) {
-        pt2 = pt[0];
-        pt3 = pt[1];
-    } else {
-        pt2 = pt[1];
-        pt3 = pt[0];
-    }
-    /* write the outputs */
-    pt[0] = pt0;
-    pt[1] = pt1;
-    pt[2] = pt2;
-    pt[3] = pt3;
-}