Changeset 25082

Timestamp:

Aug 14, 2009, 3:58:07 PM (17 years ago)

Author:

Paul Price

Message:

Fix streaksremove so that output CMF files have PHU called ext.hdr and table called ext.psf. Turned on -Wall -Werror and fixed a host of little problems.

Location:

trunk/magic

Files:

• Makefile.in (modified) (1 diff)
• remove/src/Line.c (modified) (23 diffs)
• remove/src/Makefile.simple (modified) (1 diff)
• remove/src/isdestreaked.c (modified) (2 diffs)
• remove/src/streaksastrom.c (modified) (9 diffs)
• remove/src/streakscompare.c (modified) (4 diffs)
• remove/src/streaksextern.h (modified) (1 diff)
• remove/src/streaksio.c (modified) (8 diffs)
• remove/src/streaksrelease.c (modified) (5 diffs)
• remove/src/streaksremove.c (modified) (12 diffs)
• remove/src/streaksremove.h (modified) (1 diff)
• remove/src/streaksreplace.c (modified) (4 diffs)
• remove/src/streaksutil.c (modified) (1 diff)
• remove/src/warpedpixels.c (modified) (11 diffs)

trunk/magic/Makefile.in

@@ -12 +12 @@
 
 update:
-        # get the source code (replace this with SVN interactions)
+# get the source code (replace this with SVN interactions)
         tar xvzf ~/magic.tgz
         tar xvzf ~/ssa-core-cpp.tgz

trunk/magic/remove/src/Line.c

@@ -36 +36 @@
     // If the current line is not vertical, check to see if the point
     // is inside the two endpoints independent of direction
-    
+
     if (line->begin.x != line->end.x)
     {
@@ -59 +59 @@
     ux = line->end.x - line->begin.x;
     uy = line->end.y - line->begin.y;
-    
+
     vx = x - line->begin.x;
     vy = y - line->begin.y;
-    
+
     double u_u = ux * ux + uy * uy;                         // (u . u)
     double u_v = ux * vx + uy * vy;                         // (u . v)
-    
+
     if (u_v <= 0) return vx * vx + vy * vy;                 // (v . v)
     if (u_u <= u_v)
@@ -73 +73 @@
         return wx * wx + wy * wy;                           // (w . w)
     }
-    
+
     // Compute P(b) is the base of the perpendicular dropped from tuple to
     // the line
-    
+
     b = u_v / u_u;
     px = vx - b * ux;
@@ -109 +109 @@
     double wx  = first->begin.x - second->begin.x;
     double wy  = first->begin.y - second->begin.y;
-    
+
     // Calculate the perpendicular product, dt
-    
+
     double dt = dx1 * dy2 - dy1 * dx2;
-    
+
     // Are the lines parallel?
-    
+
     if (fabs (dt) <= 1e-8)
     {
         // Check to see if they are overlap
-        
+
         if (dx1 * wy - dy1 * wx != 0 || dx2 * wy - dy2 * wx != 0) return 0;
-        
+
         // The line are coplanar, so check to see if they are degenerate
-        
+
         len1 = dx1 * dx1 + dy1 * dy1;
         len2 = dx2 * dx2 + dy2 * dy2;
-        
+
         // First, check to see if both lines are points and if they are
         // distinct tuples
-        
+
         if (len1 == 0 && len2 == 0)
         {
@@ -137 +137 @@
             return 1;
         }
-        
+
         // Check to see if the first line is a point and is inside the
         // second line
-        
+
         if (len1 == 0)
         {
@@ -148 +148 @@
                 return 0;
         }
-        
+
         // Check to see if the second line is a point and is inside the
         // first line
-        
+
         if (len2 == 0)
         {
@@ -159 +159 @@
                 return 0;
         }
-        
+
         // Since both lines are coplanar and have length, search for
         // overlap tuples
-        
+
         w2x = first->end.x - second->begin.x;
         w2y = first->end.y - second->begin.y;
-        
+
         t0 = (dx2 != 0) ? wx  / dx2 : wy  / dy2;
         t1 = (dx2 != 0) ? w2x / dx2 : w2y / dy2;
-        
+
         if (t0 > t1) SwapDouble (&t0, &t1);
         if ((inFirstSegment || inSecondSegment) && (t0 > 1 || t1 < 0))
@@ -176 +176 @@
         t0 = (t0 < 0) ? 0 : t0;             // Clip to min 0
         t1 = (t1 > 1) ? 1 : t1;             // Clip to max 1
-        
+
         // Set the first tuple of intersection
-        
+
         tuple1->x = second->begin.x + t0 * dx2;
         tuple1->y = second->begin.y + t0 * dy2;
-        
+
         // Check to see if the intersection is a single tuple
-        
+
         if (t0 == t1) return 1;
-        
+
         // Otherwise, set the second tuple of intersection as well
-        
+
         tuple2->x = second->begin.x + t1 * dx2;
         tuple2->y = second->begin.y + t1 * dy2;
         return 2;
     }
-    
+
     // The segments are skew and may intersect in a tuple.
     // Get the intersect parameter for the first line
-    
+
     double i1 = (dx2 * wy - dy2 * wx) / dt;
     if (inFirstSegment && (i1 < 0 || i1 > 1)) return 0;
-    
+
     // Get the intersect parameter for the second line
 
@@ -226 +226 @@
     vertices[3].x = 0;       vertices[3].y = numRows;
 
+    clipLine.begin.x = 0;
+    clipLine.begin.y = 0;
+    clipLine.end.x = 0;
+    clipLine.end.y = 0;
+
     for (i = 0; i < 4 && found < 2; ++i)
     {
@@ -237 +242 @@
                 ++found;
             }
-            else if (tuple1.x != clipLine.begin.x || 
+            else if (tuple1.x != clipLine.begin.x ||
                      tuple1.y != clipLine.begin.y)
             {
@@ -245 +250 @@
         }
     }
-    
+
     // If two endpoints are found, clip the line
-    
+
     if (found > 1)
     {
@@ -283 +288 @@
     vertices[3].x = minX; vertices[3].y = maxY;
 
+    clipLine.begin.x = 0;
+    clipLine.begin.y = 0;
+    clipLine.end.x = 0;
+    clipLine.end.y = 0;
+
     for (i = 0; i < 4 && found < 2; ++i)
     {
@@ -294 +304 @@
                 ++found;
             }
-            else if (tuple1.x != clipLine.begin.x || 
+            else if (tuple1.x != clipLine.begin.x ||
                      tuple1.y != clipLine.begin.y)
             {
@@ -302 +312 @@
         }
     }
-    
+
     // If two endpoints are found, clip the line
-    
+
     if (found > 1)
     {
@@ -360 +370 @@
     @param[out] pixels list of PixelPos pointers corresponding
                        based on the line settings
-    @param[in] line Line to map to pixels   
+    @param[in] line Line to map to pixels
     @param[in] numCols maximum X (columns) for the line
     @param[in] numRows maximum Y (rows) for the line            */
@@ -367 +377 @@
 {
     Line offsetLine;
-    PixelPos *pixel;
     double slope, xOffset, yOffset, xMid, yMid;
     int x, y, xBegin = numCols, yBegin = numRows, xEnd = 0, yEnd = 0;
 
     // Extract the endpoints
-    
+
     double x1 = line->begin.x;
     double y1 = line->begin.y;
     double x2 = line->end.x;
     double y2 = line->end.y;
-    
+
     // Determine the width and height of the line
-    
+
     double dx = x2 - x1;
     double dy = y2 - y1;
@@ -403 +412 @@
 
     // Step point by point based on the dominate axis
-    
+
     if (fabs (dx) > fabs (dy))
     {
@@ -411 +420 @@
         // If line is back to front, reorder endpoints and recompute
         // the line width and height
-        
+
         if (x1 > x2)
         {
@@ -421 +430 @@
 
         // Compute the slope of the line
-        
+
         slope = dy / dx;
-        
+
         // Compute the x and y offsets for the line width extent
 
-        if (xBegin > xEnd) 
+        if (xBegin > xEnd)
             SwapInt (&xBegin, &xEnd);
         else
@@ -451 +460 @@
     }
     else
-    {       
+    {
         // -------------------
         // vertical(ish) lines
@@ -457 +466 @@
         // If line is back to front, reorder endpoints and recompute
         // the line width and height
-        
+
         if (y1 > y2)
         {
@@ -465 +474 @@
             dy = y2 - y1;
         }
-        
+
         // Compute the slope of the line
-        
+
         slope = dx / dy;
-        
+
         // Compute the x and y offsets for the line width extent
 
@@ -481 +490 @@
         xMid = x1 + slope * (yBegin - y1);
         xOffset = fabs (halfWidth * dr / dy);
-        
+
         for (y = yBegin; y != yEnd; ++y)
         {

trunk/magic/remove/src/Makefile.simple

@@ -36 +36 @@
 
 STREAKSFLAGS=-DSTREAKS_COMPRESS_OUTPUT=1
-OPTFLAGS= -g -O2
+OPTFLAGS= -g -O2 -Wall -Werror
 #OPTFLAGS= -g
-CFLAGS=`psmodules-config --cflags` -std=gnu99 ${OPTFLAGS} ${STREAKSFLAGS}
-#CFLAGS=`psmodules-config --cflags` -std=gnu99 ${OPTFLAGS}
-LDFLAGS=`psmodules-config --libs`
+CFLAGS=`psmodules-config --cflags` `pslib-config --cflags` -std=gnu99 ${OPTFLAGS} ${STREAKSFLAGS}
+#CFLAGS=`psmodules-config --cflags` `pslib-config --cflags` -std=gnu99 ${OPTFLAGS}
+LDFLAGS=`psmodules-config --libs` `pslib-config --libs`
 
 PROGRAMS= streaksremove streaksreplace streakscompare streaksrelease isdestreaked

trunk/magic/remove/src/isdestreaked.c

@@ -18 +18 @@
     psFits *fits = psFitsOpen(fileName, "r");
     if (!fits) {
-        psError(PS_ERR_UNKNOWN, false, "failed to open fits file: %s\n", fileName);
-        psErrorStackPrint(stderr, "");
+        psErrorStackPrint(stderr, "failed to open fits file: %s\n", fileName);
         exit(PS_EXIT_DATA_ERROR);
     }
@@ -25 +24 @@
     psMetadata *header = psFitsReadHeader(NULL, fits);
     if (!header) {
-        psError(PS_ERR_IO, false, "failed to read header for: %s", fileName);
-        psErrorStackPrint(stderr, "");
+        psErrorStackPrint(stderr, "failed to read header for: %s", fileName);
         exit(PS_EXIT_DATA_ERROR);
     }

trunk/magic/remove/src/streaksastrom.c

@@ -93 +93 @@
     // since the transform already accounts for the chip parity we just use the cell parity
     astrom->xParity = (double) xParityCell;
-    astrom->yParity = (double) yParityCell; 
+    astrom->yParity = (double) yParityCell;
 
     astrom->numCols = numCols;
@@ -135 +135 @@
 
 void
-cellToChipInt(int *xChip, int *yChip, strkAstrom *astrom, int xCell, int yCell)
+cellToChipInt(unsigned int *xChip, unsigned int *yChip, strkAstrom *astrom, int xCell, int yCell)
 {
     if (astrom->xParity > 0) {
@@ -148 +148 @@
     }
 }
- 
+
 bool
 SkyToLocal(strkPt *outPt, strkAstrom *astrom, double ra, double dec)
@@ -156 +156 @@
 
     pmFPA *fpa   = (pmFPA *) astrom->fpa;
-    pmChip *chip = (pmChip *) astrom->chip;
 
     // find the RA,DEC coords of the 0,0 pixel for this chip:
@@ -183 +182 @@
 
     pmFPA *fpa   = (pmFPA *) astrom->fpa;
-    pmChip *chip = (pmChip *) astrom->chip;
 
     // find the RA,DEC coords of the 0,0 pixel for this chip:
@@ -251 +249 @@
     pmChip *chip = (pmChip *) astrom->chip;
     pmAstromObj *pt = (pmAstromObj *) astrom->pt;
-    
+
     pt->sky->r = ra;
     pt->sky->d = dec;
@@ -276 +274 @@
     pmChip *chip = (pmChip *) astrom->chip;
     pmAstromObj *pt = (pmAstromObj *) astrom->pt;
-    
+
     cellToChip(&pt->chip->x, &pt->chip->y, astrom, x, y);
 
@@ -294 +292 @@
 }
 
- 
+
 static bool
 readAstrometry(streakFiles *sf)
@@ -362 +360 @@
             break;
         }
-    } 
+    }
     if (!sf->chip) {
         psError(PS_ERR_UNKNOWN, true, "Failed to find chip with data.");

trunk/magic/remove/src/streakscompare.c

@@ -5 +5 @@
 int main(int argc, char *argv[])
 {
-    long i;
     bool status;
 
@@ -35 +34 @@
     int numErrors = 0;
     for (int component = 0; component < ncomponents; component++) {
-        if (component && !(psFitsMoveExtNum(file1->fits, 1, true) 
+        if (component && !(psFitsMoveExtNum(file1->fits, 1, true)
                            && psFitsMoveExtNum(file2->fits, 1, true))) {
             streaksExit("failed to advance to next extesion\n", PS_EXIT_DATA_ERROR);
@@ -45 +44 @@
         psImage *image2 = file2->image;
 
-        // TODO: do more sanity checking. For example check that extname's  (if any) match 
+        // TODO: do more sanity checking. For example check that extname's  (if any) match
         // check for matching image cubes
         if (image1 && image2) {
@@ -68 +67 @@
                         error = ! isnan(value1) && isnan(value2);
                     }
-                    
+
                     if (error) {
                         numErrors++;

trunk/magic/remove/src/streaksextern.h

@@ -54 +54 @@
             pixels
 */
-    
+
 extern StreakPixels *streak_on_component(Streaks *streaks, strkAstrom *astrom,
                                          int numCols, int numRows);

trunk/magic/remove/src/streaksio.c

@@ -12 +12 @@
 // Assumptions about the file structure of non-raw files
 // The 'image' for each file (image, mask weight) is contained in the first
-// extension. 
+// extension.
 
 
@@ -625 +625 @@
         streaksExit("", PS_EXIT_DATA_ERROR);
     }
- 
+
     bool status;
     in->extname = psMetadataLookupStr(&status, in->header, "EXTNAME");
@@ -654 +654 @@
             streaksExit("", PS_EXIT_DATA_ERROR);
         }
-        psImage *image = (psImage *) (in->imagecube->data[0]);
     }
     setDataExtent(stage, in, (stage == IPP_STAGE_RAW) && !isMask);
@@ -925 +924 @@
 replicateOutputs(streakFiles *sfiles)
 {
-    bool status = false;
-
     if (!replicate(sfiles->outImage, sfiles->inImage)) {
         psError(PM_ERR_SYS, false, "failed to replicate outImage.");
@@ -973 +970 @@
 
         if (!nebSwap(server, in->name, out->name)) {
-            psError(PM_ERR_SYS, true, "failed to swap files for: %s.",
+            psError(PM_ERR_SYS, true, "failed to swap files for %s: %s.",
                 in->name, nebErr(server));
             return false;
@@ -1097 +1094 @@
 setMaskedToNAN(streakFiles *sfiles, psU32 maskMask, bool printCounts)
 {
-        int maskedPixels = 0;
-        int nandPixels = 0;
-        int nandWeights = 0;
+        long maskedPixels = 0;
+        long nandPixels = 0;
+        long nandWeights = 0;
 
         psImage *image = sfiles->outImage->image;
@@ -1120 +1117 @@
                 psU32 maskVal;
                 if (sfiles->stage == IPP_STAGE_RAW) {
-                    int xChip, yChip;
+                    unsigned int xChip, yChip;
                     cellToChipInt(&xChip, &yChip, sfiles->astrom, x, y);
                     maskVal = psImageGet(mask, xChip, yChip);
@@ -1144 +1141 @@
         if (printCounts) {
             psLogMsg(sfiles->program_name, PS_LOG_INFO, "time to NAN mask pixels: %f\n", psTimerClear("NAN_MASKED"));
-            int totalPixels = image->numRows * image->numCols;
+            long totalPixels = image->numRows * image->numCols;
             psLogMsg(sfiles->program_name, PS_LOG_INFO, "pixels:        %10ld\n", totalPixels);
             psLogMsg(sfiles->program_name, PS_LOG_INFO, "masked pixels: %10ld %4.2f%%\n", maskedPixels, 100. * maskedPixels / totalPixels);

trunk/magic/remove/src/streaksrelease.c

@@ -3 +3 @@
 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);
 
@@ -12 +8 @@
 main(int argc, char *argv[])
 {
-    bool status;
-
     psLibInit(NULL);
     psTimerStart("STREAKSREMOVE");
@@ -25 +19 @@
     // 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 
+    psU32 maskMask = 0;             // value looked up for MASK.STREAK
 
     // Does true work here?
@@ -216 +210 @@
         // image data directly from psFits
         readImage(sf->inImage, sf->extnum, sf->stage, false);
-        
+
         // astrom struct is only needed for raw stage, and only the chip to cell parameters are used
         sf->astrom = streakSetAstrometry(sf->astrom, sf->stage, NULL, NULL, false,
@@ -304 +298 @@
             initValue = NAN;
         } else {
-            // otherwise write it to the output 
+            // otherwise write it to the output
             writeImageCube(sf->outImage, sf->inImage->imagecube, extname, sf->extnum);
             initValue = 0;

trunk/magic/remove/src/streaksremove.c

@@ -35 +35 @@
     // 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 
+    psU32 maskMask = 0;             // value looked up for MASK.STREAK
 
     psString streaksFileName = psMetadataLookupStr(NULL, config->arguments, "STREAKS");
@@ -42 +42 @@
     Streaks *streaks = readStreaksFile(streaksFileName);
     if (!streaks) {
-        psError(PS_ERR_UNKNOWN, "failed to read streaks file: %s", streaksFileName);
+        psError(PS_ERR_UNKNOWN, false, "failed to read streaks file: %s", streaksFileName);
         streaksExit("", PS_EXIT_PROG_ERROR);
     }
@@ -83 +83 @@
         psLogMsg("streaksremove", PS_LOG_INFO, "time to compute warped pixels: %f\n", cwp_t);
     }
-    
+
     if (sfiles->stage == IPP_STAGE_RAW) {
         // Except for raw stage, all of our (GPC1) files have one image extension.
@@ -97 +97 @@
     }
 
-    int totalPixels = 0;
-    int totalStreakPixels = 0;
+    long totalPixels = 0;
+    long totalStreakPixels = 0;
 
     // Iterate through each component of the input (except for raw images there is only one)
@@ -132 +132 @@
                                                         sfiles->inImage->numRows);
             psLogMsg("streaksremove", PS_LOG_INFO, "time to get streak pixels: %f\n", psTimerClear("GET_STREAK_PIXELS"));
-            
+
             // if this extension contained an image, excise the streaked pixels.
             // otherwise it contained an image cube (video cell) which is handled in the if block
@@ -172 +172 @@
                 }
 
-            } else { 
+            } else {
                 // this component contains an image cube
                 // For now excise it completely
@@ -215 +215 @@
 
     if (!replicateOutputs(sfiles)) {
-        psError(PS_ERR_UNKNOWN, false, "failed to replicate output files");
+        psErrorStackPrint(stderr, "failed to replicate output files");
         deleteTemps(sfiles);
-        psErrorStackPrint(stderr, "");
         exit(PS_EXIT_UNKNOWN_ERROR);
     }
@@ -247 +246 @@
     pmConceptsDone();
     pmModelClassCleanup();
-    streaksNebulousCleanup(); 
+    streaksNebulousCleanup();
     pmConfigDone();
     psLogMsg("streaksremove", PS_LOG_INFO, "time to run streaksremove: %f\n", psTimerClear("STREAKSREMOVE"));
@@ -356 +355 @@
             true);
     }
-    
+
     if ((argnum = psArgumentGet(argc, argv, "-keepnonwarped"))) {
         psArgumentRemove(argnum, &argc, argv);
@@ -374 +373 @@
         psArgumentRemove(argnum, &argc, argv);
     }
-        
+
     // if skycells are not provided then we have to execise all pixels  unless -keepnonwarped
     pmConfigFileSetsMD(config->arguments, &argc, argv, "SKYCELLS", "-skycell", "-skycelllist");
@@ -628 +627 @@
             initValue = NAN;
         } else {
-            // otherwise write it to the output 
+            // otherwise write it to the output
             writeImageCube(sf->outImage, sf->inImage->imagecube, extname, sf->extnum);
             initValue = 0;
@@ -799 +798 @@
     sFile *out = sfiles->outSources;
 
-    in->header = psFitsReadHeader(NULL, in->fits);
-    if (!in->header) {
-        psError(PS_ERR_IO, false, "failed to read header from %s", in->resolved_name);
-        streaksExit("", PS_EXIT_DATA_ERROR);
-    }
-
-    bool status;
-    psString extname = psMetadataLookupStr(&status, in->header, "EXTNAME");
-    if (!extname) {
-        psError(PS_ERR_IO, false, "failed to find extname in header of %s", in->resolved_name);
-        streaksExit("", PS_EXIT_DATA_ERROR);
-    }
-
-    psArray *inTable = psFitsReadTable(in->fits);
-    if (!inTable->n) {
-        psError(PS_ERR_IO, false, "table in %s is empty", in->resolved_name);
-        streaksExit("", PS_EXIT_DATA_ERROR);
-    }
-
-    psArray *outTable = psArrayAllocEmpty(inTable->n);
-    int j = 0;
-    int numCensored = 0;
-    for (int i = 0 ; i < inTable->n; i++) {
-        psMetadata *row = inTable->data[i];
-
-        psF32 x = psMetadataLookupF32 (&status, row, "X_PSF");
-        psF32 y = psMetadataLookupF32 (&status, row, "Y_PSF");
-        
-        psU32 mask = psImageGet(maskImage, x, y);
-
-        // Key the source if the center pixel is not masked with maskStreak
-        if (! (mask & maskStreak) ) {
-            psArraySet(outTable, j++, row);
+
+    // Primary header, should be "something.hdr"
+    {
+        psMetadata *header = psFitsReadHeader(NULL, in->fits);
+        if (!header) {
+            psError(PS_ERR_IO, false, "failed to read header from %s", in->resolved_name);
+            streaksExit("", PS_EXIT_DATA_ERROR);
+        }
+
+        bool status;
+        psString extname = psMetadataLookupStr(&status, header, "EXTNAME");
+        if (!extname) {
+            psError(PS_ERR_IO, false, "failed to find extname in header of %s", in->resolved_name);
+            streaksExit("", PS_EXIT_DATA_ERROR);
+        }
+        addDestreakKeyword(header);
+
+        if (!psFitsWriteBlank(in->fits, header, extname)) {
+            psError(PS_ERR_IO, false, "failed to write blank in header of %s", in->resolved_name);
+            streaksExit("", PS_EXIT_DATA_ERROR);
+        }
+        psFree(header);
+    }
+
+    // Extension with PSF fits, should be "something.psf"
+    {
+        psMetadata *header = psFitsReadHeader(NULL, in->fits);
+        if (!header) {
+            psError(PS_ERR_IO, false, "failed to read header from %s", in->resolved_name);
+            streaksExit("", PS_EXIT_DATA_ERROR);
+        }
+        psString extname = psMetadataLookupStr(NULL, header, "EXTNAME");
+        if (!extname) {
+            psError(PS_ERR_IO, false, "failed to find extname in header of %s", in->resolved_name);
+            streaksExit("", PS_EXIT_DATA_ERROR);
+        }
+
+        psArray *inTable = psFitsReadTable(in->fits);
+        if (!inTable->n) {
+            psError(PS_ERR_IO, false, "table in %s is empty", in->resolved_name);
+            streaksExit("", PS_EXIT_DATA_ERROR);
+        }
+
+        psArray *outTable = psArrayAllocEmpty(inTable->n);
+        int j = 0;
+        int numCensored = 0;
+        for (int i = 0 ; i < inTable->n; i++) {
+            psMetadata *row = inTable->data[i];
+
+            psF32 x = psMetadataLookupF32(NULL, row, "X_PSF");
+            psF32 y = psMetadataLookupF32(NULL, row, "Y_PSF");
+
+            psU32 mask = psImageGet(maskImage, x, y);
+
+            // Key the source if the center pixel is not masked with maskStreak
+            if (!(mask & maskStreak) ) {
+                psArraySet(outTable, j++, row);
+            } else {
+                numCensored++;
+            }
+        }
+
+        // get rid of unused elements (don't know if this is necessary)
+        psArrayRealloc(outTable, j);
+
+        addDestreakKeyword(header);
+        if (psArrayLength(outTable) > 0) {
+            printf("Censored %d sources\n", numCensored);
+            if (! psFitsWriteTable(out->fits, header, outTable, extname)) {
+                psError(PS_ERR_IO, false, "failed to write table to %s", out->resolved_name);
+                streaksExit("", PS_EXIT_DATA_ERROR);
+            }
         } else {
-            numCensored++;
-        }
-    }
-
-    // get rid of unused elements (don't know if this is necessary)
-    psArrayRealloc(outTable, j);
-
-    addDestreakKeyword(in->header);
-    if (psArrayLength(outTable) > 0) {
-        printf("Censored %d sources\n", numCensored);
-        if (! psFitsWriteTable(out->fits, in->header, outTable, extname)) {
-            psError(PS_ERR_IO, false, "failed to write table to %s", out->resolved_name);
-            streaksExit("", PS_EXIT_DATA_ERROR);
-        }
-    } else {
-        printf("Censored ALL %d sources\n", numCensored);
-        if (! psFitsWriteTableEmpty(out->fits, in->header, inTable->data[0], extname)) {
-            psError(PS_ERR_IO, false, "failed to write empty table to %s", out->resolved_name);
-            streaksExit("", PS_EXIT_DATA_ERROR);
-        }
-    }
+            printf("Censored ALL %d sources\n", numCensored);
+            if (! psFitsWriteTableEmpty(out->fits, header, inTable->data[0], extname)) {
+                psError(PS_ERR_IO, false, "failed to write empty table to %s", out->resolved_name);
+                streaksExit("", PS_EXIT_DATA_ERROR);
+            }
+        }
+        psFree(header);
+        psFree(outTable);
+        psFree(inTable);
+    }
+
+    // XXX Will need to update to handle extension with extended sources, etc.
 
     if (!psFitsClose(out->fits)) {

trunk/magic/remove/src/streaksremove.h

@@ -84 +84 @@
 // can't declare this in streaksastrom due to header file ordering
 extern void cellToChip(double *xChip, double *yChip, strkAstrom *astrom, double xCell, double yCell);
-extern void cellToChipInt(int *xChip, int *yChip, strkAstrom *astrom, int xCell, int yCell);
+extern void cellToChipInt(unsigned int *xChip, unsigned int *yChip, strkAstrom *astrom, int xCell, int yCell);
 
 extern bool computeWarpedPixels(streakFiles *sf);

trunk/magic/remove/src/streaksreplace.c

@@ -8 +8 @@
 int main(int argc, char *argv[])
 {
-    long i;
     bool status;
-    StreakPixels *pixels;
-    PixelPos *pixelPos;
 
     psLibInit(NULL);
@@ -92 +89 @@
 
     if (!replicateOutputs(sfiles)) {
-        psError(PS_ERR_UNKNOWN, false, "failed to replicate output files");
-        psErrorStackPrint(stderr, "");
+        psErrorStackPrint(stderr, "failed to replicate output files");
         exit(PS_EXIT_UNKNOWN_ERROR);
     }
@@ -172 +168 @@
             true);
     }
-    
+
     bool gotMask = false;
     if ((argnum = psArgumentGet(argc, argv, "-mask"))) {
@@ -353 +349 @@
         // we have an image cube
         double initValue;
-        // otherwise write it to the output 
+        // otherwise write it to the output
         writeImageCube(sf->outImage, sf->recImage->imagecube, extname, sf->extnum);
 

trunk/magic/remove/src/streaksutil.c

@@ -44 +44 @@
 // we just bail out
 void streaksExit(psString str, int exitCode) {
-    psErrorStackPrint(stderr, str);
+    psErrorStackPrint(stderr, "%s", str);
     if (ourStreakFiles) {
         deleteTemps(ourStreakFiles);

trunk/magic/remove/src/warpedpixels.c

@@ -23 +23 @@
 
     psRegion     *bounds = pmChipPixels(sf->chip);
-    
+
     int width  = bounds->x1 - bounds->x0;
     int height = bounds->y1 - bounds->y0;
@@ -83 +83 @@
 
     /* 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->astrom = streakSetAstrometry(sf->astrom, sf->stage, sf->inAstrom->fpa, sf->chip, false, NULL,
         sf->warpedPixels->numCols, sf->warpedPixels->numRows);
 
@@ -106 +106 @@
         // convert corner of skycell to sky coordinates
         if (!pmAstromWCStoSky(&sky, wcs, &pt[i])) {
-            psError(PS_ERR_IO, false, "failed to convert pt %d of %s to sky coords: %s", fileName);
+            psErrorStackPrint(stderr, "failed to convert pt %d of %s to sky coords", i, fileName);
             streaksExit("", PS_EXIT_DATA_ERROR);
         }
@@ -112 +112 @@
         // convert to chip coordinates
         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);
+            psErrorStackPrint(stderr, "failed to convert pt %d of %s to sky coords", i, fileName);
             streaksExit("", PS_EXIT_DATA_ERROR);
         }
@@ -171 +171 @@
 }
 
-// x as a function of y for the line between two points 
+// 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)
@@ -212 +212 @@
  * 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 
+ * and name the 4 points of the quad
             pt 0 is left most (bottom most corner)
             pt 1 is bottom most (right most corner)
@@ -238 +238 @@
 
                 3
-                 C       
+                 C
          ---------------2           left boundary:  line 0_1 y < pt0.y
                  B                                  line 0_3 y >= pt0.y
@@ -249 +249 @@
               C
         0----------------
-              B     
+              B
       ----------------2
               A
@@ -258 +258 @@
 
         3       2
-            
+
             B
 
@@ -265 +265 @@
 (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 
+    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
@@ -313 +313 @@
         pt[i] = pt[i+1];
     }
-            
+
     // now find the right most (top most) of the remaining 2 points
     if ((pt[0].x > pt[1].x) ||